The wind-driven halo in high-contrast images I: analysis from the focal plane images of SPHERE

Context. The wind driven halo is a feature observed within the images delivered by the latest generation of ground-based instruments equipped with an extreme adaptive optics system and a coronagraphic device, such as SPHERE at the VLT. This signature appears when the atmospheric turbulence conditions are varying faster than the adaptive optics loop can correct. The wind driven halo shows as a radial extension of the point spread function along a distinct direction (sometimes referred to as the butterfly pattern). When present, it significantly limits the contrast capabilities of the instrument and prevents the extraction of signals at close separation or extended signals such as circumstellar disks. This limitation is consequential because it contaminates the data a substantial fraction of the time: about 30% of the data produced by the VLT/SPHERE instrument are affected by the wind driven halo.Aims. This paper reviews the causes of the wind driven halo and presents a method to analyze its contribution directly from the scientific images. Its effect on the raw contrast and on the final contrast after post-processing is demonstrated.Methods. We used simulations and on-sky SPHERE data to verify that the parameters extracted with our method are capable of describing the wind driven halo present in the images. We studied the temporal, spatial and spectral variation of these parameters to point out its deleterious effect on the final contrast.Results. The data driven analysis we propose does provide information to accurately describe the wind driven halo contribution in the images. This analysis justifies why this is a fundamental limitation to the final contrast performance reached.Conclusions. With the established procedure, we will analyze a large sample of data delivered by SPHERE in order to propose, in the future, post-processing techniques tailored to remove the wind driven halo.


Introduction
Thanks to the latest generation of instruments dedicated to exoplanet and circumstellar disk imaging, the last five years have witnessed 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 are orbiting at a few astronomical units (au) from their host star, itself located at a few tens of parsecs from Earth, it is required to reach a contrast better than 10 −5 at an angular separation of 500 milliarcseconds (mas) from the star. By equipping 8-m class telescopes with dedicated instruments combining extreme adaptive optics systems (AO) 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 showing up in the focal plane images that were not visible with the first generation of HCI instruments such as Gem-ini/NICI (Artigau et al. 2008), VLT/NaCo (Rousset et al. 2003) or Keck/NIRC2 (McLean & Chaffee 2000).
The scientific region of interest is the close vicinity of the star (below 500 mas) where the detection of exoplanets is crucial to reject one or the other planet formation scenario as most giant planets are expected to be found in this region (Chauvin 2018;Nielsen et al. 2019), and where circumstellar disks are sometimes expected from the host stars infrared excess. With the latest generation of HCI instruments, the main limitations that particularly affect those inner regions by provoking leakages of the starlight from the focal plane mask element of the coronagraph are (1) the quasi-statics non-common 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 whereas not present in 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 Article number, page 1 of 17 arXiv:2003.05794v1 [astro-ph.IM] 12 Mar 2020 A&A proofs: manuscript no. aanda 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,  due to the atmospheric turbulence evolution being 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 separation 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) present a review of the contrast limitations observed in the VLT/SPHERE images and Mouillet et al. (2018) present a review of the impact of the AO performance on the SPHERE images.
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 AOcorrection, and therefore the command sent to the deformable mirror (DM) is relative to the previous correction. As there is some time delay between the WFS measurement and the setting up of the DM, if the atmospheric turbulence has varied significantly between the measure and the applied correction, the AO servolag error becomes consequential. For a fixed AO delay, the AO servolag error therefore depends on the turbulence coherence time τ 0 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 seconds 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 −3 . In the context of HCI, when using a coronagraph to suppress the coherent peak of the starlight, a raw contrast of about 10 −3 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 Madurowicz et al. 2019): the more the amplitude error and the delayed phase error are correlated, the less the asymmetry. Figure 1 shows this wind driven halo contribution in a coronagraphic image from the VLT/SPHERE-IRDIS instrument (Dohlen et al. 2008), in the case of a simulation (left, infinite exposure with a perfect coronagraph) and an on-sky image (right).
Since a significant part of the images obtained with SPHERE are affected by the WDH, ultimately, we would like to reconstruct this effect to apply correction to existing data acquired during the 5 years of SPHERE operations. For 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 as 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 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 artefacts due to the deformable mirror manufacturing technique. and metrics to analyze its contribution directly from the focal plane images (Sect. 3). We then apply this procedure on on-sky SPHERE images to highlight the impact of the WDH on the contrast after post-processing, by studying its typical spatial, temporal and spectral variations (Sect. 4). From these analysis, 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.

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) and we specify how it affects the point spread function (Sect. 2.3) and finally how it affects the raw contrast in the specific case of coronagraphic imaging (Sect. 2.4).

