# Coherent differential imaging of high-contrast extended sources with VLT/SPHERE ⋆

Axel Potier, Raphaël Galicher, Pierre Baudoz, Johan Mazoyer, Zahed Wahhaj, Ruben Tandon, Jonas G. Kühn, Laura Perez, Gael Chauvin

► To cite this version:

Axel Potier, Raphaël Galicher, Pierre Baudoz, Johan Mazoyer, Zahed Wahhaj, et al. Coherent differential imaging of high-contrast extended sources with VLT/SPHERE ⋆. Astronomy & Astrophysics - A&A, 2025, 704, pp.A98. 10.1051/0004-6361/202556606. hal-05350745

HAL Id: hal-05350745

https://hal.science/hal-05350745v1

Submitted on 12 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Coherent differential imaging of high-contrast extended sources with VLT/SPHERE*

Axel Potier<fnref n="1" />, Raphaël Galicher<fnref n="1" />, Pierre Baudoz<fnref n="2" />, Johan Mazoyer<fnref n="2" />, Zahed Wahhaj<fnref n="3" />, Ruben Tandon<fnref n="4" />, Jonas G. Kühn<fnref n="4" />, Laura Perez<fnref n="5" />, and Gael Chauvin<fnref n="6" />,<fnref n="7" />

Received 28 July 2025; accepted 19 October 2025

## ABSTRACT

**Context.** High-contrast imaging relies on advanced coronagraphs and adaptive optics (AO) to attenuate the starlight. However, residual aberrations, especially non-common path aberrations between the AO channel and the coronagraph channel, limit the instrument performance. While post-processing techniques such as spectral or angular differential imaging (ADI) can partially address those issues, they suffer from self-subtraction and inefficiencies at small angular separations or when observations are conducted far from transit.

**Aims.** We previously demonstrated the on-sky performance of coherent differential imaging (CDI), which offers a promising alternative. It allows for isolating coherent starlight residuals through speckle modulation, which can then be subtracted from the raw images during post-processing. This work aims to validate a CDI method on real science targets using VLT/SPHERE, demonstrating its effectiveness in imaging almost face-on circumstellar disks, which are typically challenging to retrieve with ADI.

**Methods.** We temporally modulated the speckle field in VLT/SPHERE images, applying small phase offsets on the AO deformable mirror while observing stars surrounded by circumstellar material: HR 4796A, CPD-36 6759, HD 169142, and HD 163296. We hence separated the astrophysical scene from the stellar speckle field, whose lights are mutually incoherent.

**Results.** Combining a dozen of data frames and reference coronagraph point spread functions through a Karhunen–Loève image projection framework, we recover the circumstellar disks without the artifacts that are usually introduced by common post-processing algorithms (e.g., self-subtraction).

**Conclusions.** The CDI method therefore represents a promising strategy for calibrating the effect of static and quasi-static aberrations in future direct imaging surveys. Indeed, it is efficient, does not require frequent telescope slewing, and does not introduce image artifacts to first order.

**Key words.** instrumentation: adaptive optics – instrumentation: high angular resolution – techniques: high angular resolution

## 1. Introduction

High-contrast imaging aims at the detection of light emitted, reflected, or scattered by any object (exoplanet, disk) in circumstellar environments. It potentially allows for a precise spectral, astrometric, and polarimetric characterization of any detected object. Such observations require cutting-edge coronagraph technologies and state-of-the-art wavefront sensing and control (WS&C) capabilities (Galicher & Mazoyer 2024) that are combined in dedicated facilities such as VLT/SPHERE (Beuzit et al. 2019), Gemini/GPI (Macintosh et al. 2014), Subaru/SCExAO (Lozi et al. 2018), and Magellan/MagAO-X (Males et al. 2020). However, the performance of these instruments remains limited by residual stellar light reaching the science detector, thereby drowning the exoplanet signal. These residuals are mainly caused by the spatiotemporal limitations of the extreme

adaptive optics (XAO) system, creating a smooth halo topped with stochastic stellar speckles in the long-exposure coronagraph image. The topping speckles, which result from unaveraged AO residuals, vary from one exposure to the other. The longer the exposure, the fainter the speckles (Singh et al. 2019). Additionally, the XAO wavefront sensor (WFS) is insensitive to the optical errors located in the science channel after the beam splitter and is biased by the aberrations in the WFS channel, leaving residual aberrations at the coronagraph. These non-common path aberrations (NCPA, Fusco et al. 2006) subsequently introduce additional stellar leakage in the form of static or quasi-static speckles. It is also worth noting that amplitude aberrations caused by beam transmission inhomogeneities or out-of-pupil phase errors are not typically corrected by conventional XAO WS&C algorithms. These amplitude aberrations also induce quasi-static stellar speckles in the coronagraphic image.

The active correction of these quasi-static speckles is a dynamic field of research. Techniques based on a focal plane WFS

<footnote n="1">¹ LIRA, Université Paris Cité, Observatoire de Paris, Université PSL, Sorbonne Université, CNRS, F-92190 MEUDON, France</footnote>

<footnote n="2">² LIRA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, 5 place Jules Janssen, 92195 Meudon, France</footnote>

<footnote n="3">³ European Southern Observatory, Alonso de Cordova 3107, Vitacura, Santiago, Chile</footnote>

<footnote n="4">⁴ Division of Space and Planetary Sciences, University of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland</footnote>

<footnote n="5">⁵ Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="6">⁶ MPIA, Max-Planck-Institut für Astronomie, Königstuhl 17, Heidelberg, D-69117, Germany</footnote>

<footnote n="7">⁷ Laboratoire J.-L. Lagrange, Université Cote d'Azur, CNRS, Observatoire de la Cote d'Azur, 06304 Nice, France
e-mail: axel.potier@obspm.fr</footnote>

