RHAPSODIE: Reconstruction of High-contrAst Polarized SOurces and Deconvolution for cIrcumstellar Environments

Polarimetric imaging is one of the most effective techniques for high-contrast imaging and characterization of circumstellar environments. These environments can be characterized through direct-imaging polarimetry at near-infrared wavelengths. The SPHERE/IRDIS instrument installed on the Very Large Telescope in its dual-beam polarimetric imaging (DPI) mode, offers the capability to acquire polarimetric images at high contrast and high angular resolution. However dedicated image processing is needed to get rid of the contamination by the stellar light, of instrumental polarization effects, and of the blurring by the instrumental point spread function. We aim to reconstruct and deconvolve the near-infrared polarization signal from circumstellar environments. We use observations of these environments obtained with the high-contrast imaging infrared polarimeter SPHERE-IRDIS at the VLT. We developed a new method to extract the polarimetric signal using an inverse approach method that benefits from the added knowledge of the detected signal formation process. The method includes weighted data fidelity term, smooth penalization, and takes into account instrumental polarization. The method enables to accurately measure the polarized intensity and angle of linear polarization of circumstellar disks by taking into account the noise statistics and the convolution of the observed objects by the instrumental point spread function. It has the capability to use incomplete polarimetry cycles which enhance the sensitivity of the observations. The method improves the overall performances in particular for low SNR/small polarized flux compared to standard methods.