The adaptive optics temporal lag
A classical on-axis single conjugated AO system is composed of three main components: (i) the wavefront sensor (WFS) analysing the incoming phase distortion, (ii) the real-time computer (RTC) calculating, from the WFS measurement, the command to be sent to the phase corrector and (iii) the deformable mirror (DM) correcting for the phase distortion. The adaptive optics system of the VLT/SPHERE instrument, SAXO, is constituted of a 40×40 sub-aperture spatially-filtered Shack-Hartmann WFS  using an EMCCD detector ; the RTC is the ESO-provided SPARTA architecture (Suárez Valles et al. 2012;Petit et al. 2014); and the correction is in two stages thanks to 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: 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 under use). For SPHERE, T read = 725 µs, including all the necessary operations needed 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, defining the maximum frame rate, is therefore 1/T read = 1/725 µs = 1380 Hz. 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 (it continues in parallel to the readout of the WFS) and finishes 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 the DM takes 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 under use. For SPHERE, T DM is very small compared to 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 follow: 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-lab as t RTC = 80 µs; t DAC is approximated by T int /2 (when assuming that, at first order, 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 fitted empirically on SPHERE, by measuring the closed loop transfer function, to be such as t DM + t com = 35 µs.
For SAXO, the main contributors are therefore the integration time T int , 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 running at 1380 Hz. In practice, according to the latest tests performed on SPHERE (in December 2018), the measurement of frame number n is mainly affected by the command number n − 2 (by 84%) and less affected by 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 1 , which affects the starlight distribution in focal plane. The trade off between the AO loop delay and the temporal evolution of the atmospheric turbulence is the critical parameter provoking 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.

The atmospheric turbulence temporal variation
At a given instant, the atmospheric turbulence state can be represented as a 2-dimensional phase screen reaching the telescope aperture. This phase screen can be described by its power spectral density (PSD) along models such as the widely used stratified von Karman model (Conan 2000). This model is parametrized by the Fried parameter (r 0 , the typical spatial extension of the turbulence cells) and the outer scale (L 0 , the largest size of the turbulence cells). During the observation sequence, this phase screen evolves in two fashions: by translation (the flow) and by the evolution of the turbulent cells shape and spatial distribution (the boiling). In a multilayer description, the flow is associated with the variations of the wind speed and direction of each layer and therefore mainly affects low spatial frequencies variations. The boiling is associated to 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 the studies Guesalaga et al. 2014;Poyneer & Macintosh 2006;Schöck & Spillar 2000, for three different sites, Cerro Pachon, Mauna Kea and Albuquerque). According to the number of subapertures of the WFS per telescope diameter (20 cm sampling for SAXO) and the AOloop delay (1.56 ms for SPHERE-SAXO), the effect of boiling can be ignored for SPHERE-like instruments. In such case, we can work under the frozen flow assumption 2 (the so-called Taylor hypothesis, Taylor 1938), 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, that is therefore not linked to the wind driven halo.
Under the frozen flow hypothesis, the temporal variation of the atmospheric turbulence is parametrized by the turbulence coherence time, τ 0 , analytically 3 defined as (Roddier 1981;Roddier et al. 1982;Hardy 1998): where r 0 is the Fried parameter (dependent upon the wavelength of observation and the zenith angle) and v e f f is the effective wind velocity, itself defined as v e f f = ∞ C 2 , with v(h) the wind velocity profile with altitude h, and C 2 n (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 2 . If τ 0 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 of 1 rad 2 , that is  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 in about 2.2 frames delay, starting in the middle of the first frame (T int /2) until the DM is effectively in the corresponding shape. to say, its Strehl ratio decreased by 63% due to the servolag error (under the small 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 three years of MASS-DIMM (Multi-Aperture Scintillation Sensor -Differential Image Motion Monitor, Kornilov et al. 2007) measurements is presented in Fig. 3 and shows a steep curve at short τ 0 (below the median of 4.5 ms). It is not possible to obtain a direct trade-off 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 τ 0 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. 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 on the C 2 n (h) profile, extracted from the Stereo-SCIDAR (SCIntillation Detection And Ranging, Vernin & Rod-dier 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 on 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) and whose wind speed can go up to 50 m/s, depending on the latitude and the season 4 . 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 e f f max(v ground , 0.4v jet−stream ), showing the impact of the jet stream and how it affects the astronomical data (Masciadri et al. 2013). In order to highlight the predominance of the jet stream layer to the effective wind speed, and therefore to 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 e f f , extracted from the Stereo-SCIDAR data. The median contribution of the ground layer to v e f f is about 10% while the median contribution of the jet stream layer to v e f f 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 e f f . By correlating the observed WDH direction within high-contrast images from GPI, Madurowicz et al. (2018) also show that the jet stream layer is indeed the main responsible for the apparition of the WDH in HCI data.
From 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, 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 of up to typically 5 m/s over an hour while the change in wind direction is likely to remain below 10 • .
As reminded in the introduction, the asymmetry of the WDH is due to the scintillation. For a given AO delay and turbu- lence coherence time, the asymmetry of the wind driven halo increases with the scintillation . 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 2 n (h), the diameter and structure of the telescope aperture D tel , the wavelength of observation 5 λ, 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 σ 2 sphere. In addition, they observed that the power of the scintillation is quite stable over 1 h timescales at Paranal observatory 7 , directly related to the jet stream speed variations, which are also shown in Fig. 6 (top). This indicates again the importance of the jet stream layer on the WDH signature.
In the following, all these aspects are taken into account to simulate SPHERE-like images.

The AO servolag error consequences in the images
The spatial variance of the AO residual phase in the pupil due to the AO servolag error varies as: For a single conjugated AO system, the power spectral density (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;: with h T = 2/(λk 2 ) 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). Note that 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 respectively). The residual phase PSD due to the AO servolag error shows a low spatial frequency structure along the effective wind direction (white arrow on Fig. 7, right) decreasing radially (until the AO correction radius) and shows no power in the direction perpendicular to the effective wind direction. In addition, due to 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).
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 τ 0 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 provoked by other typical AO errors (aliasing, chromaticity, anisoplanetism and noise propagation, but excluding NCPA or other exogenous errors).

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 7 According to the study led in Kornilov et al. (2012), the temporal behavior of the scintillation is similar among the 11 surveyed sites.
Effec%ve wind direc%on Fig. 7. Power spectral densities 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 colorbar is normalized to the maximum of the right image.
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 non-coronagraphic 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 along one given direction, varying only the wind velocity, under the median seeing conditions at Paranal observatory. To assess the impact of the WDH on the raw contrast, Fig. 8 shows the radial profiles of simulated coronagraphic images for various τ 0 , along the wind direction (solid lines, where the raw contrast is highly impacted) and along its perpendicular direction (dashed lines, where the raw contrast is less impacted), under median seeing condition (r 0 = 13.8 cm, according to MASS-DIMM measurements) and median airmass (a =1.15).
We compared these simulations with the raw contrast of a SPHERE image taken under very good observing conditions (median seeing, τ 0 > 9 ms), in which no WDH is visible (Fig. 8, grey 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 way beyond the influence of the coronagraph inner working angle), the raw contrast reached in the H-band is about 7.10 −5 (see also Vigan et al. 2015). In this on-sky image, the raw contrast is limited by the presence of speckles due to NCPA. As NCPA are always present and are the main limitation at 300 mas, this is the ultimate raw contrast we can reach under good observing conditions . From Fig. 8, the measured contrast curve (grey dash-dotted line), which was taken with a coherence time of 9 ms, is above that expected for the system with a 3 ms coherence time (green line). As such we expect the WDH to show up from a coherence time below 3 ms. Reported on the cumulative histogram in Fig. 3, it 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.
The WDH has a high contrast and is only unveiled thanks to the use of a coronagraph, as highlighted in Fig. 9, where no dif- ference can be seen in the two non-coronagraphic images (left), while the WDH is very bright in the coronagraphic image under short coherence time (bottom right). The corresponding raw contrast profiles are shown on Fig. 10, showing the impact of the WDH on real data. Hence this effect appears as an important limitation only for the latest generation of HCI instrument such as VLT/SPHERE, Gemini/GPI and Subaru/SCExAO using both extreme-AO and advanced coronagraph technology. The careful analyses of the WDH presented in this paper is thus motivated by the strong impact of the WDH on the achieved contrast with these high contrast instruments.