<footnote>* Based on observations collected at the European Southern Observatory, Chile, 113.26G2</footnote>
are being tested or implemented on ground-based instruments. Some of these methods are implemented on a daily basis, leveraging the calibration source to estimate and correct for the NCPAs (Sauvage et al. 2007; N'Diaye et al. 2016; Lamb et al. 2018). However, once on-sky, these pre-calibrations are limited by differential aberrations between the internal source unit and the science path (Vigan et al. 2019; Potier et al. 2022b) as well as by a temporal decorrelation of the NCPAs (Milli et al. 2016; Vigan et al. 2022). Therefore, the on-sky calibration (directly using the photons from the science target star) of the NCPAs was tested with various instruments (e.g., Martinache et al. 2014; Galicher et al. 2019; Vigan et al. 2019; Bos et al. 2021; Skaf et al. 2022; Xin et al. 2023; Kueny et al. 2024). In particular, a dark hole demonstration – where destructive interferences are set up on a region of a detector – has been performed in narrowband (Potier et al. 2022a) and broadband light (Galicher et al. 2024). However, in poor observing conditions (meaning unstable AO correction), the dark hole technique becomes inefficient, as wavefront sensing must be repeated several times at each iteration, requiring a stable AO correction over several minutes.

Because of these drawbacks, the current strategy to mitigate speckles relies on post-processing approaches, each relying on such particular observational techniques as angular (ADI, Marois et al. 2006), spectral (SDI, Racine et al. 1999), polarimetric (PDI, Kuhn et al. 2001), and reference-star (RDI, Lafrenière et al. 2009) differential imaging. Building on such classical strategies, algorithms such as the Locally Optimized Combination of Images (LOCI, Lafrenière et al. (2007)) and the principal component analysis (PCA, Amara & Quanz (2012)), also known as Karhunen-Loève Image Processing (KLIP, Soummer et al. (2012)), have been developed. Additionally, numerous optimizations have been introduced, such as adjusting the size of the reduction zones (Marois et al. 2014; Ren et al. 2023) or modifying the selection of linear components (Gomez Gonzalez et al. 2016). Also, alternative approaches such as the patch covariance algorithm (PACO, Flasseur et al. (2018)) have been developed, which leverage prior knowledge of speckle noise distribution to better distinguish planetary signals. However, these techniques suffer from the self-subtraction of off-axis sources, especially at small angular separations (Esposito et al. 2014). Spectral-differential imaging self-subtraction depends on the separation and on the spectral properties of the object (Rameau et al. 2015) but can be mitigated through higher spectral bandwidth (Gerard et al. 2019a). Angular-differential imaging self-subtraction can be quantified and calibrated based on the evolution of the parallactic angle (PA) and the specific parameters used in the data reduction process. But, face-on disks remain undetected through ADI because of the centrosymmetric geometry of such systems. To address this limitation, forward modeling methods have been developed for point source objects (Marois et al. 2010; Pueyo 2016) or disk analysis (Mazoyer et al. 2020). However, in the case of disks, a wide range of possible morphologies must be explored, which is computationally intensive. Moreover, the diversity of complex and irregular structures cannot be fully captured by existing models.

More recently, coherent differential imaging (CDI) techniques – which should not suffer from these limitations – have been developed, leveraging the advent of focal plane WS&C algorithms cited earlier (Bottom et al. 2017; Gerard et al. 2018; Potier et al. 2022a). These methods rely on the coherent modulation (either spatial or temporal) of the speckle field. Only the coherent signal (starlight) responds to the modulation and, therefore, can be isolated from the astrophysical signal in post-processing. Contrary to Bottom et al. (2017), our previous

CDI demonstration using the VLT/SPHERE High Order Deformable Mirror (HODM) to sequentially modulate the speckle field (Potier et al. 2022a) was not validated on a science target. In the current paper, we validate the CDI algorithm, observing four circumstellar disks with VLT/SPHERE. In Sec. 2 we describe the acquisition of CDI sequences and their combination to isolate the disk signals in post-processing. In Sec. 3 we present the resulting images and derive intensity measurements to analyze the performance of our methods. Finally, we discuss the advantages and limitations of the algorithm as well as potential improvements in Sec. 4.

## 2. Method

### 2.1. Single CDI sequence workflow

The CDI method based on pair-wise probing (PWP) has been intensively described in the literature (see, e.g., Give’On et al. 2007a; Potier et al. 2020a,b), and its equations are detailed in Appendix A. It relies on a temporal modulation of the speckle field using the HODM. The diversities (or probes) are applied to the HODM by modifying the WFS reference slopes with the XAO system running in closed-loop. Contrary to the study performed in Potier et al. (2022a), three pairs of probes are used to cover the full 360° field of view. For the reasons described in Laginja et al. (2025) (i.e., reduction of mechanical constraints and nonlinear coupling, knowledge of the HODM influence functions, and mitigation of algorithm nonlinearities), we consecutively poked a triplet of neighboring single-actuators distributed along both the HODM axes, with a peak-to-valley amplitude of 400 nm. Such an amplitude was chosen to obtain a sufficient signal-to-noise ratio on the probe images while remaining in the linear regime defined in Appendix A. For each applied diversity, one image was acquired on the science detector IRDIS (Dohlen et al. 2008); hence, six images were required to estimate the modulated signal. Their combined exposure time was arbitrarily chosen to be equal to the exposure time of the raw signal with no probe applied. Therefore, a CDI sequence took twice as long as a typical data acquisition process. We propose some workarounds in Sec. 4.

We first estimated the speckle electric field (E-field) through the PWP algorithm, as described in Appendix A. The square modulus of the field, called $I_{PWP,i}$, was then subtracted from the total intensity image (obtained with no probe) $I_{tot,i}$ of the ith CDI sequence. Fig. 1 represents an illustration of one CDI sequence. The raw total intensity image consists of: 1) a smooth halo of XAO residuals with a coherence time below 10 ms (see Appendix B), overlaid with topping speckles; 2) individual static speckles induced by NCPAs, amplitude aberrations, and the diffraction pattern of the coronagraph; and 3) the astrophysical signal to be extracted. Because the CDI algorithm operates with seconds-long probe images to average out the intensity variations caused by atmospheric turbulence (Singh et al. 2019) (and achieve a sufficient S/N), the XAO halo becomes incoherent and cannot be sensed by PWP. Static speckles alone ($I_{PWP,i}$) are therefore estimated in the center image of Fig. 1. Subsequently, the resulting image after one CDI sequence ($I_{CDI}$) contains the smooth XAO halo and the astrophysical signal with no self-subtraction.

### 2.2. Combination of CDI sequences