Introduction
With the adaptive-optics-fed high-contrast imaging instruments GPI (Macintosh et al. 2014) and SPHERE-IRDIS (Beuzit et al. 2019;Dohlen et al. 2008), we now have access to the spatial resolution and sensitivity required to observe in the near-infrared (NIR) circumstellar matter at small angular separations. Along with the Integral Field Spectrograph (IFS; Claudi et al. 2008) and the Zürich IMaging PO-Larimeter (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 so-called Dual Band Imaging (DBI) mode , or for two different polarizations, in a so-called Dual Polarimetry Imaging mode de Boer et al. 2020). Both circumstellar disks and self-luminous giant exoplanets or companion brown dwarfs can be characterized by these new instruments in direct-imaging 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 proven to be very successful for 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, in prep.) when they are surrounded by circumsubstellar disks.
Three particular types of circumstellar disks are studied: protoplanetary disks, transition disks, and debris disks. The observations of the protoplanetary and transition disks morphology linked to hydrodynamical simulations allows for the study of their formation scenario as in the study of HD 142527 (Price et al. 2018), IM Lup, RU Lup (Avenhaus et al. 2018), GSC 07396-759 (Sissa et al. 2018). Their observations are valuable because their shapes can be 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 such as for RXJ 1615, MY Lup, PDS 66 (Avenhaus et al. 2018), PDS 70 (Keppler et al. 2018(Keppler et al. , 2019Haffert 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 Article number, page 1 of 24 arXiv:2105.05546v1 [astro-ph.IM] 12 May 2021 A&A proofs: manuscript no. main disks are composed of dust and grain never accreted into 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 near-infrared and in the visible. Yet, 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 uses 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. Yet, such a method does not allow a good reconstruction of the disk morphology as it is impacted by artifacts due to self-subtraction. Moreover, the method fails when the environment is nearly rotation invariant. The DPI observations allow having 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 state-of-the-art methods to process datasets in polarimetry, apart from the calibration, are "step-by-step" 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 deal with the noise statistics. This results in some limitations in sensitivity in the case of low Signal-tonoise Ratio (SNR). Last, 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 previous processing, the results are not optimal given the available data. Still, these state-of-the-art methods have proven over the years to be sufficiently efficient to produce good 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;. Quantitative polarimetric measurements of circumstellar disks and substellar companions are currently very challenging, because existing data-reduction methods do not estimate properly 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 , calibrating the instrumental polarization effects with a sufficiently high accuracy already yield to 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 "step-by-step" 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 physical parameters estimation (Titterington 1985), mostly in adaptive optics (Borde & Traub 2006) and radio-interferometry with the well known 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 Polca-SARA (Birdi et al. 2019) in polarimetric radio-interferometry using more sophisticated tools as nonsmooth penalizations. A non-smooth method was also used for images denoising with curvelets (Starck et al. 2003). The minimization of a co-log-likelihood 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 high contrast imaging, the use of inverse problem methods is more recent. It has been used to perform auto-calibration 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. Yet, such reconstruction methods have not been used in polarimetric high contrast 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 Section 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 Section 3, we describe RHAPSODIE 1 (Reconstruction of High-contrAst Polarized Sources and Deconvolution for near-Infrared Environments), the method we developed. Section 4 is dedicated to the calibration of the detector, the instrument, and the instrumental polarization. Finally, in Section 5, we present the results obtained with RHAPSODIE on both simulated and astrophysical data.

Modeling polarimetric data
The principal parameters of interest for studying 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 param-  Fig. 1: Schematic view of the instrument ESO/VLT-SPHERE 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). eters 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 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. In the end, it is possible to reconstruct the parameters of interest I u , I p , and θ from a combination of the Stokes parameters.

Polarization effects
The four Stokes parameters S = (I, Q, U, V ) ∈ R 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 where S det j are the Stokes parameters on the detector after the left (j = 1) or right (j = 2) polarizer of the analyzer set while S are 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: M pol j for the left or right polarizer of the analyzer set, M der for the optical derotator, M HWP is for the half-wave plate (HWP), M M4 for the 4-th mirror of the telescope and M UT for the 3 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 alt-azimuthal telescope.
In the optical, only the total intensity I j det out of the Stokes parameters S j det 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 I det j , rotating the angle α of the HWP introduces a modulation of the contribution of the Stokes parameters Q and U in I det j , 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 quarter-wave 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 k-th 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 : where ψ j is the orientation angle of the left/right polarizer while α k is the HWP angle during the k-th acquisition. Table 1 lists the values of the linear coefficients ν j,k, for a typical set of HWP angles.

Parameters of interest
The model of the detected intensity given in Eq.
(3) is linear in 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: Article number, page 3 of 24 A&A proofs: manuscript no. main α k (HWP) ψ j (Analyzer) ν j,k,1 ν j,k,2 ν j,k,3 Table 1: This table lists the positions of HWP α k and the orientations of the analyzer ψ j and the corresponding values of the coefficients ν j,k, assuming no instrumental polarization and ignoring field rotation. It is computed from Eq. (4). 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 non-linear function of the parameters I u , I p and θ: I det j,k = ν j,k,1 I u + ν j,k,1 + ν j,k,2 cos(2 θ) + ν j,k,3 sin(2 θ) I p .
For a perfect SPHERE/IRDIS-like instrument and not considering field rotation, combining Eqs. (4) and (5) yields: which is the Malus law.

Accounting for the instrumental spatial PSF
In polarimetric imaging, each polarimetric parameter is a function of the 2-dimensional FOV. We consider that the Stokes parameters are represented by images of N pixels each and denote as S ,n the value of the n-th 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 write the detected intensity for a given detector pixel as: where I det j,k,m is the intensity measured during the k-th acquisition by the m-th pixel of the sub-image corresponding Polarizer of the analyser set k ∈ 1, K Data frame in the sequence ∈ 1, L Component the parameter of interest. m ∈ 1, M Pixel in data sub-image n, n ∈ 1, N Pixel in restored model maps ρ Param. to regularize (I, I u , I p , Q + U ) to the j-th 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. It follows from our assumption on spatial and polarization effects applying 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 shift-invariant and can thus be modeled by convolution with a shift-invariant PSF. Geometrical transforms and convolutions do not commute but their order may be changed provided the shift-invariant PSFs are appropriately rotated and/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 shift-invariant blurs followed by a single geometrical transform accounting for all the rotations but also possible geometrical translations and/or spatial (de)magnification 3 . Following this analysis, our model of the spatial PSF is given by: where A k : R N → R N implements the shift-invariant blur of the input model maps while T j,k : R N → R Mj performs the geometrical transform of the blurred model maps for the j-th polarizer of the analyzer set during the k-th acquisition. N is the number of pixels in the model maps and M j is the number of pixels of the detector (or sub-image) corresponding to the j-th output polarizer. In our implementation of the PSF model, the geometrical transform of images is performed using interpolation by Catmull-Rom 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 diag(p k ) implements the frequency-wise multiplication byp k = F p k , the discrete Fourier transform of p k , the shift-invariant PSF. Note that p k must be specified in the same reference frame as the FOV. If the shift-invariant PSF is calibrated from empirical images of, e.g., the host star, acquired by the detector, the inverse (or pseudo-inverse) of T j,k must be applied to the empirical images.

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 • }. Besides, 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, sub-images, 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 pre-processing 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 I det j,k,m is given in Eq. (9) for k ∈ 1, K the index of the acquisition, j ∈ {1, 2} indicating the left/right polarizer of the analyzer, and m ∈ 1, M the pixel index in the corresponding left/right sub-image.
The ≈ sign in Eq. (12) is to account for an unknown random perturbation term due to the noise. Noise in the pre-processed images can be assumed centered and independent between two different pixels or frames because the pre-processing 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 read-out noise, 3 If the pixel size of the polarimetric maps is not chosen to be equal to the angular size of the detector pixels etc. For most 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 non-uniform Gaussian distribution whose mean is given by the right-hand-side term of Eq. (12) and whose variance Σ j,k,m = Var(d j,k,m ) is estimated in the calibration stage (see subsection 4.1).
It is worth noticing 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, non-linear 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 Section 3.2.