Analysis of the wind driven halo in the focal plane images of SPHERE
In this section, we present a method to analyze 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 analysis are conducted directly within the focal plane images as the results are more reliable in the context of high-contrast performance than using external data such as AO-telemetry or turbulence profiling. In addition to the analysis presented in this paper on how the WDH affects the contrast after post-processing, we intend to use this WDH analysis procedure for two studies: (i) making a statistical analysis of how the SPHERE data are affected by the WDH and correlate it with the AO telemetry and profiling data and (ii) estimate the WDH to remove it from the data. To illustrate and verify our approach, we use four types of multispectral coronagraphic images, following the integral field spectrograph (IFS, Antichi et al. 2009) of SPHERE in the YHband (from 0.96 to 1.66 µm, with a spectral resolution R ∼ 30): -Case 1: Simulated images obtained as a temporal stack of short exposures, containing only the AO-residuals due to fitting and servolag errors (produced using our analytical AO simulator) and using an apodized Lyot coronagraph (APLC, Martinez et al. 2009) as the one used on SPHERE (see Fig. 11, left);  Fig. 10. Raw contrast as a function of the separation to the star for the two images with (solid lines) and without (dotted-dash line) WDH presented in Fig. 9. The red line is 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 affected by the coronagraph (left) and outside the AO correction zone (right).