The CDI algorithm can be performed in parallel to the dark-hole algorithm, but only on an individual sequence, since the
<figure><img src="image_2.png" /><figcaption>Fig. 1. Image decomposition of HD 163296 through CDI. Left: Coronagraph image with no probe (i.e., total intensity image). Center: PWP-estimated speckle field. Right: CDI result (estimated astrophysical scene subtracted from the total intensity image) after one sequence of PWP. The negative values, indicated by black-colored regions, show local over-subtraction in areas where PWP cannot properly estimate the speckle field.</figcaption></figure>

speckle intensity distribution varies at each iteration (Potier et al. 2022a). However, standalone CDI (with a static speckle pattern, which is not the case during a dark hole convergence loop) enables the combination of reference images by relying on previous methods that improve reference differential imaging algorithms (Lafrenière et al. 2007; Soummer et al. 2012).

The method can first be generalized using the following mathematical notations. For N CDI sequences, the reference image $I_{ref,i}$ – containing the static speckle signal – is reconstructed and then subtracted from $I_{tot,i}$ in order to retrieve the astrophysical signal ($I_{CDI,i}$). To ensure the speckles remain quasi-static during observations, the command of the instrument derotator is set such that the pupil does not rotate. Pupil tracking mode also guarantees that the calibrated model of the instrument, represented by the coronagraph operator $C$ in Appendix A, remains consistent along all the sequences. Then, for the sake of increasing the S/N, $I_{CDI,i}$ is derotated by an angle $\theta_i$ corresponding to the PA at sequence $i$. The result is then averaged (the median can also be used) over the full sequence,

$$I_{CDI} = \frac{1}{N} \sum_{i=1}^{N} R_{\theta_i} [I_{tot,i} - I_{ref,i}], \quad (1)$$

where $R$ corresponds to the rotation matrix, while $I_{CDI}$ is the result of the algorithm. The main objective of the method is to determine $I_{ref,i}$ as accurately as possible with respect to the static stellar residuals, using the sequence of $I_{PWP,i}$.

### 2.2.1. Batch process

Since PWP aims to retrieve the stellar speckle E-field alone, the most straightforward algorithm is a simple batch process, where the reference image at each sequence is simply the squared modulus of the PWP E-field estimation for the same sequence:

$$I_{ref,i} = I_{PWP,i}. \quad (2)$$

To illustrate its limitations, the final result of this option is shown in Fig. 2 (top right) and compared to more advanced approaches (see Sec. 2.2.2 and 2.2.3). Although the image quality is improved with respect to the noADI case ($I_{ref,i}$ =), speckle residuals are not perfectly canceled. We suspect that the inaccuracy of $I_{PWP,i}$ is due to PWP model uncertainties and nonlinearities,

<figure><img src="image_3.png" /><figcaption>Fig. 2. $I_{CDI}$ results for different reference estimations. Top left: $I_{ref,i}$ = 0, no CDI post-processing applied (equivalent to noADI). Top right: $I_{ref,i}$ = $I_{PWP,i}$, a simple batch process used. Bottom left: Karhunen-Loève coherent component projection. Bottom right: Karhunen-Loève image projection. All images have been high-passed filtered, as explained in Sec. 2.3.</figcaption></figure>

photon noise, as well as strong variations in atmospheric turbulence – not in a steady state – that bias the recording of probe images (see, e.g., the bright estimated speckles at small angular separations that over-subtract the total intensity image in Fig. 1). As for ADI or RDI, more advanced sequence combinations can therefore be key to mitigating the latter two noise components in the reference image. The number of required CDI sequences is discussed in Sec. 4.2.
### 2.2.2. Karhunen–Loève coherent component projection

We now assume a library of reference coherence intensity components obtained with PWP, whose singular value decomposition is given by

$$M_{PWP} = [I_{PWP,1}, I_{PWP,2}, \dots, I_{PWP,N}] = U \Sigma V^T, \quad (3)$$

where $U$ and $V$ are two square matrices with orthogonal vectors, called the left and right singular vectors, respectively, and $\Sigma$ is a rectangular diagonal matrix containing the singular values $\sigma_i$ of $M_{PWP}$, sorted in ascending order. In this framework, the principal components or eigenvectors of $M$, described in the pixel space, are the vectors of $V^T$. They are shown in Fig. 3 for the observation of HD 163296 (16 PWP sequences, see Tab. 1). As expected the first eigenvectors exhibit a speckle-shaped intensity distribution, while the last eigenvectors display a noisier distribution at the pixel scale. It means that truncating the highest $N - K_{klip}$ eigenvector components in the $I_{PWP,i}$ preserves the speckle intensity distribution displayed in most of the PWP estimates while reducing high-frequency noise. Mathematically, using the notations from Soummer et al. (2012), the reference images are the projection of the PWP estimate onto the truncated PWP eigenvector basis:

$$I_{ref,i}(K_{klip}) = \sum_{k=1}^{K_{klip}} < I_{PWP,i}, V_k^T > V_k^T = (U \Sigma K_{klip} V^T)|_i, \quad (4)$$

where

$$\Sigma_{K_{klip}, i,j} = \begin{cases} \Sigma_{ij} & \text{if } i < K_{klip} \\ 0 & \text{otherwise.} \end{cases} \quad (5)$$

The $K_{klip}$ threshold value is then optimized with respect to the current dataset and observing conditions. We note that such an analysis can also be directly performed on a library of probe image differences.

The former projection aims at tracking and calibrating speckle intensity modes that vary at minute timescales along the CDI sequences. However, if pixel-scale noise is significant, the solutions can simply be averaged over all PWP sequences to be used as the image reference:

$$I_{ref,i} = \frac{1}{N} \sum_{i=1}^{N} (U \Sigma_{K_{klip}} V^T)|_i. \quad (6)$$

Such a solution is more robust to noise, but it only targets perfectly static speckles throughout the sequences.

### 2.2.3. Karhunen–Loève image projection

Instead of projecting the eigenvectors computed in Eq. 3 onto the coherent signal, the standard approach is the Karhunen–Loève image projection (KLIP, Soummer et al. 2012), where the eigenvectors are directly projected onto each total intensity image:

$$I_{ref,i}(K_{klip}) = \sum_{k=1}^{K_{klip}} < I_{tot,i}, V_k^T > S V_k^T, \quad (7)$$

where $S$ corresponds to the search area. However, since this method aims to project static speckle components onto images that also contain the astrophysical source and a bright AO halo, it may result in an over-subtraction of the static speckle in order to match the halo brightness. Hence, $S$ may be chosen in areas

