Issue 
A&A
Volume 653, September 2021



Article Number  A138  
Number of page(s)  22  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202039618  
Published online  24 September 2021 
RHAPSODIE: Reconstruction of HighcontrAst Polarized SOurces and Deconvolution for cIrcumstellar Environments
^{1}
Univ Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon, UMR5574, 69230 SaintGenisLaval, France
email: maud.langlois@univlyon1.fr
^{2}
Univ. Lyon, ENS de Lyon, Univ. Claude Bernard Lyon 1, CNRS, Laboratoire de Physique, 69342 Lyon, France
Received:
8
October
2020
Accepted:
11
May
2021
Context. Polarimetric imaging is one of the most effective techniques for the highcontrast imaging and characterization of circumstellar environments. These environments can be characterized through directimaging polarimetry at nearinfrared wavelengths. The SpectroPolarimetric Highcontrast Exoplanet REsearch (SPHERE)/IRDIS instrument, installed on the Very Large Telescope (VLT) in its dualbeam polarimetric imaging mode, offers the capability to acquire polarimetric images at high contrast and high angular resolution. However, dedicated image processing is needed to eliminate the contamination from the stellar light, instrumental polarization effects, and blurring from the instrumental point spread function.
Aims. We aim to reconstruct and deconvolve the nearinfrared polarization signal from circumstellar environments.
Methods. We used observations of these environments obtained with the highcontrast imaging infrared polarimeter SPHEREIRDIS at the VLT. We developed a new way to extract the polarimetric signal using an inverse approach method that benefits from the additional knowledge of the detected signal formation process. The method includes a weighted data fidelity term and smooth penalization, and it takes the instrumental polarization into account.
Results. This method enables us to accurately measure the polarized intensity and angle of linear polarization of circumstellar disks by taking into account the noise statistics and the convolution by the instrumental point spread function. It has the capacity to use incomplete polarimetry cycles, which enhance the sensitivity of the observations. The method improves the overall performances in particular for instances of both low signaltonoise (S/N) and small polarized flux compared to standard methods.
Conclusions. By increasing the sensitivity and including deconvolution, our method will allow for more accurate studies of these disks morphology, especially in the innermost regions. It also will enable more accurate measurements of the angle of linear polarization at low S/N, which would lead to indepth studies of dust properties. Finally, the method will enable more accurate measurements of the polarized intensity, which is critical for the construction of scattering phase functions.
Key words: methods: numerical / techniques: polarimetric / methods: data analysis / techniques: high angular resolution / protoplanetary disks / techniques: image processing
© L. Denneulin et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
With the adaptiveopticsfed highcontrast imaging instruments GPI (Macintosh et al. 2014) and SPHEREIRDIS (Beuzit et al. 2019; Dohlen et al. 2008), we now have access to the spatial resolution and sensitivity required to observe in the nearinfrared (NIR) circumstellar matter at small angular separations. Along with the Integral Field Spectrograph (IFS; Claudi et al. 2008) and the Zürich IMaging POLarimeter (ZIMPOL; Schmid et al. 2018) which can also be used to observe circumstellar environments in polarimetry, the IRDIS instrument is one of the three SPHERE instruments (Beuzit et al. 2019). SPHERE/IRDIS is able to acquire two simultaneous images at two different wavelengths, in a socalled Dual Band Imaging (DBI) mode (Vigan et al. 2014), or for two different polarizations, in a socalled Dual Polarimetry Imaging mode (Langlois et al. 2014; de Boer et al. 2020). Both circumstellar disks and selfluminous giant exoplanets or companion brown dwarfs can be characterized thanks to these new instruments in directimaging processes at these wavelengths.
The NIR polarimetric mode of SPHERE/IRDIS at the Very Large Telescope (VLT), which is described in Beuzit et al. (2019), de Boer et al. (2020), van Holstein et al. (2020), has been proven to be very successful in the detection of circumstellar disks in scattered light (Garufi, et al. 2017) and shows much promise for the characterization of brown dwarfs (van Holstein et al. 2017) and exoplanets (van Holstein et al. 2021, and in prep.) when they are surrounded by circumsubstellar disks.
There are three particular types of circumstellar disks that are the subject of recent studies: protoplanetary disks, transition disks, and debris disks. Observations of the protoplanetary and transition disks, in terms of their morphology and linked to hydrodynamical simulations, allow for the study of their formation scenario, as in studies of HD 142527 (Price et al. 2018), IM Lup, RU Lup (Avenhaus et al. 2018), and GSC 07396759 (Sissa et al. 2018). Their observations are valuable because their shapes can serve as the signposts for the formation of one or several exoplanets. In fact, during their formation, the planets “clean” the dust off their orbits, creating gaps without dust, as in the case of RXJ 1615, MY Lup, PDS 66 (Avenhaus et al. 2018), and PDS 70 (Keppler et al. 2018, 2019; Haffert et al. 2019). The planet formation scenario can be explained with hydrodynamical simulations as the cases of HL Tau (Dipierro et al. 2015), HD 163296 (Pinte et al. 2019). Due to gravity, exoplanets can also create spiral arms, as in RY Lup (Langlois et al. 2018) and MWC 758 (Benisty et al. 2015) where the presence of exoplanets is predicted by hydrodynamical simulations. Debris disks are the oldest step in the evolution of circumstellar disks, when there is already one or several planets in the system and the gas is almost completely consumed. These disks are composed of dust and grain that never accreted onto the planets, as in the case of HR 4796A (Perrin et al. 2015; Milli et al. 2019).
These environments can be observed with SPHERE in the nearinfrared and in the visible. However, such observations are difficult because of the high contrast between the light of the environment and the residual light from the host star. As a result, when acquiring images, the light of the environment is contaminated by the diffraction stains from the host star. Two methods can be used to disentangle the light of the disk from that of the star: angular differential imaging (ADI; Marois et al. 2006) and differential polarimetric imaging (DPI; van Holstein et al. 2020). The ADI technique makes use of the fact that the stellar residuals are fixed in the pupil plane and the object of interest artificially rotates. The resulting diversity makes it possible to disentangle the light of the object of interest from the residual light of the star. Still, such a method does not allow for a good reconstruction of the disk morphology as it is impacted by artifacts due to selfsubtraction. Moreover, the method fails when the environment is nearly rotationinvariant. The DPI observations allow for access to the morphology of the disks, without the artifacts, by using the difference of polarization states between the light scattered by the environment and the light of the host star.
The stateoftheart methods to process datasets in polarimetry, apart from the calibration, are “stepbystep” methods. First, the data are transformed with the required translations and rotations to be easier to process and the bad pixels are interpolated. Such interpolations introduce correlations that are not taken into account in the following processing. Second, the interpolated data are reduced to the Stokes parameters, which are directly related to the different polarization states. These reductions can be done with the double difference or the double ratio (Tinbergen 2005; Avenhaus et al. 2014). If the double ratio takes into account the possibility of multiplicative instrumental effects, none of these methods are equipped to deal with the noise statistics. This results in some limitations in sensitivity in the case of a low signaltonoise ratio (S/N). Lastly, a deconvolution may be performed to get rid of the blurring by the instrumental point spread function. As this is done without accounting for the noise statistics after all the earlier processing, the results are not optimal given the available data. Still, these stateoftheart methods have proven over the years to be sufficiently efficient to produce quality results. However, studies of circumstellar disks are often limited to analyses of the orientation (position angle and inclination) and morphology (rings, gaps, cavities, and spiral arms) of the disks (Muto et al. 2012; Quanz et al. 2013; Ginski et al. 2016; de Boer et al. 2016a). Quantitative polarimetric measurements of circumstellar disks and substellar companions are currently very challenging because existing datareduction methods do not properly estimate the sources of the errors from both noise and detector calibration. They also require complete polarimetric cycles and do not account for the instrument convolution. For observations of circumstellar disks (van Holstein et al. 2020), calibrating the instrumental polarization effects with a sufficiently high accuracy has already yielded several improvements.
Over the last decades in image processing, it has been proven that the reconstruction of parameters of interest benefits from a global “inverse problems” approach, taking into account the noise statistics and all instrumental effects, rather than “stepbystep” procedures. In this framework, image restoration methods rely on a physically grounded model of the data as a function of the parameters of interest and express the estimated parameters of interest as the constrained minimum of an objective function. This objective function is generally the sum of a data fidelity term and of regularization terms introduced to favor known priors. Depending on the convexity and on the smoothness of the objective function, several numerical algorithms with guarantees of convergence may be considered to seek the minimum (Nocedal & Wright 1999; Combettes & Pesquet 2011; Pustelnik et al. 2016). Such methods have been used over decades in astrophysics for the estimation of the physical parameters (Titterington 1985), mostly in adaptive optics (Borde & Traub 2006) as well as radiointerferometry with the wellknown algorithm CLEAN (Högbom 1974). This last method has been the starting point of a wide variety of algorithms, such as the algorithms SARA (Carrillo et al. 2012) and PolcaSARA (Birdi et al. 2020) in polarimetric radiointerferometry using more sophisticated tools as nonsmooth penalizations. A nonsmooth method was also used for images denoising with curvelets (Starck et al. 2003). The minimization of a cologlikelihood was also used in the blind deconvolution of images convolved by an unknown PSF with aberrations (Thiébaut & Conan 1995). Learning methods have been used more recently for the estimation of the CMB (Adam et al. 2016) and the imaging of the supermassive blackhole (Akiyama et al. 2019). In highcontrast imaging, the use of inverse problem methods is more recent. It has been used to perform autocalibration of the data (Berdeu et al. 2020) with the IFS/SPHERE. It has also been used to reconstruct extended objects in total intensity by using ADI data (Pairet et al. 2019; Flasseur et al. 2019) with the SPHERE instrument. However, such reconstruction methods have not been used in polarimetric highcontrast direct imaging.
In the present work, we describe in details the method and the benefits of the use of an inverse problem formalism for the reconstruction of circumstellar environments observed in polarimetry with the instrument ESO/VLT SPHERE IRDIS. In Sect. 2, we develop the physical model of the data obtained with the ESO/VLT SPHERE IRDIS instrument. This includes the polarimetric parametrization, the convolution by the PSF and the observing sequence. In Sect. 3, we describe Reconstruction of HighcontrAst Polarized Sources and Deconvolution for nearInfrared Environments (RHAPSODIE^{1}), the method we developed. Section 4 is dedicated to the calibration of the detector, the instrument, and the instrumental polarization. Finally, in Sect. 5, we present the results obtained with RHAPSODIE on both simulated and astrophysical data.
2. Modeling polarimetric data
The principal parameters of interest for studies of circumstellar environments in polarimetry are the intensity, I^{p}, of the linearly polarized light and the corresponding polarization angle θ which are caused by the reflection of the stellar light onto the circumstellar dust. By modulating the orientation of the instrumental polarization, these parameters can be disentangled from the intensity I^{u} of the unpolarized light received from the star and its environment. Without ADI observations, it is not possible to unravel the contributions by the star and by its environment from the unpolarized light, I^{u}.
The estimation of the parameters (I^{u}, I^{p}, θ) from polarimetric data is the objective of the present contribution. Stokes parameters are, however, more suitable in striving to account for the effects of the instrument on the observable polarization as the model of the data happens to be a simple linear combination of these parameters. Stokes parameters account for the total light, the linearly polarized light, and the circularly polarized light. Since circular polarization is mostly generated by magnetic interactions and double scattering, it is often negligible in the case of circumstellar environments and thus not measured by the SPHERE or GPI instruments. Ultimately, it is possible to reconstruct the parameters of interest I^{u}, I^{p}, and θ from a combination of the Stokes parameters.
2.1. Polarization effects
The four Stokes parameters S = (I, Q, U, V) ∈ ℝ^{4} describe the state of polarized light: I is the total intensity accounting for the polarized and unpolarized light, Q and U are the intensities of the light linearly polarized along 2 directions at 45° to each other, and V is the intensity of the circularly polarized light. Under this formalism, polarization effects by an instrument like SPHERE/IRDIS (see Fig. 1) can be modeled by (see Eq. (17) of van Holstein et al. 2020):
Fig. 1. Schematic view of the instrument ESO/VLTSPHERE IRDIS showing the various optical parts that can induce polarization effects. The same notations as in the text are used (e.g., n is the pixel index in the restored maps, k is the sequence index, and j is the polarizer index of the analyzer set). 
where are the Stokes parameters on the detector after the left (j = 1) or right (j = 2) polarizer of the analyzer set while S represents the Stokes parameters at the entrance of the telescope. In the above equation, M denotes a Mueller matrix accounting for the polarization effects of a specific part of the instrument: for the left or right polarizer of the analyzer set, M^{der} for the optical derotator, M^{HWP} is for the halfwave plate (HWP), M^{M4} for the fourth mirror of the telescope and M^{UT} for the three mirrors (M1 to M3) constituting the telescope. The term T(θ) denotes a rotation matrix of the polarization axes by an angle θ: θ^{der} is the derotator angle, α is the HWP angle, θ^{alt} is the altitude angle, and θ^{par} is the parallactic angle of the pointing of the altazimuthal telescope.
In the optical, only the total intensity out of the Stokes parameters can be measured by any existing detector. It then follows from Eq. (1) that the quantity measured by the detector is a simple linear combination of the input Stokes parameters:
where, for every ℓ ∈ ⟦1, 4⟧, S_{ℓ} denotes the ℓth component of the input Stokes parameters S and ν_{j, ℓ}(Θ), for j ∈ {1,2}, are real coefficients depending on all involved angles Θ = (θ^{der}, θ^{alt}, θ^{par}, α).
Even though the observables are restricted to the component , rotating the angle α of the HWP introduces a modulation of the contribution of the Stokes parameters Q and U in , which can be exploited to disentangle the Stokes parameters I, Q and U. The Stokes parameter V characterizing the circularly polarized light cannot be measured with an instrument such as SPHERE/IRDIS (a modulation by a quarterwave plate would have been required to do so). In the following, we therefore neglect the circularly polarized light and only consider the unpolarized and linearly polarized light characterized by the Stokes parameters S = (I, Q, U). As a direct simplification, the sum in the right hand side of Eq. (2) is reduced to its first three terms. For a sequence of acquisitions with different angles of the HWP, the detected intensities follow:
where ν_{j, k, ℓ} = ν_{j, ℓ}(Θ_{k}) with Θ_{k} the set of angles during the kth acquisition.
The instrumental polarization effects being reduced to those caused by the analyzers and the HWP and ignoring the rotation due to the altitude and parallactic angles, the detected intensity writes (van Holstein et al. 2020):
where ψ_{j} is the orientation angle of the left/right polarizer while α_{k} is the HWP angle during the kth acquisition. Table 1 lists the values of the linear coefficients ν_{j, k, ℓ} for a typical set of HWP angles.
Positions of HWP α_{k} and orientations of the analyzer ψ_{j} and the corresponding values of the coefficients ν_{j, k, ℓ} assuming no instrumental polarization and ignoring field rotation.
2.2. Parameters of interest
The model of the detected intensity given in Eq. (3) is linear with regard to the Stokes parameters S = (I, Q, U), which makes its formulation suited to inverse problem solving. However, to study circumstellar environments in polarimetry, the knowledge of the linearly polarized light I^{p} and the polarization angle θ is crucial. Both set of parameters are related as it follows:
and conversely by^{2}:
We assume that these relations hold independently at any position of the field of view (FOV). From Eqs. (3) and (5), the direct model of the detected intensity is a nonlinear function of the parameters I^{u}, I^{p}, and θ:
For a perfect SPHERE/IRDISlike instrument and not considering field rotation, combining Eqs. (4) and (5) yields:
which is the Malus law.
2.3. Accounting for the instrumental spatial PSF
In polarimetric imaging, each polarimetric parameter is a function of the twodimensional (2D) FOV. We consider that the Stokes parameters are represented by images of N pixels each and denote as S_{ℓ, n} the value of the nth pixel in the map of the ℓth Stokes parameter.
Provided that polarization effects apply uniformly across the FOV and that the instrumental spatial point spread function (PSF) does not depend on the polarization of light, all Stokes parameters of a spatially incoherent source are independently and identically affected by the spatial PSF (e.g., Birdi et al. 2018; Smirnov 2011; Denneulin 2020). As the effects of the spatial PSF are linear, we can give the detected intensity for a given detector pixel as:
where is the intensity measured during the kth acquisition by the mth pixel of the subimage corresponding to the jth polarizer of the analyzer set. The coefficients ν_{j, k, ℓ} (accounting for the instrumental polarization) are defined in Eq. (3) and H_{j, k, m, n} denotes a given entry of the discretized spatial PSF of the instrument. Here, j ∈ {1,2}, k ∈ ⟦1, K⟧ and m ∈ ⟦1, M⟧. Table 2 summarizes the main notations used in this paper.
Notations.
It follows from our assumption on spatial and polarization effects applied independently that the discretized spatial PSF does not depend on the polarization index ℓ. Consequently, the spatial and polarization effects in Eq. (9) mutually commute.
Figure 1 and Eq. (1) provide a representation of the instrument from which we build a model of the spatial effects of the instrument. Accordingly, an image representing the spatial distribution of the light as the input of the instrument should undergo a succession of image transformations before reaching the detector. These transformations are either geometrical transformations (e.g., the rotation depending on the parallactic angle) or blurring transformations (e.g., by the telescope). Except in the neighborhood of the coronographic mask, the effects of the blur can be assumed to be shiftinvariant and can thus be modeled by convolution with a shiftinvariant PSF. Geometrical transforms and convolutions do not commute but their order may be changed provided the shiftinvariant PSFs are appropriately rotated or shifted. Thanks to this property, and without loss of generality, we can model the spatial effects of the instrument by a single convolution accounting for all shiftinvariant blurs followed by a single geometrical transform accounting for all the rotations but also possible geometrical translations or spatial (de)magnification^{3}. Following this analysis, our model of the spatial PSF is given by:
where A_{k} : ℝ^{N} → ℝ^{N} implements the shiftinvariant blur of the input model maps while T_{j, k} : ℝ^{N} → ℝ^{Mj} performs the geometrical transform of the blurred model maps for the jth polarizer of the analyzer set during the kth acquisition. N is the number of pixels in the model maps and M_{j} is the number of pixels of the detector (or subimage) corresponding to the jth output polarizer.
In our implementation of the PSF model, the geometrical transform of images is performed using interpolation by CatmullRom splines. The blurring due to the instrument and the turbulence is applied by:
where F denotes a fast Fourier transform (FFT) operator of suitable size and implements the frequencywise multiplication by , the discrete Fourier transform of p_{k}, the shiftinvariant PSF. We note that p_{k} must be specified in the same reference frame as the FOV. If the shiftinvariant PSF is calibrated from empirical images of, for instance, the host star, acquired by the detector, the inverse (or pseudoinverse) of T_{j, k} must be applied to the empirical images.
2.4. Polarimetric data
During a sequence of observations with SPHERE/IRDIS in DPI mode, the HWP is rotated several times along a given cycle of angles α ∈ {0° ,45° ,22.5° ,67.5° }. In addition, the two polarizers of the analyzer set of SPHERE/IRDIS are imaged on two disjoint parts of the same detector. The resulting dataset consists in K frames composed of two (left and right) subimages, each with a different position of the HWP. Typical values of K can go from 32 to more than 512 depending of the observed target. After preprocessing of the raw images to compensate for the bias and the uneven sensitivity of the detector and to extract the two subimages, the available data are modeled by:
where is given in Eq. (9) for k ∈ ⟦1, K⟧ the index of the acquisition, j ∈ {1, 2} indicating the left and right polarizer of the analyzer, and m ∈ ⟦1, M⟧ the pixel index in the corresponding left and right subimage.
The ≈ sign in Eq. (12) is to account for an unknown random perturbation term due to the noise. Noise in the preprocessed images can be assumed to be centered and independent between two different pixels or frames because the preprocessing suppresses the bias and treats pixels separately thus introducing no statistical correlations between pixels. There are many sources of noise: shot noise for the light sources and the dark current, detector readout noise, etc. For most of the actual data, the shot noise is the most important contribution and the number of electrons (photoelectrons plus dark current) integrated by a pixel is large enough to approximate the statistics of the data by an independent nonuniform Gaussian distribution whose mean is given by the righthandside term of Eq. (12) and whose variance Σ_{j, k, m} = Var(d_{j, k, m}) is estimated in the calibration stage (see Sect. 4.1).
It is worth noting that since the proposed model is not valid in the neighborhood of the coronographic mask, the data pixels in this region have to be discarded. Moreover, the detector contains defective pixels (e.g., dead pixels, nonlinear pixels, saturated pixels) which must be also discarded. This is achieved by assuming that their variance is infinite, which amounts to setting their respective weights to zero in the data fidelity term of the objective function described in Sect. 3.2.
3. RHAPSODIE : Reconstruction of HighcontrAst Polarized SOurces and Deconvolution for cIrcumstellar Environments
3.1. Inverse problems approach
In polarimetric imaging, one is interested in recovering sampled maps of the polarimetric parameters, which can be the Stokes parameters (I, Q, U) or the intensities of unpolarized and linearly polarized light and the angle of the linear polarization (I^{u}, I^{p}, θ) or some mixture of these parameters. In order to remain as general as possible, we denote by X ∈ ℝ^{N × L} the set of parameters of interest to be recovered and by X_{ℓ} ∈ ℝ^{N}, with ℓ ∈ ⟦1, L⟧, the ℓth parameter component which is a Npixel map.
Given the direct model of the preprocessed data developed in the previous section, we propose to recover the parameters of interest X by a penalized maximum likelihood approach. This approach is customary in the solving of inverse problems (Titterington 1985; Tarantola 2005) and amounts to defining the estimated parameters as the ones that minimize a given objective function f(X) possibly under constraints expressed as X ∈ 𝒞 with 𝒞 the set of acceptable solutions. The objective function takes the form of the sum of a datafidelity term f_{data}(X) and of regularization terms f_{ρ}(X):
where λ_{ρ} ≥ 0 (∀ρ) are socalled hyperparameters introduced to tune the relative importance of the regularization terms. The datafidelity term f_{data}(X) imposes the direct model be as close as possible to the acquired data while the regularization terms f_{ρ}(X) enforce the components of the model to remain regular (e.g., smooth). Regularization must be introduced to lift degeneracies and avoid artifacts caused by the data noise and the illconditioning of the inverse problem. Additional strict constraints may be imposed on the sought parameters via the feasible set 𝒞, e.g., to account for the requirement that intensities are nonnegative quantities. These different terms and constraints are detailed in the following subsections.
3.2. Data fidelity
With knowledge of the appropriate statistics for the preprocessed data, the agreement of the model with the data is properly insured by the cologlikelihood of the data (Tarantola 2005), or equivalently by the following criterion:
where denotes the Mahalanobis (1936) squared norm, d_{j, k} = (d_{j, k, 1}, …, d_{j, k, M})^{⊤} ∈ ℝ^{M} collects all the pixels (e.g., in lexicographic order) of the jth subimage in the kth acquisition as defined in Eq. (12). Similarly, where the terms are given by the model in Eq. (9) applied to the Stokes parameters as a function of the parameters of interest X, S = S(X). For instance, if , then S(X) is obtained by Eq. (5). In the expression of the data fidelity term given by Eq. (14), W_{j, k} is the precision matrix of the data. The precision matrix is diagonal because pixels are considered as mutually independent. To account for the nonuniform noise variance and for invalid data (see Sect. 2.4), we define the diagonal entries of the precision matrix as:
with Σ_{j, k, m} = Var(d_{j, k, m}) the variance of a valid datum d_{j, k, m}. Invalid data include dead pixels, pixels incorrectly modeled by our direct model because of saturation or of the coronograph, missing frames for a given HWP angle, and unusable frames due to overly strong atmospheric turbulence or improper coronograph centering.
We note that in Eq. (14), the Mahalanobis squared norms arise from our Gaussian approximation of the statistics, while the simple sum of these squared norms for each subimage in each frame is justified by the fact that all frames and all subimages are mutually independent.
3.3. Regularization
The problem of recovering the polarimetric parameters from the data is an illconditioned inverse problem mainly due to the instrumental blur. Furthermore, the problem may also be illposed if there are too many invalid data. In the case of an illconditioned inverse problem, the maximum likelihood estimator of the parameters of interest, that is the parameters which minimize the data fidelity term f_{data}(X) defined in Eq. (14) alone, cannot be used because it is too heavily corrupted by noise amplification. Explicitly requiring that the sought parameters be somewhat regular is mandatory to avoid this (Titterington 1985; Tarantola 2005). In practice, this amounts to adding one or more regularization terms f_{ρ}(X) to the datafidelity as assumed by the objective function defined in Eq. (13).
3.3.1. Edgepreserving smoothness
We expect that the light distribution of circumstellar environments be mostly smooth with some sharp edges, hence, using an edgepreserving smoothness regularization (Charbonnier et al. 1997) appears to be the most suited choice to this kind of light distribution. When considering the recovering of polarimetric parameters, such a constraint can be directly imposed to the unpolarized intensity I^{u}, to the polarized intensity I^{p} or to the total intensity I by the following regularization terms:
where we denote in boldface sampled maps of polarimetric parameters, for instance the image of the unpolarized intensity or I^{u}(X) to make explicit that it is uniquely determined by the sought parameters X. In the above expressions, μ_{ρ} > 0 models the smoothing threshold and D_{n} : ℝ^{N} → ℝ^{2} is a linear operator which yields an approximation of the 2D spatial gradient of its argument around the nth pixel. This operator is implemented by means of finite differences; more specifically applying D_{n} to a sampled map u of a given parameter is expressed as:
where (n_{1}, n_{2}) denotes the row and column indices of the nth pixel in the map. At the edges of the support of the parameter maps, we simply assume flat boundary conditions and set the spatial gradient to zero there.
The regularizations in Eqs. (16)–(18) implement a hyperbolic version of a pseudonorm of the spatial gradient of a given component of the light distribution which behaves as an L_{2}norm (i.e., quadratically) for gradients much smaller than μ and as an L_{1}norm (i.e., linearly) for gradients much larger than μ. Hence imposing smoothness for flat areas where the spatial gradient is small while avoiding strong penalization for larger spatial gradients at edges of structures.
It has been shown (Lefkimmiatis et al. 2013; Chierchia et al. 2014) that grouping different sets of parameters in regularization terms that are subL_{2} norm of the gradient, as with the last one in Eqs. (16)–(18), yields solutions in which strong changes tend to occur at the same locations in the sets of parameters. In order to encourage sharp edges to occur at the same places in the Stokes parameters Q and U, we also consider using the following regularization for these components:
Many other regularizations implementing smoothness constraints can be found in the literature from the simple quadratic one (Tikhonov 1963) to the very popular total variation (TV; Rudin et al. 1992). Quadratic regularizations tend to yield strong ripples which, owing to the contrast of the recovered maps, are an unacceptable nuisance while TV yields maps affected by a socalled cartoon effect (i.e., piecewise flat images), which is not appropriate for realistic astronomical images. We however note that the hyperbolic edgepreserving regularization with a threshold, μ_{ρ}, set to a very small level can be seen as a relaxed version of TV and has been widely used as a differentiable approximation of this regularization. Our choice of a differentiable regularization is also motivated by the existence of efficient numerical methods to minimize nonquadratic but differentiable objective functions of many (millions or even billions) variables possibly with additional strict constraints (Thiébaut 2002). We refer to Denneulin et al. (2019, 2020), Denneulin (2020) for a comparison of possible advanced regularizers.
3.3.2. Tuning of the hyperparameters
In the regularization function f_{prior}(X), the terms defined in Eqs. (16)–(18) and Eq. (20) can be activated (or inhibited) by choosing the corresponding λ_{ρ} > 0 (or λ_{ρ} = 0). It is also required to tune the threshold level μ_{ρ} > 0 in addition to the λ_{ρ} multipliers. All these hyperparameters have an incidence on the recovered solution: the higher λ_{ρ} the smoother the corresponding regularized component and a lowering of the threshold, μ_{ρ}, allows us to capture sharper structures. A number of practical methods have been devised to automatically tune the hyperparameters: Stein’s Unbiased Risk Estimator (SURE: Stein 1981; Eldar 2008; Ramani et al. 2008; Deledalle et al. 2014), generalized crossvalidation (GCV: Golub et al. 1979), the Lcurve (Hansen & O’Leary 1993), and hierarchical Bayesian strategies (Molina 1994).
Unsupervised tuning of the hyperparameters with GSURE (Eldar 2008), in a prediction error formulation, has been considered in this context. Figure 2 shows the value of GSURE for several reconstructions of RXJ 1615 (Avenhaus et al. 2018). Since it is expected that there are relatively few sharp edges, ∥D_{n}X_{ρ}∥≪μ_{ρ} should hold for most pixels n in the component X_{ρ}. As a result, for most pixels n, the regularization penalty behaves as a quadratic Tikhonov (1963) smoothness weighted by λ_{ρ}/μ_{ρ}:
Fig. 2. Comparison of the nonlinear RHAPSODIE reconstructions of RXJ 1615 (Avenhaus et al. 2018) for different values of λ_{Q + U} and μ_{Q + U}. The value of the GSURE criterion is indicated for each reconstruction. 
The strength of the blurring imposed by the regularization is therefore mostly controlled by the value of λ_{ρ}/μ_{ρ} (e.g., from top to bottom in Fig. 2) while the sharp edges are controlled by the threshold μ_{ρ} (e.g., from left to right in Fig. 2).
In Fig. 2, we can observe that an automatic selection of the hyperparameters with GSURE would lead to an over regularized solution. In other highcontrast contexts, we also observed this tendency of GSURE to oversmooth the result. We think that devising a good unsupervised method for tuning the hyperparameters deserve further stud; in addition, for the results presented in this paper, we tuned the hyperparameters by hand by visually inspecting several reconstructions under different settings, as presented in Fig. 2.
The polarized parameter I^{p} contribute to I, Q, and U (see Eq. (5)). To avoid a contamination of the extracted polarized parameters I^{p} by the unpolarized component I^{u}, we regularize I^{u} and I^{p} separately. Hence, the regularization of the unpolarized component I^{u} should be done via defined in Eq. (16) rather than via f_{I} defined in Eq. (18). For the polarized light, the joint regularization of Q and U by f_{Q + U} defined in Eq. (20) is more effective than the regularization of I^{p} alone by defined in Eq. (17) which does not constrain the angle θ of the linear polarization. In Eq. (13), we therefore take the Set 1 or Set 2 of hyperparameters presented in Table 3. The latter combination is preferable as discussed previously.
Set of hyperparameters used in the present work.
3.4. Imposing the positivity of the intensities
Imposing nonnegativity constraints on the restored intensities has proven its efficiency for astronomical imaging where large parts of the images consist in background pixels whose value should be zero (Biraud 1969). Whatever the choice of the parametrization, the intensities I, I^{u}, and I^{p} should all be everywhere nonnegative.
Since I = I^{u} + I^{p}, it is sufficient to require that I^{u} and I^{p} be nonnegative. Hence, for the set of parameters , the positivity constraint is expressed as:
As given for the Stokes parameters X = (I, Q, U), the positivity yields an epigraphical constraint:
Such a constraint can be found in Birdi et al. (2018), but it has not yet been implemented in high contrast polarimetric imaging.
Since (for all pixels n), the positivity of the intensity I^{p} of the polarized light automatically holds if the parameters X = (I^{u}, Q, U) are considered. It is then sufficient to impose the positivity of the intensity I^{u} of the unpolarized light as expressed by the following feasible set:
3.5. Choice of the polarimetric parameters
Our method expresses the recovered parameters as the solution of a constrained optimization problem specified in Eq. (13). As explained in Sect. 3.3.2, the weights of imposed regularization is chosen via the values of the multipliers λ_{ρ}.
The constraints can be implemented by the feasible 𝒞 for different choices of the parameters X. More specifically, X = (I^{u}, I^{p}. θ), X = (I, Q, U), or X = (I^{u}, Q, U) can be chosen. Whatever the choice for X, the relations given in Eqs. (5) and (6) can be used to estimate any parameter of interest given the recovered . These relations can also be used to compute the objective function f(X) which require the Stokes parameters needed by the direct model of the data Eq. (9) and various polarimetric component depending on the choice of the regularization.
With X = (I, Q, U), the positivity constraints take the form of an epigraphic constraint that is more difficult to enforce as being not separable in the parameters space. To solve the problem in Eq. (13) with such a constraint, an epigraphic projection is required, leading to the use of a forwardbackward scheme (Combettes & Wajs 2005), reduced in this context to a standard projected gradient descent. A description of the method for such a minimization problem can be found in Denneulin et al. (2020).
The choice of X = (I, Q, U) may avoid some degeneracies because it ensures the convexity of the problem in Eq. (13). With X = (I^{u}, I^{p}, θ) or X = (I^{u}, Q, U), the positivity constraints amounts to applying simple separable bound constraints on some parameters. Since the objective is differentiable, a method such as VMLMB (Thiébaut 2002) can be used to solve the problem in Eq. (13). The VMLMB algorithm is a quasiNewton method with limited memory requirements and able to account for separable bound constraints. VMLMB is applicable to large size problems and only requires the provision of the bounds and a numerical function to compute the objective function and its gradient.
In the following, we compare the performances of RHAPSODIE for the polarimetric parameters X = (I, Q, U) and X = (I^{u}, Q, U) on simulated synthetic datasets. We refer to linear RHAPSODIE when we reconstruct X = (I, Q, U) and to nonlinear RHAPSODIE when we reconstruct X = (I^{u}, Q, U). Both parametrizations allow us to access the parameters of interest I^{p} and θ using Eq. (6). The best choice of polarimetric parameters is then used to process astrophysical datasets.
4. Data calibrations
4.1. Detector calibration
Before the application of a reconstruction method, the calibration of the data is essential to account for the noise and the artifacts linked to the measurement. It allows for the estimation and correction of any pollution induced by the sky background or the instrument as well as the detector behavior in terms of errors on pixel values.
We use an inverse method to calibrate the raw data from these effects: the quantity required for the calibration are jointly estimated from the likelihood of the calibration data direct model (Denneulin 2020). In this method, all the calibration data are expressed as a function of the different contributions (ı.e., flux, sky background, instrumental background, gain, noise, and quantum efficiency). All these quantities are then jointly estimated by the minimization of the cologlikelihood of the data. Calibrated data are corrected for contribution of the bias and the background and for nonuniform sensitivity and throughput. The calibration also provides associated weights computed according to the estimated variance, see Eq. (15). Finally, the defective pixels are detected by crossing several criteria, such as their linearity, their covariance compared to that of the other pixels, or the values of their likelihood in the calibration data. This calibration method produces calibrated data outputs (d_{k})_{k ∈ ⟦1, K⟧} and their weights (W_{k})_{k ∈ ⟦1, K⟧}.
4.2. Instrumental calibration
The instrumental calibration is a required step for the estimation, from dedicated data, of the instrumental PSF and of the star centers on each side of the detector. Since the star is placed behind the coronagraphic mask, a simultaneous PSF estimation is not possible. To estimate the PSF, we use a dedicated flux calibration (STARFLUX) that is acquired just before and after the science exposure by offsetting the telescope to about 0.5 arcsec with respect to the coronagraphic mask by using the SPHERE tip/tilt mirror (Beuzit et al. 2019). Consequently, the PSF is recorded with similar atmospheric conditions (listed in Table 4) as the science observations. When performing this calibration, suitable neutral density filters are inserted to avoid detector saturation. It has been shown in Beuzit et al. (2019) that these neutral density filters do not affect the PSF shape and thus its calibration. This instrumental calibration leads to the estimation of the PSF modeled through the operators A_{k}. The PSF model does not include the spiders so that it may remain rotationinvariant (see Fig. 3).
Fig. 3. Example of the fitting of the true PSF on the empirical SPHERE/IRDIS data of the target HD 106906 reduced with RHAPSODIE. Observed PSF in H band (upper left) and different parametrization (Airy, Gaussian, and Moffat). Intensities are in logscale to enhance the faint diffraction patterns. 
Information of the datasets used for the target reconstruction.
In the case of synthetic observations, the assumed 2D PSF is extracted from the real data RY Lup (e.g., similar to the top left image in Fig. 3). For astrophysical observations, we fit a 2D PSF model on prereduced PSF data for each observed target (e.g., Gaussian, Airy, and Moffat fits in Fig. 3).
During the coronographic observing sequence, the star point spread function peak is hidden by the coronagraphic mask and its position was determined using a special calibration (STARCENTER) where four faint replicas of the star image are created by giving a bidimensional sinusoidal profile to the deformable mirror (see Beuzit et al. 2019). The STARCENTER calibration was repeated before and after each science observation, and the resulting center estimations were averaged. In addition, we used the derotator position and the true north calibration from Maire et al. (2016) to extract the angle of rotation of the north axis. These instrumental calibration steps lead to the estimation of the transformation operators, T_{j, k}, which rotate and translate the maps of interest to make the centers and the north axes coincide with those in the data.
4.3. Polarization calibration
When the light is reflected by the optical devices in the instrument, some instrumental polarization is introduced, resulting in a loss of polarized intensity and crosstalk contamination between the Stokes parameters Q and U.
The classical method to compensate for this instrumental polarization is to employs the azimuthal Q_{ϕ} and U_{ϕ} parameters to reduce the noise floor in the image (Avenhaus et al. 2014). Because this method is limited to faceon disks, we instead use the method developed by van Holstein et al. (2020) which rely on the precomputed calibration of the instrumentation polarization as a function of the observational configurations. The pipeline IRDAP (van Holstein et al. 2020) yields the possibility to determine and correct the instrumental polarization in the signal reconstruction, yet this reconstruction requires to estimate first the Q and U parameters and then to perform the instrumental polarization correction. After computing the double difference, IRDAP uses a Mueller matrix model of the instrument carefully calibrated using real on sky data to correct for the polarized intensity I^{p} (created upstream of the HWP) and crosstalk of the telescope and instrument to compute the modelcorrected Q and U images.
Using IRDAP, we estimate the instrumental transmission parameters (ν_{j, k, l})_{j ∈ {1, 2},k ∈ ⟦1, K⟧, ℓ ∈ ⟦1, L⟧} from the Mueller matrix to calibrate the instrumental polarization in our datasets.
IRDAP also determines the corresponding uncertainty by measuring the stellar polarization for each HWP cycle individually and by computing the standard error of the mean over the measurements. Finally, IRDAP creates an additional set of Q and U images by subtracting the measured stellar polarization from the modelcorrected Q and U images. We use a similar method to correct for the stellar polarization which is responsible for strong light pollution at small separation in particular for faint disks, such as debris disks. The contribution of this stellar light is estimated as different factors of I^{u} in both Q and U, respectively, ε_{Q} and ε_{U}. We estimate this stellar contribution by using a pixel annulus Ω located at a separation where there is no disk signal (either near the edge of the coronograph or at the separation corresponding to the adaptive optics cutoff frequency). Both correction factors ε_{Q} and ε_{U} are expressed as follows:
where N_{Ω} is the number of pixels of the ring Ω. We then compute Q^{cor} = Q − ε_{Q}I^{u} and U^{cor} = U − ε_{U}I^{u} in order to create an additional set of Q^{cor} and U^{cor} images by subtracting the measured stellar polarization from the modelcorrected Q and U images.
5. Applications on highcontrast polarimetric data
In this section, we compare two RHAPSODIE configurations. The first configuration is defined as “without deconvolution”, meaning that in the data fidelity term (14) does not include the convolution, leading to the simplification of the Eqs. (9) and (10):
(i.e., A_{k} is the identity). Such a configuration of RHAPSODIE aims to be comparable to stateoftheart methods (i.e., Double Difference and IRDAP). These methods do not include any deconvolution and we aim to show that RHAPSODIE also performs well for such a case.
The second configuration is defined as “with deconvolution” and states for the full RHAPSODIE capabilities. Such a configuration of RHAPSODIE is aimed at showing the benefits of the global model compared to an a posteriori deconvolution with the aim of improving the angular resolution.
5.1. Application on synthetic data
The performance of the linear and nonlinear methods are first evaluated on synthetic datasets, without (resp. with) the deconvolution displayed in Figs. 4–6 (resp. Figs. 7–9). These datasets are composed of unpolarized residual stellar flux, mixed with unpolarized and polarized disk flux. We produce several synthetic datasets following steps given in Appendix A, for different ratios of the polarization of the disk, called τ^{disk} and defined in Eq. (A.1).
Fig. 4. Visual comparison of the reconstructed polarized intensity I^{p} with the stateoftheart double difference and the RHAPSODIE methods without deconvolution. 
Fig. 5. Maps of errors of the reconstructions displayed in Fig. 4. These errors are obtained as the difference between the true and the reconstructed images. 
Fig. 6. Comparison of the MSE between the true map and the estimated map of the polarized intensity I^{p} and of the angle of polarization θ for the double difference and the linear and nonlinear RHAPSODIE methods without deconvolution. 
Fig. 7. Visual comparison of the reconstructed polarized intensity I^{p} with the stateoftheart double difference and the RHAPSODIE methods with deconvolution. 
Fig. 8. Maps of errors of the reconstructions displayed in Fig. 7. These errors are obtained as the difference between the true and the reconstructed images. 
Fig. 9. Comparison of the MSE between the true map and the estimated map of the polarized intensity I^{p} and of the angle of polarization θ for the double difference and the linear and nonlinear RHAPSODIE methods with deconvolution. 
The results of the RHAPSODIE methods are compared to the results obtained with the classical double difference method (Tinbergen 2005; Avenhaus et al. 2014). The double difference is applied on recentered and rotated datasets with the bad pixels interpolated. For the comparison with deconvolution, the results of the Double Difference are deconvolved after the reduction. In order to provide fair comparisons with RHAPSODIE, we propose to use an inverse approach rather than applying a highpass spatial filtering. The deconvolved double difference reconstructions are obtained by solving the following problem:
where and are the Stokes parameters reconstructed with the double difference and A represents the convolution by the PSF. This deconvolution method performs the deconvolution of Stokes parameters Q and U jointly in order to sharpen the polarized intensity image, but does not recover the polarization signal lost by averaging over closeby polarization signals with the opposite sign and may even introduce artificial structures that are not present in the original source.
For the reconstruction with RHAPSODIE, the hyperparameters of regularization, λ_{ρ} and μ_{ρ}, are chosen in order to minimize the total mean square error (MSE), i.e., the sum of the MSE between each estimated parameter , and , obtained from from Eq. (6), and the ground truth , and θ_{gt}. The total MSE is given by:
where, , , and N_{θ} are the number of pixels with a signal of interest in the respective I^{u}, I^{p}, and θ maps. For the deconvolution of the double difference results, and are chosen to minimize only the sum of the MSE on I^{u} and θ.
For the reconstructions without the deconvolution, Figs. 4 and 5 show that the reconstructions are less noisy with RHAPSODIE. The inner circle is always better reconstructed. Moreover, the nonlinear reconstruction (i.e., minimization on I^{u}, Q and U) is better than the linear when τ^{disk} grows. In both configurations of RHAPSODIE (nonlinear and linear), the thin ring is not as well reconstructed as with the double difference. It is possible to recover such sharper structure with RHAPSODIE, by reducing the regularization weight (i.e., reducing the hyperparameter λ_{Q + U}). However, if λ_{Q + U} is too small, the data will be overfitted and the noise in the reconstructed images will be amplified. It is thus necessary to keep a good tradeoff between a smooth solution and a solution close to the data. Classically, minimizing the MSE is a good tradeoff between underfitting and overfitting. The MSE curves in Fig. 6 are coherent with the observations. When τ^{disk} = 25%, the nonlinear RHAPSODIE MSE and double difference MSE are equivalent, but on the reconstructions, we can see that if the thin circle is better reconstructed with the double difference, the inner circle is better reconstructed with RHAPSODIE. Moreover, we can see that the angle error is more than threetime smaller with RHAPSODIE compared to the double difference error.
For reconstructions with the deconvolution, the RHAPSODIE methods have a better approach to the reconstructions of the polarized intensity, as shown on the reconstructions in Fig. 7, mostly for the inner ring, even if the background is less biased with the proposed deconvolved double difference (as shown on the error maps in Fig. 8). In fact, as shown on the MSE displayed in Fig. 9, for τ^{disk} = 3%, RHAPSODIE delivers a more accurate polarized intensity estimation. The errors of the reconstructions are larger for the double difference than for the RHAPSODIE methods. The RHAPSODIE methods also allow us to achieve a better reconstruction of the angle of polarization, mostly with the nonlinear estimation. The nonlinear reconstruction appears to be more efficient as being not polluted by the artifacts of deconvolution of the unpolarized point source companion as seen in Fig. 8. Even if according to the MSE in Fig. 9, for τ^{disk} = 25%, the MSE is smaller for the linear RAPSODIE. in fact, the space between the outer ring and the thin ring is better reconstructed with such a configuration (see Figs. 7 and 8).
According to these results, the RHAPSODIE methods are better than the stateoftheart methods, in particular the nonlinear RHAPSODIE method. The benefits of our contribution is clearly visible in the case of faint disks and structures, which are the most common on astrophysical data. This is why in the following section dedicated to the astrophysical data, we select the nonlinear RHAPSODIE method with manual selection of the hyperparameters.
5.2. Astrophysical data
RHAPSODIE was applied to several IRDIS datasets dedicated to protoplanetary, transition, and debris disks to test the efficiency of our method and to compare it to the stateoftheart method. For each reconstruction, we present the map of the projected intensities, by using the standard azimuthal Stokes parameters Q_{ϕ} and U_{ϕ} estimated from:
with where (n_{1},n_{2}) denotes the row and the column indices of the nth pixel in the map and those of the pixel center. We also present the maps of polarized intensities I^{p} and the angles of polarization θ estimated by the different methods.
First, the reconstructions of the target TW Hydrae with the double difference and RHAPSODIE without deconvolution are compared in Fig. 10, without and with the correction of the instrumental polarization. The images are scaled by the square of the separation to account for the dropoff of stellar illumination with distance.
Fig. 10. Reconstruction of the Q (a) and U (b) parameters of the target TW Hydrae, for the first three cycle of HWP rotation, without (upper row) and with (lower row) the correction of the polarization. Without the correction, both Q and U are rotated and attenuated. The Q_{ϕ} and U_{ϕ} images reconstructed from the entire dataset are presented in (c). All the reconstructions are done without deconvolution to demonstrate mainly the efficiency of the instrumental polarization correction and the benefits of RHAPSODIE. (c) Azimutal Stokes parameters displayed in arcseconds: (i) Double Difference, (ii) Double Difference with the instrumental polarization correction, (iii) RHAPSODIE which includes by default the instrumental polarization correction. The intensities are multiplied in each pixel by the distance to the star r^{2}. 
The instrumental polarization in this dataset introduces a polarization rotation and an attenuation (loss of polarization signal) of the intensity with varying time during the observations, which can be monitored easily because the disk is faceon. If uncorrected, the combination of data from multiple polarimetric cycles will result in very poorly constrained polarimetric intensity measurements. When we correct the instrumental polarization (see Fig. 10c) by using the IRDAP method, these effects are compensated and the disk reconstructed is more accurate (Fig. 10c (iii)). As a result, the contamination of the U_{ϕ} signal from crosstalk is decreased and becomes negligible as seen in (Fig. 10c (iii)).
The comparison with van Boekel et al. (2017), de Boer et al. (2020) shows that our method is less impacted by the bad pixels and improves the disk S/N in the area where the signal is low. On the other hand, it is slightly more sensitive to detector flat calibration accuracy. It is worth noticing that IRDAP has a dedicated correction of the detector response nonuniformity (flat field variation between the detector column) for the various amplifiers that is not implemented in RHAPSODIE. However, the flat calibrations of observations more recent than 2016 have improved and do allow for a better calibration and, as a consequence, they do not impact our method efficiency anymore.
The efficiency of the method we have developed is demonstrated on the Figs. 11–17. Figures 11 and 12 present the double difference and RHAPSODIE reconstructions of the protoplanery disks TW Hydrae (van Boekel et al. 2017), IM Lupus (Avenhaus et al. 2018) and MY Lupus (Avenhaus et al. 2018). Figures 13, 14, and Fig. 17 present RHAPSODIE reconstructions of transition disks RY Lupus (Langlois et al. 2018) and T Chamaeleontis (Pohl et al. 2017). Figures 13 and 14 also present RHAPSODIE reconstructions of transition disk RXJ 1615 (de Boer et al. 2016a). Figures 15 and 16 present the RHAPSODIE reconstructions of the debris disks HD 106906 (Kalas et al. 2015; Lagrange et al. 2016), HD 61005 (Olofsson et al. 2016), and AU Mic (Boccaletti et al. 2018). These reconstructions are all normalized in contrast to the unpolarized stellar flux (estimated when registering the PSF offcentered from the coronograph) by taking into account the transmission of the neutraldensity filter used when registering the PSF to prevent saturation. The hyperparameters of regularization and the PSF models used for the reconstructions are listed in Table 5.
Fig. 11. Reconstructions of the polarized intensity I^{p} of the protoplanetary disks TW Hydrae, IM Lupus, and MY Lupus. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 
Fig. 12. Reconstructions of the polarization angle θ of the protoplanetary disks TW Hydrae, IM Lupus, and MY Lupus. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 
Fig. 13. Reconstructions of the polarized intensity I^{p} of the transition disks RY Lupus, T Chae, and RXJ 1615. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 
Fig. 14. Reconstructions of the polarization angle θ of the transition disks RY Lupus, T Chae, and RXJ 1615. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 
Fig. 15. Reconstructions of the polarized intensity I^{p} of the debris disks HD 106906, HD 61005, and AU MIC. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 
Fig. 16. Reconstructions of the polarized intensity I^{p} of the debris disks HD 106906, HD 61005, and AU MIC. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIEwithout deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 
Fig. 17. Q_{ϕ} and U_{ϕ} projections of the reconstructions of RY Lupus and T Chae. The maps are displayed in symmetrical logarithmic scale and normalized in contrast to the unpolarized stellar flux. The coronograph is masked in white. 
Information related to the parameters used for the target reconstruction presented in this section.
The contrast between the selected disks polarized light and their host stars unpolarized light responsible for the level of stellar pollution are also very different (ranging from 3 × 10^{−2} to 8 × 10^{−6}). The achievable contrast is more favorable for highly inclined disk, such as MY Lup in Fig. 11, because in such a case, the star likely shines partially through the disk, which is dimming the starlight and thus decreasing the contrast between the star and the disk.
In all cases, our method produces highquality reconstructions of the disk polarized signal and minimizes the artifacts from bad pixels. The comparison of these reductions with the stateoftheart methods confirms the benefits of our method, but it is more difficult to quantify the gain on real data without knowing the reality; and, thus, this approach should be strengthened through the use of numerical models for these objects. In addition, our method can produce deblurred results, which clearly enhances the angular resolution and is beneficial for interpreting the disk morphology and for studying its physical properties. In all cases, the deconvolution sharpens the polarized intensity image, helps to recover the polarization signal lost by averaging over closeby polarization signals with opposite sign, and does not introduce artificial structures that are not present in the original source (see Fig. 13 (RY Lup)). As mentioned earlier in this paper, the hyperparameters could be further tuned to adjust the smoothing of the deconvolved images according to the required noiseandangular resolution tradeoff.
For several of these disks (IM Lup, MY Lup in Fig. 11), the outer edge of the disk and thus the lower disk surface have been detected in Avenhaus et al. (2018) and are confirmed by our method. The deconvolution of these datasets allows to further highlight these fainter features, to enhance or reveal the midplane gaps in the case of T Cha (in Fig. 13) which was not identified by Pohl et al. (2017). The reason for this is that without deconvolution, the PSF smears light from the disk upper and lower sides into the midplane gap.
Except for Au Mic, the debris disks presented in Fig. 15 are much fainter in contrast than the protoplanetary or the transition disks in polarized intensity. As a consequence, these datasets have required careful stellar polarization compensation. The HD 106906 debris disk is viewed close to edgeon in polarized light as reported in van Holstein et al. (2020), Esposito et al. (2020). The image clearly shows the known eastwest brightness asymmetry of the disk, which was detected in total intensity (Kalas et al. 2015; Lagrange et al. 2016). Thanks to the deep dataset and good reconstruction, we also detect the backwardscattering far side of the disk to the west of the star, just south of the brighter near side of the disk. This feature is further highlighted by the deconvolution.
We may remark that the RHAPSODIE deblurred reconstructions of RXJ 1615 (in Fig. 13) and HD 106906 (in Fig. 15) appear noisier than those without deconvolution. This is due to the hyperparameters that were purposely chosen so as to achieve the best angular resolution with the side effect of slight noise amplification. Thus, the noise is not negligible in the reconstruction. As shown in Fig. 2, increasing the regularization weight to smooth the background would lead to a loss of the thin structures.
We also analyzed two datasets (HD 61005 and Au Mic), which were taken under very bad conditions. For these two datasets, our method is also proven to be efficient in using incomplete polarization cycles to recover the disk polarized signal and to deconvolve this signal despite strong artifacts additionally produced by the rotating spiders (i.e., unmasked by the Lyot stop) when observing in field stabilized. The difference in the spider position during the polarimetric cycle results in an artificial polarimetric signal when using standard data reduction techniques. It is worth noticing that the strength of our method to deal with these artifacts comes from its ability to use weighted maps to account for them. For instance, the spiders are weighted by their variance during their rotation frame by frame, in addition to the inclusion of a static bad pixel ponderation. As a result, the contribution of the spiders to the polarimetric signal when using our method is decreased compared to the other methods. The noise created by these spiders remains, as seen in Fig. 15, and can generate artifacts in the deconvolution as seen for HD 61005. Counteracting the artifacts caused by the spiders can be efficiently done by performing DPI observation in pupil tracking mode as proposed in van Holstein et al. (2017). We have also validated the efficiency of RHAPSODIE in such observations which allows a further gain in the reduction of the instrumental artifacts from the telescope spiders.
Another advantage of our method is its ability to use incomplete polarimetric cycles, which are discarded by the classical methods. This capability of our method leads to an increase in S/N, which was quantified more precisely using our model of the data in Denneulin (2020). To benefit from this improvement, the instrumental polarization has to remain small because incomplete polarimetric cycles do not benefit from the instrumental polarization compensation performed by estimating both couples Q and −Q (and U and −U, respectively).
6. Conclusion
In this work, we develop a new method for extracting the polarimetric signal using an inverse problems approach that exploits a model of the measured signal formation process. The method includes a weighted data fidelity term, which takes into account the blur and the polarization due to the instrument, and effectively disentangles polarized signal of interest from stellar contamination. In order to avoid noise amplification in the minimization of the data fidelity term, an edge preserving smoothing penalization was added allowing to favor smooth estimates almost everywhere. The associated minimization problem is solved by standard optimization techniques. Our method enables to accurately measure the polarized intensity and angle of linear polarization of circumstellar disks by taking into account the noise propagation and the observed objects convolution. It has the capability to use incomplete polarimetry cycles (when the instrumental polarization is small), which enhances the sensitivity of these observations. It also takes proper account of bad pixels by using weighted maps instead of interpolating them. These bad pixels can cause systematic errors of several tenths of a percent in the polarization measurements as shown by van Holstein et al. (2021) In addition, the effect of bad pixel interpolation could also have some impact when reaching 0.1% polarimetric accuracy.
We validated the method on both simulated and archive data from SPHERE/IRDIS and compared its performance with stateoftheart methods. We implemented the method in an endtoend dataanalysis package called RHAPSODIE. The method we developed improves the overall performances in particular at low S/Nandsmall polarized flux compared to standard methods.
By increasing the sensitivity and including deconvolution, this method will allow for more accurate studies of the orientation and morphology of the disks, especially in the innermost regions. It also will enable more accurate measurements of the angle of linear polarization at low S/N, which would allow for more indepth studies of dust properties. Finally, the method will enable more accurate measurements of the polarized intensity which is critical to construct the scattering phase functions.
RHAPSODIE is the first regularized inverse approach implemented for highcontrast polarimetric imaging. It demonstrates the benefits of advanced signal processing methods in this domain.
The code of RHAPSODIE is available online at https://github.com/LaurenceDenneulin/Rhapsodie.jl
Acknowledgments
We acknowledge Rob Van Holstein for his help with the instrumental polarization correction based on careful instrumental calibrations he performed which are implemented in the IRDAP tool he developed. We also thank the anonymous referee for her/his careful reading of the manuscript as well as her/his insightful comments and suggestions. This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon (OSUL/CRAL). This work is supported by the French National Programms PNP and the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU cofunded by CNES. The authors are grateful to the LABEX Lyon Institute of Origins (ANR10LABX0066) of the Université de Lyon for its financial support within the program Investissements d’Avenir (ANR11IDEX0007) of the French government operated by the National Research Agency (ANR). This paper is based on observations made with ESO Telescopes at the La Paranal Observatory under programme ID: 598.C0359, 095.C0273, 0102.C0916, 096.C0523, 097.C0523, 096.C0248.
References
 Adam, R., Ade, P. A., Aghanim, N., et al. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Akiyama, K., Alberdi, A., Alef, W., et al. 2019, ApJ, 875, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Avenhaus, H., Quanz, S. P., Schmid, H. M., et al. 2014, ApJ, 781, 87 [Google Scholar]
 Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berdeu, A., Soulez, F., Denis, L., Langlois, M., & Thiébaut, É. 2020, A&A, 635, A90 [CrossRef] [EDP Sciences] [Google Scholar]
 Beuzit, J.L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Biraud, Y. 1969, A&A, 1, 124 [NASA ADS] [Google Scholar]
 Birdi, J., Repetti, A., & Wiaux, Y. 2018, MNRAS, 478, 4442 [CrossRef] [Google Scholar]
 Birdi, J., Repetti, A., & Wiaux, Y. 2020, MNRAS, 492, 3509 [Google Scholar]
 Boccaletti, A., Sezestre, E., Lagrange, A. M., et al. 2018, A&A, 614, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borde, P. J., & Traub, W. A. 2006, ApJ, 638, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Carrillo, R. E., McEwen, J. D., & Wiaux, Y. 2012, MNRAS, 426, 1223 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Charbonnier, P., BlancFéraud, L., Aubert, G., & Barlaud, M. 1997, IEEE Trans. Image Process., 6, 298 [Google Scholar]
 Chierchia, G., Pustelnik, N., PesquetPopescu, B., & Pesquet, J.C. 2014, IEEE Trans. Image Process., 23, 5531 [Google Scholar]
 Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, in Groundbased and Airborne Instrumentation for Astronomy II, Int. Soc. Opt. Photonics, 7014, 70143E [Google Scholar]
 Combettes, P., & Pesquet, J. C. 2011, Fixedpoint algorithms for inverse problems in science and engineering, (Springer) 185 [Google Scholar]
 Combettes, P. L., & Wajs, V. R. 2005, Multiscale Model. Simul., 4, 1168 [CrossRef] [Google Scholar]
 de Boer, J., Salter, G., Benisty, M., et al. 2016, A&A, 595, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Boer, J., Langlois, M., van Holstein, R. G., et al. 2020, A&A, 633, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deledalle, C.A., Vaiter, S., Fadili, J., & Peyré, G. 2014, SIIMS, 7, 2448 [Google Scholar]
 Denneulin, L. 2020, PhD Thesis, Lyon, France [Google Scholar]
 Denneulin, L., Langlois, M., Pustelnik, N., & Thiébaut, E. 2019, in GRETSI LILLE 2019 [Google Scholar]
 Denneulin, L., Pustelnik, N., Langlois, M., Loris, I., & Thiébaut, E. 2020, iTwist, Nantes, France, Dec. 2–4, 2020 [Google Scholar]
 Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in Groundbased and Airborne Instrumentation for Astronomy II, Int. Soc. Opt. Photonics, 7014, 70143L [Google Scholar]
 Eldar, Y. C. 2008, IEEE Trans. Signal Process., 57, 471 [Google Scholar]
 Esposito, T. M., Kalas, P., Fitzgerald, M. P., et al. 2020, AJ, 160, 24 [Google Scholar]
 Flasseur, O., Denis, L., Thiébaut, É., Olivier, T., & Fournier, C. 2019, in 2019 27th European Signal Processing Conference (EUSIPCO) (IEEE), 1 [Google Scholar]
 Garufi,, A., Benisty, M., Stolker, T., et al. 2017, The Messenger, 169, 32 [Google Scholar]
 Ginski, C., Stolker, T., Pinilla, P., et al. 2016, A&A, 595, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Golub, G. H., Heath, M., & Wahba, G. 1979, Technometrics, 21, 215 [CrossRef] [Google Scholar]
 Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
 Hansen, P. C., & O’Leary, D. P. 1993, SISC, 14, 1487 [Google Scholar]
 Högbom, J. 1974, A&AS, 15, 417 [Google Scholar]
 Kalas, P. G., Rajan, A., Wang, J. J., et al. 2015, ApJ, 814, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lagrange, A. M., Langlois, M., Gratton, R., et al. 2016, A&A, 586, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Langlois, M., Dohlen, K., Vigan, A., et al. 2014, in Groundbased and Airborne Instrumentation for Astronomy V, Int. Soc. Opt. Photonics, 9147, 91471R [Google Scholar]
 Langlois, M., Pohl, A., Lagrange, A.M., et al. 2018, A& A, 614, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lefkimmiatis, S., Ward, J. P., & Unser, M. 2013, IEEE Trans. Image Process., 22, 1873 [Google Scholar]
 Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proc. Natl. Acad. Sci., 111, 12661 [NASA ADS] [CrossRef] [Google Scholar]
 Mahalanobis, P. C. 1936, National Institute of Science of India [Google Scholar]
 Maire, A. L., Langlois, M., Dohlen, K., et al. 2016, in Groundbased and Airborne Instrumentation for Astronomy VI, Int. Soc. Opt. Photonics, 9908, 990834 [Google Scholar]
 Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
 Milli, J., Engler, N., Schmid, H. M., et al. 2019, A&A, 626, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Molina, R. 1994, IEEE Trans. Pattern Anal. Mach. Intell., 16, 1122 [Google Scholar]
 Muto, T., Grady, C., Hashimoto, J., et al. 2012, ApJ, 748, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Nocedal, J., & Wright, S. J. 1999, Numerical Optimization, Springer Series in Operations Research (New York: Springer) [Google Scholar]
 Olofsson, J., Samland, M., Avenhaus, H., et al. 2016, A&A, 591, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pairet, B., Jacques, L., & Cantalloube, F. 2019, Proceedings of SPARS’19, 1, 1 [Google Scholar]
 Perrin, M. D., Duchene, G., MillarBlanchaer, M., et al. 2015, ApJ, 799, 182 [Google Scholar]
 Pinte, C., van der Plas, G., Menard, F., et al. 2019, Nat. Astron., 3, 1109 [Google Scholar]
 Pohl, A., Sissa, E., Langlois, M., et al. 2017, A&A, 605, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Price, D. J., Cuello, N., Pinte, C., et al. 2018, MNRAS, 477, 1270 [Google Scholar]
 Pustelnik, N., BenazzaBenhayia, A., Zheng, Y., & Pesquet, J.C. 2016, Wiley Encyclopedia of Electrical and Electronics Engineering (Hoboken, NJ, USA: John Wiley & Sons, Inc.), 1 [Google Scholar]
 Quanz, S. P., Avenhaus, H., Buenzli, E., et al. 2013, ApJ, 766, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Ramani, S., Blu, T., & Unser, M. 2008, IEEE Trans. Image Process., 17, 1540 [Google Scholar]
 Rudin, L. I., Osher, S., & Fatemi, E. 1992, Physica D: Nonlinear Phenomena, 60, 259 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Schmid, H. M., Bazzon, A., Roelfsema, R., et al. 2018, A&A, 619, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sissa, E., Olofsson, J., Vigan, A., et al. 2018, A&A, 613, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smirnov, O. M. 2011, A&A, 527, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Starck, J.L., Donoho, D. L., & Candès, E. J. 2003, A&A, 398, 785 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stein, C. M. 1981, Ann. Stat., 1135 [Google Scholar]
 Tarantola, A. 2005, Inverse Problem Theory and Methods for Model Parameter Estimation (SIAM) [Google Scholar]
 Thiébaut, E. 2002, in Astronomical Data Analysis II, Int. Soc. Opt. Photonics, 4847, 174 [Google Scholar]
 Thiébaut, E., & Conan, J.M. 1995, JOSA A, 12, 485 [Google Scholar]
 Tikhonov, A. N. 1963, Sov. Math. Dokl. [Google Scholar]
 Tinbergen, J. 2005, Astronomical Polarimetry (Cambridge University Press), googleBooksID: SAS4JzAaMxkC [Google Scholar]
 Titterington, D. M. 1985, A&A, 144, 381 [NASA ADS] [Google Scholar]
 van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132 [NASA ADS] [CrossRef] [Google Scholar]
 van Holstein, R. G., Snik, F., Girard, J. H., et al. 2017, Proc. SPIE, 10400, 1040015 [Google Scholar]
 van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020, A&A, 633, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Holstein, R. G., Stolker, T., JensenClem, R., et al. 2021, A&A, 647, A21 [EDP Sciences] [Google Scholar]
 Vigan, A., Langlois, M., Dohlen, K., et al. 2014, in Groundbased and Airborne Instrumentation for Astronomy V, Int. Soc. Opt. Photonics, 9147, 91474T [Google Scholar]