Raw contrast
-Case 2: Simulated images like before, additionally including NCPA upstream and downstream the coronagraph focal plane mask. In order to have realistic simulations, close to the images actually obtained with SPHERE, we used the upstream phase estimated by the ZELDA mask (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 direction for each pupil fragment separated by the spider of the telescope (see Fig. 11, middleright); Article number, page 7 of 17 A&A proofs: manuscript no. aanda -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 are stacked to obtain a long-exposure image. The injected wind direction is 125 degrees. For each multispectral cube, Fig. 11 shows only the shortest wavelength of the cube, at 0.966 µm.

Extracting the WDH contribution in 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 atmospheric turbulence temporal variation (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 DIT are long enough so that the WDH shows as a smooth structure. Therefore, 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) is to spatially filter the data in the Fourier domain to keep only its 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 apply an annular binary mask to cover the coronagraph signature (below 2λ/D) and the seeing-limited area (beyond 20λ/D). As SPHERE images show two artefacts on its correction ring, which are due to the DM square grid (see Fig. 1 right image, encircled with yellow dashed line), these two spots are 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 fraction (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 too smooth halo whereas from a filtering fraction of 15% on, it shows sharper structures due to the capture of speckles. For a filtering fraction of 5% the low-pass filtered image (bottom row) still shows a slight elongation along the wind speed direction at short separation, which completely disappears from 15% on. When the filtering fraction is too high, we can visualize very faint structures such as the grid of microlenses from the IFS design. 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.
To check whether this filtering procedure indeed brings up mainly, if not only, the WDH component, we applied it on the three simulated data cases and compared with the exact same simulated cases but without the servolag error included. 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 the ones 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 non-filtered 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 from a filtering fraction of 22% it reaches a plateau at the expected absolute values. When adding NCPA and LWE or LOR, 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 the case 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%. Note that this rule of thumb is only valid for SPHERE data and the optimum fraction must be determined for a given instrument.

Direction of the WDH
In order to assess the direction of the WDH elongation, we developed the following procedure. In a first step, we apply a discrete Radon transform 8 (Radon 1917(Radon , 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 average 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 perform a Gaussian fit around the maximum value of this Radon profile in order to extract the preferential direction of the WDH. Indeed, using directly the maximum value of the Radon profile is not the most robust option since, depending on the number of pixels within the masked image, the profile might be irregular. After testing different possibilities and data sets, we found out that performing a Gaussian fit around the maximum value of the Radon profile is the most robust method that does not require tuning any user-parameter. The uncertainties on the estimated direction are therefore the uncertainties on the Gaussian fit performed as, in an ideal case, the Radon profile is purely sinusoidal and the Gaussian fit around the maximum value is quite accurate. Thanks to 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 to extract the direction of the WDH on the first simulated data-set. Since in this data set only the WDH contribution is present, there is 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 systematic errors (such as centering of the raw image and plate scale) are dominating the error. In the following we will therefore ignore the uncertainty on the estimated direction.
As a second step, we apply 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), show the obtained Radon profiles (left) and the extracted WDH contour plots with the fitted direction overlaid (right). On the simulations including NCPA (case 2), the estimated direction is 125.1 degrees, as in the case without NCPA. The NCPA are mainly inducing high spatial frequencies in the focal plane (see Vigan et al. 2019), they do not affect greatly the extraction of the WDH 8 The discrete Radon transform is a linear operator that transforms a given 2D image into a 2D map showing the intensity along lines over the image, as a function of the angle (for a square image of dimension N img × N img , the Radon map is of dimension (π.N img ) × (2.N img − 1)).  and therefore of its direction. On the simulations including additionally LOR and LWE (case 3), the estimated direction is 124.5 degrees. Indeed these low spatial frequencies slightly offset the WDH, as shown in the Radon profile that 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 perfectly the inner part or the outer part, but is a good fit overall.
As a last step, we checked the robustness of the estimated direction with the filtering fraction 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, even without filtering, the estimation is off from the injected value by only two degrees. We repeated this experiment on on-sky data (for which the real wind direction is not known) showing different amount of WDH. Except for the case of very low WDH contribution, the estimated direction is stable to 1 degree for a filtering fraction between 5% and 25%. In the case of very low WDH, the estimation is dominated by residual tip-tilt errors and varies of up to 7 degrees between a filtering fraction of 5% to 17% but from 17%, a plateau is reached. Note that this experiment can be used to determined the optimal filtering fraction for a given HCI instrument. In the context of the present study and our two future applications, we would like to estimate the wind direction within a few degrees accuracy. The procedure to extract the direction of the WDH is therefore robust enough to further analyze the WDH structure.

Strength of the WDH
When no WDH is present in the image, the latter procedure provides a random direction. In order to check whether some WDH is present in the data, after the filtering and extraction of the Radon profile, we compute the standard deviation of the Radon profile. If this standard deviation is less or equal to the contrast expected without WDH, then there is no WDH. This threshold C τ is obtained empirically via the simulated PSD of the residual Fig. 13. Amount of light within the WDH structure extracted by spatial filtering, as a function of the filtering fraction (starting from 5%, at the vertical dotted black line). Red solid line is the first case containing only the WDH contribution, orange dashed line is the second case with NCPA and green dotted-dash 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 (resp. 73%, 66% and 64%). 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 is set to C τ = 7.10 −4 ). We checked this procedure on SPHERE on-sky data and it efficiently sorted the data without WDH from the data showing WDH.
More generally, in order to assess the amount of starlight scattered into the WDH within the corrected area, we define its strength, S WDH , which is 100% if all the starlight in the AOcorrected zone belongs to the WDH and 0% if there is no WDH: where σ(R WDH (θ)) is the standard deviation of the Radon profile R WDH extracted previously. Note that 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 applying this metric 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). From a filtering fraction of 10% on, this metric reaches a plateau at a strength of 95%. On simulated data containing only WDH, this metric provides a strength of 99.6%, and on simulated data without WDH, it provides a strength of 0.05%.
In order to check that the defined metric makes sense on real data, we sorted out six images within the on-sky data cube of 51 Eri, showing visually 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 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 obtained for the six cases is consistent with what is observed in the image.

Asymmetry of the WDH
To quantitatively characterize the asymmetry of the WDH, we define its asymmetry factor, A WDH , which is 0% when the WDH Estimation of the wind direction on 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 the case 1, a filtering fraction of 15% has been used to isolate the WDH. For the first three cases simulated, the green line shows the simulated wind direction (125 degree).
is perfectly symmetric and 100% when the WDH shows only one wing:  where I(r, θ) = I + (r, θ)+ I − (r, θ) is the total intensity contained in the WDH. The two wings, I + and I − , are obtained by cutting the image along the perpendicular direction 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%, whatever the filtering fraction. To verify that this metric makes sense, we simulated images exactly in the same conditions as for the third case (including NCPA, LOR and LWE) but with a varying amount of asymmetry (Fig. 17, bottom). We compare the extracted values from the images with the values obtained by analysing the PSD of the AO-residuals used to produce the simulated images (solid dashed lines in Fig. 17,   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, A WDH , defined at Eq. 5, as a function of the filtering fraction (starting at 5%, black dotted line). The blue solid line is for the image simulated without asymmetry. Bottom: Corresponding simulated images showing no (left) to large asymmetry of the WDH (right). The dashed lines correspond to the value extracted directly from the PSD. 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, from 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.
We checked how our estimation of the direction is affected by the asymmetry of the WDH. Due 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 less than 10 degrees).
This section confirms that the described procedure is valid to extract the WDH parameters directly 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 with stable results), in order to analyze the spatial, temporal and spectral behavior of the WDH.