<figure><img src="image_4.png" /><figcaption>Fig. 3. Illustration of the $M_{PWP}$ eigenvector decomposition. Top row: Principal component 1, 2, and 3. Bottom row: Principal component 5, 10, and 15.</figcaption></figure>

where the intensity of the static speckles dominates, i.e., where the turbulence halo is the faintest typically in the outer part of the HODM influence function. In this paper, we chose $S$ to be an annulus, ranging from 551 to 919 mas for all the disks. On the other hand, KLIP is convenient because it enables the use of such widely employed advanced tools and software as VIP (Gomez Gonzalez et al. 2017; Christiaens et al. 2023).

### 2.3. High-pass filter

We explained why the technique can neither target the XAO smooth halo nor its topping speckles. As mentioned the topping speckle intensity is minimized while recording long exposures. The XAO halo remains in the post-processed images. To extract the astrophysical signal, this halo could potentially be estimated through AO telemetry and then filtered in the future. Instead, assuming the halo to be smoother than the astrophysical signal, the processed image $I_{CDI}$ can be high-pass filtered, for instance using a Gaussian kernel,

$$I_f = \mathcal{G}[I_{CDI}(K_{klip}), \sigma_{hp}], \quad (8)$$

where $\mathcal{G}$ represents the high-pass Gaussian filter operator and $\sigma_{hp}$ is the standard deviation of the Gaussian kernel. Contrary to the previous steps of the reduction, this last procedure causes the subtraction of the astrophysical signal, especially for large structures, resulting in a reduced throughput. $\sigma_{hp}$ therefore must be chosen wisely, depending on the structures one aims to image. Sec. 4.1 discusses its repercussions more quantitatively for our dataset. It is worth noting that alternative CDI methods using faster science detectors may be able to calibrate fast speckles induced by atmospheric turbulence (Kühn et al. 2018; Gerard et al. 2019b; Thompson et al. 2022).

## 3. Results of disk imaging

Coherent differential imaging was performed on VLT/SPHERE on June 4, 2024. The observing conditions are described in Appendix B. Four targets with known circumstellar and almost face-on disks were observed throughout the night. Their respective data acquisition is detailed in Tab. 1. In addition to the set of N-acquired CDI sequences, $N_{filt}$ CDI sequences were down-selected by eye for post-processing, getting rid of the few sequences where bursts of seeing occurred and led to either a low-
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of disk images obtained with ADI (left column) and CDI (center column). For each disk, the ADI and CDI images are displayed using the same color scale; however, the scales vary from one disk to another. ADI and CDI results were obtained through the same dataset. We use K-klip ADI (Soummer et al. 2012) with $K_{klip} = N_{filt}/2$. Right column: Comparison of disk intensities, for pixels inside the region encircled in green, obtained with both ADI and CDI as a function of the position angle.</figcaption></figure>

quality total intensity image or a misestimate of the coronagraph
point spread function (PSF) through PWP.

The final post-processed images obtained with Eqs. 7 and 8
are shown in Fig. 4 (central column) together with the ADI re-
sults obtained through the same dataset, discarding probe images
<table><caption>**Table 1.** Observations from 2024-06-04</caption><thead><tr><th>System</th><th>Obs. date<br>(UTC)</th><th>t<sub>exp</sub><br>(s)</th><th>t<sub>p</sub><br>(s)</th><th>N</th><th>N<sub>filt</sub></th><th>ΔPA<br>(degree)</th></tr></thead><tbody><tr><td>HR 4796A</td><td>23:46:13</td><td>96</td><td>16</td><td>16</td><td>9</td><td>57.7</td></tr><tr><td>CPD-36 6759</td><td>01:12:46</td><td>96</td><td>16</td><td>20</td><td>15</td><td>52.5</td></tr><tr><td>HD 169142</td><td>03:36:10</td><td>192</td><td>32</td><td>16</td><td>11</td><td>30.4</td></tr><tr><td>HD 163296</td><td>06:04:23</td><td>96</td><td>16</td><td>16</td><td>14</td><td>14.7</td></tr></tbody></table>

**Notes.** The columns give the host star name, the time the observation began, the exposure time $t_{exp}$ of the no-probe coronagraphic images, the exposure time $t_p$ of the probe images, the number $N$ of CDI sequences, the $N_{filt}$ of selected CDI sequences for data processing, and $\Delta PA$, the total change of PA throughout the observing sequences.

for comparison (left column). On one hand, in addition to the speckle intensity being minimized with respect to noADI (see in Fig. 2), CDI results show that the disk features are well recovered by the algorithm, while self-subtraction is absent to first order – even for face-on disks. On the other hand, because ADI tends to subtract the signal of interest (Milli et al. 2012; Ruane et al. 2019; Xie et al. 2022), almost no signal is detected for three of our targets using ADI, and the HR 4796 disk image is impacted by self-subtraction particularly at small angular separations.

To quantitatively compare the performance of CDI with respect to ADI in retrieving disk structures, we measured the normalized intensity of each pixel within the green area overplotted in Fig. 4, which corresponds to the estimated disk locations. These intensities were then plotted as a function of their respective position angles in Fig. 4 (right column). The results show that in all cases, ADI suffers from self-subtraction — significantly reducing the disk signal — although CDI remains largely unaffected. For example, the mean intensity ratio between CDI and ADI across these pixels reaches 216 for HD 163296. In the case of HR 4796A, although the average ratio is lower (approximately 7) due to the disk’s inclination, which makes it less vulnerable to ADI self-subtraction, and poorer observing conditions, values as high as ~100 are observed for pixels closest to the disk center.

Regardless, all these results demonstrate that CDI is a promising alternative to ADI, mostly free of any self-subtraction artifacts, especially for spatially extended sources such as disks. Coherent differential imaging is also a serious contender to RDI, as it requires neither extra telescope slewing and exposure time on a nonscience target (Wahhaj et al. 2021), nor the need to build up an archival PSF library (Ruane et al. 2019; Xie et al. 2022; Ren et al. 2023). However, without further observations of targets devoid of extended structures, deriving reliable and general contrast curve comparisons with other post-processing methods is difficult. As a result, it remains uncertain whether CDI could serve as a viable replacement for ADI in point-source detection. This question will be addressed in future follow-up studies. Also, a deep comparison with alternative post-processing techniques such as RDI (e.g., Xie et al. 2022), PDI (e.g., Gratton et al. 2019; Ren et al. 2023), or advanced ADI pipelines designed to reduce self-subtraction of face-on disks (e.g., Juillard et al. (2024); Flasseur et al. (2024)), as well as the detailed analysis of disk structures or single sources in these images is beyond the scope of this paper.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Effect of the high-pass filter on the off-axis source throughput for various Gaussian kernel standard deviations.</figcaption></figure>