Appendix A: Synthetic dataset simulation
In order to evaluate the performance of the RHAPSODIE method, synthetic data have been created. These synthetic data are designed to reproduce astrophysical cases. First, the truth N = 128 × 128 maps , et θ_{gt} are created. Such a value of N pixels fits the main Region Of Interest (ROI) size.
The hardest disk structures to reconstruct are faint, small, and lightly polarized structures, and consequently have high contrast with the unpolarized stellar intensity and a low S/N. This is why the synthetic environment we generated a disk with three rings of equal brightness but with a different contrast with the unpolarized stellar intensity. This disk is partially polarized with a linearly polarization I^{p} and a polarization angle θ ∈ ] − π, π] and an unpolarized component I^{u disk}. The disk polarization ratio between both component is given by:
These synthetic images are then combined as maps of Stokes parameters I, Q, and U.
The I^{u disk} component is mixed with the unpolarized I^{u star} stellar components and a point source companion (star close to the host star of different brightness). The unpolarized intensity is represented as I^{u} = I^{u star} + I^{u disk}. It is important to keep in mind that unmixing both the disk and stellar unpolarized components is not possible from the DPI data without the diversity introduced by ADI. The τ^{disk} value used to synthetise datasets is thus inaccessible in practice from observational polarimetric datasets. To assert the case difficulty, one can then use the total polarization, given for all pixel n ∈ {1, …, N} by:
and the S/N, given for all pixel n ∈ {1, …, N} by:
where is the readout noise variance. The difference between τ^{total} and τ^{disk} is that the last one does not take into account the unpolarized star residuals. Figure A.1 presents the S/N maps and the maps of total polarization ratio of the synthetic parameters generated for different τ^{disk}. At the center, where the unpolarized star residual are the brightest, the S/N and the total polarization ratio are the weakest, especially in the case of small I^{u disk}. Yet the S/N grows with the separation from the star center (as with the stellar contribution or when the polarized contribution of the disk increases). Figure A.2 present the true simulated maps for τ^{disk} = 10%.
Fig. A.1. Total ratios of polarization τ^{total} maps, with respect to the total intensities, and S/N maps for the different values of disk polarized fractions τ^{disk}. 
Fig. A.2. Schematic illustration describing the process of the data simulation. Starting from synthetic maps , , and θ_{gt} based on a disc model, we generate artificial calibrated and preprocessed datasets d and associated weights W, illustrated here for τ^{disk} = 10%. 
Finally, to generate synthetic calibrated data, the Stokes maps are combined following the expression of the data physical model (12), with K × M noise realizations from a fixed random seed. The values of 10% of the pixels, chosen at random, are replaced by zeros to mimic bad pixel behaviours. The weights related to each acquisition are simulated at the same time following Eq. (15). The datasets are composed of 8 HWP cycles, with two acquisitions per positions in each cycles, leading to 16 images per position of the halfwave plate (HWP). Thus, the total number of images per datasets is K = 64.
Several datasets were created with τ^{disk} ∈ {3%, 7%, 10%, 15%, 25%}, corresponding to difficult cases for τ^{disk} < 10% and less difficult brighter cases above this threshold. Before producing a K × 2N noise realization for each dataset, the random seed is reset to the same value. This allow the reproducibility of the results. In fact this realization is obtained by the multiplication of the standard deviation of the pixel to a Gaussian, centered and reduced gaussian realization. Since the random seed is the same for each dataset, the realization is the same for the given pixel, only the standard deviation changes.
In order to compare the results of the double difference to the results of the RHAPSODIE methods the dataset are preprocessed. The bad pixels are interpolated; then the left and right part of the images are cut, recentered, and rotated.
All Tables
Positions of HWP α_{k} and orientations of the analyzer ψ_{j} and the corresponding values of the coefficients ν_{j, k, ℓ} assuming no instrumental polarization and ignoring field rotation.
Information related to the parameters used for the target reconstruction presented in this section.
All Figures
Fig. 1. Schematic view of the instrument ESO/VLTSPHERE IRDIS showing the various optical parts that can induce polarization effects. The same notations as in the text are used (e.g., n is the pixel index in the restored maps, k is the sequence index, and j is the polarizer index of the analyzer set). 