Effect of the WDH on the final contrast after post-processing
Within the coronagraphic images sequence delivered by the latest generation of HCI instruments, the starlight residuals in the corrected area have a contrast level ranging from 10 −3 to 10 −5 .
To carve the starlight towards deeper contrast, advanced postprocessing techniques have been developed, relying 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 varying fast enough to be smoothly averaged and removed by an appropriate filtering. Those QSS are originating from NCPA and, for short exposure time, atmospheric residuals (e.g. Cantalloube et al. 2019).
Consequently the post-processing techniques that have been developed are aiming at removing these QSS but are not tailored for the WDH, the LWE or other source of errors whose temporal, spectral or spatial behavior differ from the QSS. As a consequence the contrast reached after post-processing is usually worse than expected from simulations. After 5-years of SPHERE operations, the median contrast reached at 500 mas is 2.10 −5 (Langlois et al., in prep.), instead of the 10 −6 expected from simulations before the commissioning of the instrument (Vigan et al. 2010).
The mainstream post-processing techniques used today rely on differential imaging that consist in: (1) estimating the QSS, (2) subtracting it from the image and (3) combining all the subtracted images available to increase the signal-to-noise-ratio of potential companions. In order to estimate the QSS, most HCI instruments working in 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. 2006) technique. If the instrument additionally provides simultaneous images at two or more wavelengths, one can use the spectral diversity 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. At last, if a sufficient number of images on various targets are available from the instrument, one can apply Multiple Reference Differential Imaging (MRDI, Lafrenière et al. 2009;Xuan et al. 2018;Ruane et al. 2019, with space-based instruments and with ground-based instruments resp.) using the correlation between the images to be processed and the library of images from the instrument.
For these three solutions, the estimated QSS is usually not perfectly estimated, which yields a high amount of speckle residuals, mostly at close separation to the star where the behavior of speckles is non-linear due to the coronagraphic device (for 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 a signal completely removed 9 . For point source detection, one can perform a high-pass spatial filtering to remove the light from the WDH, at a cost of a slight loss of signal that can be modeled and accounted for to characterize the candidate (Cantalloube et al. 2015). This is however not possible for extended sources, in which case most of the signal would be removed, in particular for low surface brightness disks seen face on. As a consequence, when comparing the number of disks detected in scattered light (about 40) to the number of potentially detectable disks with an infrared excess greater than 10 −4 from SPITZER data (Chen et al. 2014), we find a detection rate of ∼ 25% (Milli et al., in prep.). It is still unclear if 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 see how it affects the post-processed contrast and to understand better how it could be removed by a different approach.

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 lasts 75 min and the reduced cube consist of 64 frames that are each made of 4 binned images of 16 sec integration time (the total number of raw images is therefore 256). Four frames have been automatically rejected by the SPHERE reduction pipeline (Delorme et al. 2017) due to their bad quality.
The direction of the WDH is rotating 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 accurately its direction), we find that the slope is not perfectly equal to one but rather 1.3. This slight discrepancy is expected if the jet stream layer wind direction changes over the timescale of the observing sequence, as shown in Fig. 6 (bottom). The strength is varying erratically, from 60% to 88%, which corresponds to what we can observe by visualizing the data cube 10 . The asymmetry factor is decreasing slowly from 17% to almost no asymmetry, without 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 is always varying significantly from one frame to another. This prevents the ADI-based algorithm from capturing correctly the level of this feature in the model of the starlight residuals, resulting in the presence of 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. Thus, if one simply rotates and median combines the temporal data cube, the WDH signature co-aligns 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 applying a principal component analysis (PCA, Soummer et al. 2012;Amara & Quanz 2012), the intensity of the WDH residuals depends on the number of principal components kept to build the model. As shown in Fig. 19 (right), a more aggressive PCA removes low spatial frequencies but also considerably reduce the signal sought, mostly for extended features, and still leaves high residuals at close separation to the star. We additionally processed the 10 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 via these two external measurements. F. Cantalloube et al.: The wind-driven halo in high-contrast images I: analysis from the focal plane images of SPHERE same data set with locally optimized combination of images 11 (LOCI, Lafrenière et al. 2007) and Non-negative matrix factorization 12 (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 −3 at 300 mas, compared to 5.10 −6 when the WDH is filtered out.

Spectral variations of the WDH
To investigate the spectral variation of the WDH, we run 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 is shown in Fig. 11). Each cube is constituted 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 varying at most by 2 degrees for all of the four tested cases. As mentioned in Sect. 2, the turbulence coherence time τ 0 , through Eq. 1, scales with the wavelength as r 0 , that is to say as λ 6/5 . The variance of the AO residual phase due to the servolag error, through Eq. 2, scales as τ −5/3 0 . Therefore, the WDH strength scales with the wavelength as λ −2 , which is indeed what we observe in Fig. 20 (middle). The asymmetry factor is expected to steadily increase with the wavelength, as demonstrated in . This is indeed what is observed in Fig. 20 (right) where, in the four cases, the asymmetry of the WDH follows the wavelength variation.
The widely accepted model to perform SDI is, assuming that the aberrations are achromatic and that the small phase approximation is valid, that 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), via the square of the ratio between the two wavelengths: s λ 2 ( λ 2 λ 1 r) = ( λ 1 λ 2 ) 2 × s λ 1 (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). Since SDI assumes that the QSS varies as λ 2 whereas the WDH 11 As implemented in the SpeCal pipeline dedicated to post-process SPHERE infrared data, described in Galicher et al. (2018). 12 As implemented in the VIP open-source pipeline gathering various state-of-the-art post-processing methods, described in Gonzalez et al. (2017). varies as λ −2 and its asymmetry as λ, a SDI post-processing 13 leaves a strong asymmetric WDH imprint in the residual map, as shown in Fig. 21. This imprint gets worse as the reference channel is further from the science channel since the rescaling square law deviates more from the actual WDH spectral variation. After classical SDI, the typical contrast of the WDH residual (as in Fig. 21 left) is 10 −4 at 300 mas, compared to 6.10 −5 when the WDH is filtered out.