# 4. Discussion

## 4.1. Influence of the high-pass filter

The main performance limitation of the current algorithm is its inability to tackle speckles whose lifetimes are smaller than the total exposure time of one CDI sequence. In particular, the dynamic speckles induced by the atmospheric turbulence are averaged into a smooth halo in long exposure images and filtered with a high-pass filter (see in Sec. 2.3). However, this final process also filters out the astrophysical signal that was unaffected until this stage. Fig. 5 shows the influence of the high-pass filter on the entire algorithm transmission – calculated with VIP – for the HR 4796A disk-free dataset. Indeed, to prevent the disk signal from interfering with the contrast measurement, the disk was forward-modeled and subtracted from the dataset, according to the disk parameters detailed in Milli et al. (2017) and as described in Appendix C. While the planet transmission is scarcely impacted at any separation without high-pass filtering – indicating the absence of self-subtraction with CDI – the potential planet throughput decreases with smaller $\sigma_{hp}$. For instance, it decreases by ~60% for $\sigma_{hp} = 3$ pixels, which is close to the PSF size. Other approaches to estimate and subtract this halo – such as using WFS telemetry – would therefore be essential to preserve the main benefit of CDI. SPHERE/IRDIS has a minimal exposure time of ~1s, but other CDI strategies using faster science detectors have been proposed to enable quicker estimation of the speckle field, targeting dynamic speckles caused by atmospheric turbulence (Kühn et al. 2018; Gerard et al. 2019b).

## 4.2. Result versus number of CDI sequences

It can be argued that the acquisition of the PWP library is a time-intensive process to the detriment of astrophysical signal acquisition. We show in Fig. 6 the final combined-CDI images obtained throughout the $N$ CDI sequences, which are compared with the noADI case ($I_{ref,i} = 0$ in eq. 1, Galicher et al. (2018)). It demonstrates that CDI leads to better performance for any number of sequences, particularly during the first few, when the static speckles do not average with diurnal motion. This means the CDI strategy can be quickly applied on targets of interest – much faster than ADI, which requires the planet image to rotate in the field of view (FoV), especially at low angular separations – such as during surveys when many systems must be explored within a limited time.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** NoADI (top row), CDI through KL image projection (center row), and single-PWP scheme CDI (bottom row) final images with respect to the number of considered frames: one, three, six, and nine sequences.</figcaption></figure>

Additionally, only a few CDI schemes were acquired at the beginning of the observation sequence, followed by more standard acquisitions, assuming the speckle field remains static over time. This option, called single-PWP scheme CDI, is also illustrated in Fig. 6, where it is assumed that only one PWP intensity estimate is performed at the beginning of the sequence to retrieve the astrophysical scene (i.e., $I_{ref,i} = I_{PWP,i}$ for all i). First, we see that this high-speed speckle calibration constantly outperforms the noADI case. Also, the final image after one sequence is identical to the KL-CDI result. Although the speckle subtraction performance degrades afterward, it saves 6N probe acquisitions, which in our case can be utilized to double the exposure time dedicated to image acquisition.

### 4.3. CDI self-calibration

In addition to the former point, CDI sequences can self-calibrate, with the astrophysical signal extracted from the probe images themselves, similarly to phase closure in interferometry. In brief, after the speckle E-field $E_S$ is estimated from eq. A.9, it can be added to or subtracted from the modeled probe E-field $iC[A\psi_m]$ to calibrate both the positive and negative probe images without the astrophysical signal $|E_S \pm iC[A\psi_m]|^2$. Subsequently from Eq. A.1, the astrophysical image $I_a$ can be retrieved via

$$I_a = I_{m\pm} - |E_S \pm iC[A\psi_m]|^2, \quad (9)$$

where $I_{m\pm}$ are the probe images. The full mathematical derivation proposed in this paper to combine many CDI sequences is described in Appendix D. Such a strategy is demonstrated in Fig. 7, where the probe images were self-calibrated to extract the HD 163296 disk. We, however, observe an additional amount of speckle noise in the final image caused by misedtimated speckles that rotate in the FoV. This can be explained either by a model error occurring in $iC[A\psi_m]$ and/or by additional shot-noise error contained in the probe images. Nevertheless, CDI self-calibration can be a workaround to compensate for its slowness. It could also be used in parallel to dark hole algorithms to maximize science return in the context of space-based coronagraphs such as RST/CGI.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Comparison of signal extraction methods: the signal obtained from probed images (right, i.e. self-calibration using eq. D.6) versus conventional CDI (center, using eq. 7) and noADI (left). All results are high-pass filtered, as explained in Sec. 2.3.</figcaption></figure>

### 4.4. Current PWP downsides

#### 4.4.1. Nonlinearities

Reducing the duty cycle is not the only limitation of the proposed strategy. Coherent differential imaging is also altered by several other drawbacks. First, PWP is inherently limited by nonlinearities arising when the probe amplitude is too high (Groff et al. 2015). Here, we chose high-amplitude probes, such as maximizing the S/N in the probe images, while the dark hole (DH) correction performed in Potier et al. (2022a); Galicher et al. (2024) still converged after a few iterations. However, in the case of standalone CDI, the post-processing might directly be biased by the nonlinear factor, limiting its contrast gain. Subsequent work will aim to quantitatively assess the best choice of probe amplitude with respect to S/N while preserving a high duty cycle. Nevertheless, previous studies have shown that single-actuator probes are less affected by nonlinearities than other shapes (Laginja et al. 2025).

#### 4.4.2. Chromatism