In the text 
Fig. 2. Comparison of the nonlinear RHAPSODIE reconstructions of RXJ 1615 (Avenhaus et al. 2018) for different values of λ_{Q + U} and μ_{Q + U}. The value of the GSURE criterion is indicated for each reconstruction. 

In the text 
Fig. 3. Example of the fitting of the true PSF on the empirical SPHERE/IRDIS data of the target HD 106906 reduced with RHAPSODIE. Observed PSF in H band (upper left) and different parametrization (Airy, Gaussian, and Moffat). Intensities are in logscale to enhance the faint diffraction patterns. 

In the text 
Fig. 4. Visual comparison of the reconstructed polarized intensity I^{p} with the stateoftheart double difference and the RHAPSODIE methods without deconvolution. 

In the text 
Fig. 5. Maps of errors of the reconstructions displayed in Fig. 4. These errors are obtained as the difference between the true and the reconstructed images. 

In the text 
Fig. 6. Comparison of the MSE between the true map and the estimated map of the polarized intensity I^{p} and of the angle of polarization θ for the double difference and the linear and nonlinear RHAPSODIE methods without deconvolution. 

In the text 
Fig. 7. Visual comparison of the reconstructed polarized intensity I^{p} with the stateoftheart double difference and the RHAPSODIE methods with deconvolution. 