Spatial variations of the WDH
As discussed in Sect. 2, the shape, direction and intensity of the WDH signature depends on various external atmospheric turbulence parameters, on the observed target and on the observation set-up. In the context of RDI, 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 observe again strong WDH residuals in RDI post-processed images, as illustrated in Fig. 22. The images shown were processed using an algorithm implementing the PCA, as described in . 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 which did not include the science target. This library was build from targets part of the SHARDDS program survey for debris disk 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). We see that WDH residuals still remain in the RDI-PCA postprocessed images, mostly at large distance where it is known that RDI reaches a turnover (Ruane et al. 2019) and shows an equal or higher level of residuals compared to 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 how the reference to be subtracted is picked within the library of images. Usually the highest correlation is set on the QSS and not the WDH. Also using larger library should leave more chances to get a reference that also contains a similar WDH signature. Using a basic RDI implementation, the contrast loss due to the WDH residual signature is of about an order of magnitude in the affected regions.
In a next paper we will propose a method to reduce the impact of the WDH on the contrast after differential imaging, without affecting close-in and/or extended signals. The idea is to use Fig. 19. Post-processed image of the 51 Eri SPHERE data cube (as in Fig. 11, left) using ADI-based algorithms. Left: classical-ADI (cADI Marois et al. 2006) for which the temporal median is used as a QSS field model. Middle-left: Principal component analysis (PCA Amara & Quanz 2012;, for which a linear combination of the 3 first principal components is used as a QSS field model. Middle-right: PCA using 5 components. Right: PCA using 10 components. About one order of magnitude in contrast is lost due to the WDH residuals.  Post-processed images of the 51 Eri data using SDI. Left: the reference is the 10 th channel (1.12 µm), which is rescaled and removed from the frame in the first channel 1.00 µm. Right: the reference is this time the 20 th channel (1.31 µm). About half an order of magnitude in contrast is lost due to the WDH residuals. the analysis procedure described in Sect. 3 to estimate a model of the WDH and subtract it so as to recover a better contrast, intrinsically limited by the NCPA.

Conclusions and perspectives
In this paper, we conducted a study of the wind driven halo 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 are playing a role in building up this WDH, from the instrument design to its interaction with the observing conditions and op-erations. 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 refined by an upcoming detailed analysis of the SPHERE-SHINE survey data sample.
To that end, we proposed a procedure to analyze the WDH contribution directly from the focal plane images. We checked, thanks to various simulations, that this procedure is able to extract the relevant parameters to describe 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 correlate the results to the atmospheric turbulence parameters, essentially the turbulence coherence time and the scintillation measured by the MASS-DIMM and Stereo-SCIDAR at Paranal observatory, and by the AO telemetry data. Such a study will bring a deeper understanding of this specific limitation, towards proposing a way to alleviate this frequent contrast limitation on the data already 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+, 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 to analyze the WDH contribution, we highlighted its effect on the final contrast after postprocessing. 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 de- 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. tection of close planets and/or circumstellar disks. As such, 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, function of the three parameters discussed in the present paper, to fit the WDH signature directly 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 COF-FEE Herscovici-Schiller et al. 2019) are aiming at estimating the phase responsible for the observed coronagraphic image (using respectively the spectral diversity of an IFU or induced phase diversity). The AO residuals are taken into account in the model via the turbulence phase structure function. In that context, studying and understanding the WDH is crucial to obtain a correct model of the AO residuals for such algorithms to be operational on sky. These type of algorithms are relying on a full model of the HCI instrument and reach theoretical contrasts close to the photon noise limit, making them potentially the next big step to post-processing techniques.