Another downside is chromatism. Pair-wise probing usually works with narrowband filters, used subsequently, limiting the telescope duty cycle. For example, RST/CGI takes advantage of three ~ 18 nm filters that are employed sequentially to perform wavefront sensing before applying a larger 56.5 nm (10% bandwidth) imaging filter Cady et al. (2025). In the SPHERE context, the dual-band imager IRDIS enables the use of the H2 filter in parallel to H3 to increase the bandwidth to 7.8% at no cost. But, the current version of the PWP algorithm cannot directly work in high-bandwidth filters such as H-band filters (17.8%) with IRDIS. Future work will thus rely on the SPHERE integral field spectrograph (IFS) to expand the working bandwidth up to 41% (Mesa et al. 2015).

#### 4.4.3. Field of view

The final downside is the limited obtained FoV. While adding a third probe in the algorithm enabled reaching a 360° post-processed image, the outer working angle of the corrected region is limited to the HODM influence function. This is acceptable since it corresponds to the AO correction region, where planets will be searched for with future XAO systems (Boccaletti et al. 2020; Kasper et al. 2021).

## 5. Conclusion

In this paper, we validated the CDI method relying on the temporal modulation of the VLT/SPHERE XAO HODM on astrophysical targets. Proposing a framework similar to those used
in more classical post-processing algorithms (ADI and RDI), we combined a dozen sequences of reference images, each representing the speckle intensity distribution estimated by PWP and subtracted from raw images. By retrieving four almost face-on disks, we demonstrated that the post-processing technique is robust and does not suffer from self-subtraction, unlike commonly used ADI, while improving the image quality compared to noADI. Additionally, we discussed potential workarounds to increase the instrument duty cycle in future surveys with CDI and plans for future studies regarding the nonlinearity and limited bandwidth of the algorithm. We believe CDI is an efficient approach for direct imaging and advocate its use in the context of future facilities equipped with cutting-edge XAO instruments such as SAXO+ (Boccaletti et al. 2020; Mazoyer et al. 2024). It is worth noting that CDI should be even more efficient in stable conditions, such as those encountered in space-based instruments. Coherent differential imaging could then increase the duty cycle of such instruments, either by being applied while digging a DH directly on the target star (Laginja et al. 2025), and/or while stabilizing the DH (Redmond et al. 2024), for the coronagraph instrument on the Roman Space Telescope, or the Habitable Worlds Observatory.