In the text 
Fig. 8. Maps of errors of the reconstructions displayed in Fig. 7. These errors are obtained as the difference between the true and the reconstructed images. 

In the text 
Fig. 9. Comparison of the MSE between the true map and the estimated map of the polarized intensity I^{p} and of the angle of polarization θ for the double difference and the linear and nonlinear RHAPSODIE methods with deconvolution. 

In the text 
Fig. 10. Reconstruction of the Q (a) and U (b) parameters of the target TW Hydrae, for the first three cycle of HWP rotation, without (upper row) and with (lower row) the correction of the polarization. Without the correction, both Q and U are rotated and attenuated. The Q_{ϕ} and U_{ϕ} images reconstructed from the entire dataset are presented in (c). All the reconstructions are done without deconvolution to demonstrate mainly the efficiency of the instrumental polarization correction and the benefits of RHAPSODIE. (c) Azimutal Stokes parameters displayed in arcseconds: (i) Double Difference, (ii) Double Difference with the instrumental polarization correction, (iii) RHAPSODIE which includes by default the instrumental polarization correction. The intensities are multiplied in each pixel by the distance to the star r^{2}. 

In the text 
Fig. 11. Reconstructions of the polarized intensity I^{p} of the protoplanetary disks TW Hydrae, IM Lupus, and MY Lupus. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 