RHAPSODIE : Reconstruction of High-contrAst
Polarized SOurces and Deconvolution for cIrcumstellar Environments

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. To remain as general as possible, we denote by X ∈ R N ×L the set of parameters of interest to be recovered and by X ∈ R N , with ∈ 1, L , the -th parameter component which is a N -pixel map. Given the direct model of the pre-processed 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 X as the ones that minimize a given objective function f (X) possibly under constraints expressed as X ∈ C with C the set of acceptable solutions. The objective function takes the form of the sum of a data-fidelity term f data (X) and of regularization terms f ρ (X): where λ ρ ≥ 0 (∀ρ) are so-called hyperparameters introduced to tune the relative importance of the regularization terms. The data-fidelity 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 ill-conditioning of the inverse problem. Additional strict constraints may be imposed on the sought parameters via the feasible set C, e.g. to account for the requirement that intensities are nonnegative quantities. These different terms and constraints are detailed in the following sub-sections.

Data fidelity
Knowing the sufficient statistics for the pre-processed data, agreement of the model with the data is properly insured by the co-log-likelihood of the data (Tarantola 2005) or equivalently by the following criterion: where · 2 W = ·, W· denotes Mahalanobis (1936) squared norm, d j,k = (d j,k,1 , . . . , d j,k,M ) ∈ R M collects all the pixels (e.g., in lexicographic order) of the j-th sub-image in the k-th acquisition as defined in Eq. (12). Similarly, I det j,k (X) = (I det j,k,1 , . . . , I det j,k,M ) ∈ R M where the terms I det j,k,m 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 X = (I u , I p , θ), 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 non-uniform noise variance and for invalid data (see 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 because of too strong atmospheric turbulence or unproper coronograph centering. 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 sub-images are mutually independent.

Regularization
The problem of recovering the polarimetric parameters from the data is an ill-conditioned inverse problem mainly due to the instrumental blur. Furthermore, the problem may also be ill-posed if there are too many invalid data. In the case of an ill-conditioned 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 data-fidelity as assumed by the objective function defined in Eq. (13).

Edge-preserving smoothness
We expect that the light distribution of circumstellar environments be mostly smooth with some sharp edges, hence edge-preserving 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 I u ∈ R N 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 : R N → R 2 is a linear operator which yields an approximation of the 2D spatial gradient of its argument around the n-th pixel. This operator is implemented by means of finite differences; more specifically applying D n to a sampled map u of a given parameter writes: where (n 1 , n 2 ) denotes the row and column indices of the n-th 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 pseudo-norm 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 sub-L 2 norm of the gradient like the last one in Eq. (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 so-called cartoon effect (i.e., piecewise flat images) which is not appropriate for realistic astronomical images. We however note that the hyperbolic edge-preserving 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 non-quadratic but differentiable objective functions of many (millions or even billions) variables possibly with additional strict constraints (Thiébaut 2002). See Denneulin et al. (2019Denneulin et al. ( , 2020; Denneulin (2020) for a comparison of possible advanced regularizers.

Tuning of the hyperparameters
In the regularization function f prior (X), the terms defined in Eq. (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 Cross-Validation (GCV: Golub et al. 1979), the L-curve (Hansen & O'Leary 1993), or 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. Fig. 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 λ ρ /µ ρ : 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 high contrast contexts, we also observed this tendency of GSURE to over-smooth the result. We think that devising a good unsupervised method for tuning the hyperparameters deserve further studies, and, 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 po- µ Q+U = 10 −3.5 µ Q+U = 10 −2.5 µ Q+U = 10 −1.5 regularization of I p alone by f I p 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.

Imposing the positivity of the intensities
Imposing non-negativity 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 non-negative. 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 X = (I u , I p , θ), the positivity constraint writes: Expressed for the Stokes parameters X = (I, Q, U ), the positivity yields an epigraphical constraint: Article number, page 7 of 24 A&A proofs: manuscript no. main Such a constraint can be found in Birdi et al. (2018), but it has not yet been implemented in high contrast polarimetric imaging. Since I p n = Q 2 n + U 2 n (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:

Choice of the polarimetric parameters
Our method expresses the recovered parameters X as the solution of a constrained optimization problem specified in Eq. (13). As explained in Sec. 3.3.2, the weights of imposed regularization is chosen via the values of the multipliers λ ρ . The constraints can be implemented by the feasible C 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 Eq. (5) and Eq. (6) can be used to estimate any parameter of interest given the recovered X. 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 forward-backward 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 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 VMLM-B (Thiébaut 2002) can be used to solve the problem in Eq. (13). The VMLM-B algorithm is a quasi-Newton method with limited memory requirements and able to account for separable bound constraints. VMLM-B is applicable to large size problems and only require to provide 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 non-linear 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.

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 calibration data are expressed as a function of the different contributions (i.e., flux, sky background, instrumental background, gain, noise, and quantum efficiency). All these quantities are then jointly estimated by the minimization of the co-log-likelihood of the data. Calibrated data are corrected for contribution of the bias and the background and for non-uniform 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 .

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, simultaneous PSF estimation is not possible. To estimate the PSF, we use a dedicated flux calibration (STAR-FLUX) 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 to remain rotation invariant (see Fig. 3).
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 on Fig. 3). For astrophysical observations, we fit a 2D PSF model on pre-reduced PSF data for each observed target (e.g. Gaussian, Airy, and Moffat fits on 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 (STAR-CENTER) where four faint replicas of the star image are created by giving a bi-dimensional sinusoidal profile to the deformable mirror (see Beuzit et al. 2019). The STAR-CENTER calibration was repeated before and after each science observation, and resulting center estimations were Intensities are in log-scale to enhance the faint diffraction patterns.
averaged. In addition we use 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.

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 cross-talk contamination between 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 face-on disks, we instead use the method developed by van Holstein et al. (2020) which rely on the pre-computed 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, 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 model-corrected Q and U images.
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 model-corrected 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 cut-off frequency). Both corrections factors ε Q and ε U as follow: 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 model-corrected Q and U images.

Applications on high contrast polarimetric data
In this section, we will compare two RHAPSODIE configurations. The first configuration defined as without deconvolution, means that I det j,k,m in the data fidelity term (14) does not include the convolution, leading to the simplification of the equations (9) and (10): (i.e., A k is the identity). Such a configuration of RHAP-SODIE aims to be comparable to the state-of-the-art methods (i.e., Double Difference and IRDAP). These methods do not include any deconvolution and we aim to show that RHAPSODIE also performs well in such a case. The second configuration defined as with deconvolution states for the full RHAPSODIE capabilities. Such a configuration of RHAPSODIE aims to show the benefits of the global model compared to an a posteriori deconvolution to improve the angular resolution.

Application on synthetic data
The performance of the linear and non-linear methods are first evaluated on synthetic datasets, without (resp. with) the deconvolution displayed on Fig. 4, Fig. 5 and Fig. 6 (resp. Fig. 7, Fig. 8 and Fig. 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 the equation (A.1).
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 high-pass spatial filtering. The deconvolved Double Difference reconstructions are obtained by solving the following problem: where Q D. D. and U D. D. 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 close-by polarization signals with 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 I u , I p and θ, obtained from X from Eq. (6), and the ground truth I u gt , I p gt and θ gt . The total MSE is given by: where, N I u , N I p , and N θ are the number of pixels with signal of interest in the respective I u , I p , and θ maps. For the deconvolution of the Double Difference results, λ (Q+U ) D. D. and µ (Q+U ) D. D. are chosen to minimize only the sum of the MSE on I u and θ.
For the reconstructions without the deconvolution, Fig. 4 and Fig. 5 show that the reconstructions are less noisy with RHAPSODIE. The inner circle is always better reconstructed. Moreover, the non-linear reconstruction (i.e., minimization on I u , Q and U ) is better than the linear when τ disk grows. In both configurations of RHAP-SODIE (non-linear 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 trade-off between a smooth solution and a solution close to the data. Classically, minimizing the MSE is a good trade-off between underfitting and overfitting. The MSE curves in Fig. 6 are coherent with the observations. When τ disk = 25%, the non-linear 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 three-time smaller with RHAPSODIE compared to the Double Difference error.
For the reconstructions with the deconvolution, the RHAPSODIE methods better 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 RHAP-SODIE methods also allow us to achieve a better reconstruction of the angle of polarization, mostly with the nonlinear estimation. The non-linear reconstruction appears to be more efficient as being not polluted by the artifacts of deconvolution of the unpolarized point source companion as seen on Fig. 8. Even if according to the MSE in Fig. 9, for τ disk = 25%, the MSE is smaller for the linear RAP-SODIE. in fact, the space between the outer ring and the thin ring is better reconstructed with such a configuration (see Fig. 7 and Fig. 8).
According to these results, the RHAPSODIE methods are better than the state-of-the-art methods, in particular the non-linear 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 non-linear RHAPSODIE method with manual selection of the hyperparameters.      Table 4: Information of the datasets used for the target reconstruction: the name of the target, the date of the observation, the filter used, the exposition time ∆ t for one acquisition, the number of cycles of HWP, the number of frame missing or removed K miss , the total exposition time of the observation ∆ tot t considering K miss , the seeing, and observation conditions τ 0 . The total number of frames in each dataset is given by  , λ no-dec Q+U , µ no-dec (resp. λ dec I u , λ dec Q+U , µ dec ) used for the reconstruction with RHAPSODIE without deconvolution (resp. with deconvolution), and the parametrization model of the PSF used. map of the projected intensities, by using the standard azimuthal Stokes parameters Q φ and U φ estimated from: with ψ n = arctan n1−n center 1 n2−n center 2 where (n 1 , n 2 ) denotes the row and the column indices of the n-th pixel in the map and (n center 1 , n center 2 ) 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 the 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 drop-off of stellar illumination with distance.
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 face-on. 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 IR-DAP 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 cross-talk 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 SNR 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 calibration of observations more recent than 2016 have improved and do allow better calibration and as a consequence do not impact our method efficiency anymore.
The efficiency of the method we have developed is demonstrated on the Figures 11,12,13,14,15,16,and Fig. 17. Fig. 11 and Fig. 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). Fig. 13, 14, and Fig. 17 present RHAPSODIE reconstructions of transition disks RY Lupus (Langlois et al. 2018) and T Chamaeleontis (Pohl et al. 2017). Fig. 13, 14 also present RHAPSODIE reconstructions of transition disk RXJ 1615 . Fig. 15 and Fig. 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), andAU Mic (Boccaletti et al. 2018).These reconstructions are all normalized in contrast to the unpolarized stellar flux (estimated when registering the PSF off-centered from the coronograph) by taking into account the transmission of the neutral-density filter used when registering the PSF to prevent saturation.
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 on 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 high-quality reconstructions of the disk polarized signal and minimizes the artifacts from bad pixels. The comparison of these reductions with the state-of-the-art confirms the benefits of our method but it is more difficult to quantify the gain on real data not knowing the truth and should be strengthened by using numerical models of 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 close-by 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 before the hyperparameters could be further tuned to adjust the smoothing of the deconvolved  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 (on Fig. 13) which was not identified by Pohl et al. (2017). The reason 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 disk presented on known East-West 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 backward-scattering 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. The reader may remark that the RHAPSODIE deblurred reconstructions of RXJ 1615 (on Fig. 13) and HD 106906 (on Fig. 15) seems noisier than without deconvolution. This is due to hyperparameters 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) taken under very bad conditions. For these two datasets, our method also proves to be efficient in using incomplete polarization cycles to recover the disk polarized signal, and to deconvolve this signal despite strong artifacts produced in addition 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 which is created by these spiders remains, as seen on 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 SNR 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).

Conclusion
We developed a new method to extract 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 am-plification in the minimization of the data fidelity term, an edge preserving smoothing penalization is 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 for bad pixels by using weighted maps instead of interpolating them. These bad pixels can cause systematic errors of several tenths of 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 have validated the method on both simulated and archive data from SPHERE/IRDIS and compared its performances with the state-of-the-art methods. We have implemented the method in an end-to-end data-analysis package called RHAPSODIE. The method we developed improves the overall performances in particular at low SNR/small polarized flux compared to standard methods.
By increasing the sensitivity and including deconvolution, this method will allow for more accurate studies of the A&A proofs: manuscript no. main τ total τ disk = 3% τ disk = 7% τ disk = 10% τ disk = 15% τ disk = 25% τ disk = 50% Starting from synthetic maps I u gt , I p gt , and θ gt based on a disc model, we generate artificial calibrated and pre-processed datasets d and associated weights W, illustrated here for τ disk = 10%.