**Acknowledgements.** The authors thank the anonymous referee for their valuable comments and suggestions, which have significantly improved this manuscript. This study was supported by the IdEx Université Paris Cité, ANR-18-IDEX-0001, and the Physics department of Université Paris Cité, and the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES as well as an ECOS-CONICYT grant (#C20U02). L.P., J.M., G.C. acknowledge support from Programa de Cooperación Científica ECOS-ANID ECOS200049. L.P. acknowledges support by the ANID BASAL project FB210003 and ANID FONDECYT Regular 1221442. J. G. Kühn and R. Tandon have received funding from the Swiss State Secretariat for Education, Research, and Innovation as a SERI-Funded ERC 2021 Consolidator Grant, project RACE-GO # M822.00084, following the discontinued participation of Switzerland to Horizon Europe.

## References

* Amara, A. & Quanz, S. P. 2012, MNRAS, 427, 948

* Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

* Boccaletti, A., Chauvin, G., Mouillet, D., et al. 2020, arXiv e-prints, arXiv:2003.05714

* Bordé, P. J. & Traub, W. A. 2006, ApJ, 638, 488

* Bos, S. P., Bottom, M., Ragland, S., et al. 2021, in SPIE, Vol. 11823, 118231E

* Bottom, M., Wallace, J. K., Bartos, R. D., Shelton, J. C., & Serabyn, E. 2017, MNRAS, 464, 2937

* Cady, E., Bowman, N., Greenbaum, A. Z., et al. 2025, JATIS, 11, 021408

* Christiaens, V., Gonzalez, C., Farkas, R., et al. 2023, JOSS, 8, 4774

* Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in SPIE, Vol. 7014, 70143L

* Esposito, T. M., Fitzgerald, M. P., Graham, J. R., & Kalas, P. 2014, ApJ, 780, 25

* Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, A&A, 618, A138

* Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2024, MNRAS, 535, 689

* Fusco, T., Rousset, G., Sauvage, J. F., et al. 2006, Optics Express, 14, 7515

* Galicher, R., Baudoz, P., Delorme, J. R., et al. 2019, A&A, 631, A143

* Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92

* Galicher, R. & Mazoyer, J. 2024, CRPhys, 24, 133

* Galicher, R., Potier, A., Mazoyer, J., et al. 2024, A&A, 686, A54

* Gerard, B. L., Marois, C., Currie, T., et al. 2019a, AJ, 158, 36

* Gerard, B. L., Marois, C., & Galicher, R. 2018, AJ, 156, 106

* Gerard, B. L., Marois, C., Galicher, R., et al. 2019b, arXiv e-prints, arXiv:1910.04554

* Give’On, A., Belikov, R., Shaklan, S., & Kasdin, J. 2007a, OE, 15, 12338

* Give’On, A., Kern, B., Shaklan, S., Moody, D. C., & Pueyo, L. 2007b, Proc. SPIE, 6691, 66910A

* Gomez Gonzalez, C. A., Absil, O., Absil, P. A., et al. 2016, A&A, 589, A54

* Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, AJ, 154, 7

* Gratton, R., Ligi, R., Sissa, E., et al. 2019, A&A, 623, A140

* Groff, T. D., Riggs, A. J. E., Kern, B., & Kasdin, N. J. 2015, JATIS, 2, 1

* Juillard, S., Christiaens, V., Absil, O., Stasevic, S., & Milli, J. 2024, A&A, 688, A185

* Kasper, M., Cerpa Urra, N., Pathak, P., et al. 2021, The Messenger, 182, 38

* Kueny, J. K., Van Gorkom, K., Kautz, M., et al. 2024, in SPIE, Vol. 13097, 1309761

* Kühn, J., Patapis, P., Lu, X., & Arikan, M. 2018, in SPIE, Vol. 10706, 107062N

* Kuhn, J. R., Potter, D., & Parise, B. 2001, ApJ Lett., 553, L189

* Lafrenière, D., Marois, C., Doyon, R., & Barman, T. 2009, ApJ, 694, L148

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

* Laginja, I., Baudoz, P., Mazoyer, J., et al. 2025, A&A, 698, A130

* Lamb, M., Norton, A., Macintosh, B., et al. 2018, in SPIE, Vol. 10703, 107035M

* Lozi, J., Guyon, O., Jovanovic, N., et al. 2018, in SPIE, Vol. 10703, SPIE, 1070359

* Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, PNAS, 111, 12661

* Males, J. R., Close, L. M., Guyon, O., et al. 2020, in SPIE, Vol. 11448, 114484L

* Marois, C., Correia, C., Galicher, R., et al. 2014, in SPIE, Vol. 9148, 91480U

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

* Marois, C., Macintosh, B., & Véran, J.-P. 2010, in SPIE, Vol. 7736, 77361J

* Martinache, F., Guyon, O., Jovanovic, N., et al. 2014, PASP, 126, 565

* Mazoyer, J., Arriaga, P., Hom, J., et al. 2020, in SPIE, Vol. 11447, 1144759

* Mazoyer, J., Goulas, C., Vidal, F., et al. 2024, in SPIE, Vol. 13096, 130969D

* Mesa, D., Gratton, R., Zurlo, A., et al. 2015, A&A, 576, A121

* Milli, J., Banas, T., Mouillet, D., et al. 2016, in SPIE, Vol. 9909, 99094Z

* Milli, J., Mouillet, D., Lagrange, A. M., et al. 2012, A&A, 545, A111

* Milli, J., Vigan, A., Mouillet, D., et al. 2017, A&A, 599, A108

* N'Diaye, M., Vigan, A., Dohlen, K., et al. 2016, A&A, 592, A79

* Potier, A., Baudoz, P., Galicher, R., Singh, G., & Boccaletti, A. 2020a, A&A, 635, A192

* Potier, A., Galicher, R., Baudoz, P., et al. 2020b, A&A, 638, A117

* Potier, A., Mazoyer, J., Wahhaj, Z., et al. 2022a, A&A, 665, A136

* Potier, A., Wahhaj, Z., Galicher, R., et al. 2022b, in SPIE, Vol. 12185, 1218568

* Pueyo, L. 2016, ApJ, 824, 117

* Racine, R., Walker, G. A. H., Nadeau, D., Doyon, R., & Marois, C. 1999, PASP, 111, 587

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

* Redmond, S. F., Pueyo, L. A., Por, E. H., et al. 2024, JATIS, 10, 049003

* Ren, B. B., Benisty, M., Ginski, C., et al. 2023, A&A, 680, A114

* Ruane, G., Ngo, H., Mawet, D., et al. 2019, AJ, 157, 118

* Sauvage, J.-F., Fusco, T., Rousset, G., & Petit, C. 2007, JOSAA, 24, 2334

* Singh, G., Galicher, R., Baudoz, P., et al. 2019, A&A, 631, A106

* Skaf, N., Guyon, O., Gendron, É., et al. 2022, A&A, 659, A170

* Soummer, R. 2005, ApJ, Letters, 618, L161

* Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, Letters, 755, L28

* Thompson, W., Marois, C., Singh, G., et al. 2022, in SPIE, Vol. 12185, 121852C

* Vigan, A., Dohlen, K., N'Diaye, M., et al. 2022, A&A, 660, A140

* Vigan, A., N'Diaye, M., Dohlen, K., et al. 2019, A&A, 629, A11

* Wahhaj, Z., Milli, J., Romero, C., et al. 2021, A&A, 648, A26

* Xie, C., Choquet, E., Vigan, A., et al. 2022, A&A, 666, A32

* Xin, Y., Xuan, J. W., Mawet, D., et al. 2023, JATIS, 9, 035001
## Appendix A: PWP algorithm

Building on the idea developed by Bordé & Traub (2006) and first described by Give’On et al. (2007b), the Pair-Wise probing (PWP) algorithm assumes low-amplitude aberrations and diversities relative to the wavelength. Let $\psi_m$ be the diversity introduced by a deformable mirror located in the pupil plane. The total intensity on the detector $I_m$ can then be written as

$$I_m = (E_S + iC[A\psi_m]) (E_S + iC[A\psi_m])^* + I_a, \quad (A.1)$$

where $E_S$ is the speckle E-field to be estimated, $C$ the linear operator that transforms a complex E-field in the pupil plane to a complex E-field in the focal plane through a coronagraph, and $A$ is the complex amplitude in the entrance pupil plane. $I_a$ represents the intensity of the circumstellar scene and the AO incoherent residuals. It fluctuates weakly with the applied probe and is incoherent with the stellar signal. If we record the pair of images $I_{m+}$ and $I_{m-}$, with respective diversities $+\psi_m$ and $-\psi_m$, the difference then becomes

$$I_{m+} - I_{m-} = 4 [\Re(E_S)\Re(iC[A\psi_m]) + \Im(E_S)\Im(iC[A\psi_m])], \quad (A.2)$$

where $\Re$ and $\Im$ denote the real and imaginary parts of a complex number, respectively. We note that this image difference explicitly removes the astrophysical component $I_a$, ensuring no circumstellar scene in the final reconstructed coherent image. This equation can be rewritten in matrix form:

$$I_{m+} - I_{m-} = 4 \begin{bmatrix} \Re(iC[A\psi_m]) & \Im(iC[A\psi_m]) \end{bmatrix} \begin{bmatrix} \Re(E_S) \\ \Im(E_S) \end{bmatrix}. \quad (A.3)$$

This system is undetermined and admits infinitely many solutions due to degeneracy between the real and imaginary parts of $E_S$. Generalizing to $j$ pairs of images, for each pixel $(u, v)$ on the detector,

$$D_{(u,v)} = 4M_{(u,v)}F_{(u,v)}, \quad (A.4)$$

where

$$D = \begin{bmatrix} I_1^+ - I_1^- \\ \vdots \\ I_j^+ - I_j^- \end{bmatrix}, \quad (A.5)$$

$$M = \begin{bmatrix} \Re(iC[A\psi_1]) & \Im(iC[A\psi_1]) \\ \vdots & \vdots \\ \Re(iC[A\psi_j]) & \Im(iC[A\psi_j]) \end{bmatrix}, \quad (A.6)$$

$$F = \begin{bmatrix} \Re(E_S) \\ \Im(E_S) \end{bmatrix}. \quad (A.7)$$

To reconstruct the E-field, we minimize the cost function:

$$J_{PW} = \min_F \|D - 4M \cdot F\|^2. \quad (A.8)$$

This inverse problem is linear. It can be solved using the Singular Value Decomposition (SVD), provided that the matrix $M$ is invertible, i.e., its determinant is nonzero:

$$\tilde{F}_{(u,v)} = \frac{1}{4} M_{(u,v)}^{\dagger} D_{(u,v)}. \quad (A.9)$$

Here, $\tilde{F}$ is the approximate reconstruction of $F$, and $M^\dagger$ is the Moore-Penrose pseudoinverse of $M$. The matrix $M$ is invertible at each pixel $(u, v)$, if there exists a pair of probes $(m, n)$ such that

$$\begin{aligned} \Re(iC[A\psi_m])_{(k,l)} \cdot \Im(iC[A\psi_n])_{(k,l)}, \\ - \Re(iC[A\psi_n])_{(k,l)} \cdot \Im(iC[A\psi_m])_{(k,l)} \neq 0. \end{aligned} \quad (A.10)$$

This condition means that at least two of the probes induce different E-fields at that location.

<figure><img src="image_9.png" /><figcaption>Fig. B.1. Observing conditions provided by the ESO UT3 visitor telescope report on June 4-5, 2024.</figcaption></figure>

## Appendix B: Observing sequences

The observations were conducted during the night of June 4, 2024. According to the Paranal DIMM-Seeing monitor (see in Fig. B.1), the atmospheric seeing ranged from 0.245” to 0.834” at 500 nm, improving rapidly to below 0.5” after 00:00 UTC. As a result, the conditions were less favorable during the acquisition of HR 4796A data compared to those for HD 163296. These values reflect exceptional observing conditions at Paranal. The atmospheric coherence time increased steadily over the course of the night, ranging from approximately 5 ms to nearly 10 ms during the full observing sequence. Wind conditions were stable, originating from the north with speeds remaining below 5 m.s⁻¹. Observations were carried out using the H3 filter ($\lambda_0$ = 1667 nm, $\Delta\lambda$ = 54 nm) in combination with an Apodized Pupil Lyot Coronagraph (APLC, Soummer 2005) in the APO1/ALC2 configuration, featuring a focal plane occulting mask with a diameter of 185 mas.

## Appendix C: Forward modeling of HR 4796A disk

Coherent-differential imaging performance gain is hard to assess in images with extended structures as obtained during our experiment. We therefore decided to produce a synthetic image cube of the disk around HR 4796 at specific position angles to simulate a real observation. This cube is then used to numerically subtract the disk signal from $I_{tot}$. We used the ray-tracing approach from the ScatteredLightDisk class in vip_hci. Disks parameters from Milli et al. (2017) have been used to feed the function and have been slightly adjusted to minimize the disk signal by eye in the resulting image. Final images with and without CDI obtained through such process are shown in Fig. C.1 before and after disk subtraction. The HR 4796A disk-free images are then used to quantify the throughput performance of the method.
<figure><img src="image_10.png" /><figcaption>Fig. C.1. Forward modeling of HR 4796A. Shown are NoADI (left column) and CDI (right column) results before (bottom row) and after (top row) synthetic HR 4796A disk subtraction. We note a slight over-subtraction on the forward side.</figcaption></figure>

## Appendix D: Framework for probe self-calibration

Coherent-differential imaging self-calibration sequences can also be combined to increase the S/N, in a similar framework as in Sec. 2.2.2. Here we assume three libraries containing the N difference images for each probes m, whose singular value decomposition is given by

$$ \mathcal{M}_m = [\Delta I_{m,1}, \Delta I_{m,2}, \dots, \Delta I_{m,N}] = \mathcal{U}_m Z_m \mathcal{V}_m^T, \qquad (D.1) $$

where $\Delta I_{m,j}$ represents the m$^{th}$ probe image difference for the j$^{th}$ CDI sequence. $\mathcal{U}_m$ and $\mathcal{V}_m$ are the left and right singular vectors, while $Z_m$ are the singular values in ascending orders. This decomposition can be used to filter out the noise components in the $\Delta I_m$, assuming noise to be represented in the $(N - K_{klip})$ last principle components of $\mathcal{M}_m$. Mathematically, $\Delta \mathcal{I}_{m,j}(K_{klip})$, the m$^{th}$ filtered probe image difference for the j$^{th}$ CDI sequence then becomes

$$ \Delta \mathcal{I}_{m,j}(K_{klip}) = \sum_{k=1}^{K_{klip}} \langle \Delta I_{m,j}, \mathcal{V}_{m,k}^T \rangle \mathcal{V}_{m,k}^T = (\mathcal{U}_m Z_{m,K_{klip}} \mathcal{V}_m^T)|_j. \qquad (D.2) $$

As in Appendix A, the matrix of filtered probe image difference $\mathcal{D}_j(K_{klip})$ is computed via

$$ \mathcal{D}_j(K_{klip}) = \begin{bmatrix} \Delta \mathcal{I}_{1,j}(K_{klip}) \\ \vdots \\ \Delta \mathcal{I}_{m,j}(K_{klip}) \end{bmatrix}, \qquad (D.3) $$

to perform a new E-field estimation $\mathcal{E}_j(K_{klip})$ at each sequence j, leveraging the same M probe model matrix used in usual PWP:

$$ \mathcal{E}_j(K_{klip})_{(u,v)} = \frac{1}{4} M_{(u,v)}^\dagger \mathcal{D}_j(K_{klip})_{(u,v)}. \qquad (D.4) $$

Such an E-field can now be added to the probe E-field model to calibrate for any probe images taken during the N CDI sequences and retrieve the astrophysical signal $I_{a,mjl}$ from the probe image $I_{mjl}$:

$$ I_{a,mjl} = I_{mjl} - |\mathcal{E}_j(K_{klip}) + (-1)^l iC[A\psi_m]|^2, \qquad (D.5) $$

where $I_{mjl}$ is the m$^{th}$ probe of the j$^{th}$ CDI sequence, whose sign is defined with the argument $l$ ($l=1$ (resp. $l=2$) meaning negative (resp. positive) probe). We then chose to derotate each obtained solution while all N sequences of six probe images are averaged out to form the astrophysical scene $I_a$ via

$$ I_a = \frac{1}{6N} \sum_{j=1}^{N} \sum_{m=1}^{3} \sum_{l=1}^{2} R_{\theta_j} I_{a,mjl} . \qquad (D.6) $$