In the text 
Fig. 12. Reconstructions of the polarization angle θ of the protoplanetary disks TW Hydrae, IM Lupus, and MY Lupus. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 

In the text 
Fig. 13. Reconstructions of the polarized intensity I^{p} of the transition disks RY Lupus, T Chae, and RXJ 1615. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 

In the text 
Fig. 14. Reconstructions of the polarization angle θ of the transition disks RY Lupus, T Chae, and RXJ 1615. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 

In the text 
Fig. 15. Reconstructions of the polarized intensity I^{p} of the debris disks HD 106906, HD 61005, and AU MIC. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIE without deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 

In the text 
Fig. 16. Reconstructions of the polarized intensity I^{p} of the debris disks HD 106906, HD 61005, and AU MIC. From the left to the right, the reconstructions have been obtained with the Double Difference, RHAPSODIEwithout deconvolution and RHAPSODIE with deconvolution. The maps are displayed in logarithmic scale and normalized in contrast to the unpolarized stellar flux. The pixels lying underneath the coronograph are masked in black. North is up and east is to the left in all frames. 

In the text 
Fig. 17. Q_{ϕ} and U_{ϕ} projections of the reconstructions of RY Lupus and T Chae. The maps are displayed in symmetrical logarithmic scale and normalized in contrast to the unpolarized stellar flux. The coronograph is masked in white. 

In the text 
Fig. A.1. Total ratios of polarization τ^{total} maps, with respect to the total intensities, and S/N maps for the different values of disk polarized fractions τ^{disk}. 

In the text 
Fig. A.2. Schematic illustration describing the process of the data simulation. Starting from synthetic maps , , and θ_{gt} based on a disc model, we generate artificial calibrated and preprocessed datasets d and associated weights W, illustrated here for τ^{disk} = 10%. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.