L Denneulin

M Langlois

É Thiébaut

N Pustelnik

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

Context. 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 nearinfrared wavelengths. The Spectro-Polarimetric High-contrast Exoplanet REsearch (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. Aims. We aim to reconstruct and deconvolve the near-infrared polarization signal from circumstellar environments. Methods. We use observations of these environments obtained with the high-contrast imaging infrared polarimeter SPHERE-IRDIS at the Very Large Telescope (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 Results. 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 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. 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 SNR, which would lead to in-depth studies of dust properties. Finally, the method will enable more accurate measurements of the polarized intensity which is critical to construct scattering phase functions.

Introduction

With the adaptive-optics-fed high-contrast imaging instruments GPI (Macintosh et al. 2014) and SPHERE-IRDIS (Beuzit et al. 2019;[START_REF] Dohlen | Ground-based and Airborne Instrumentation for Astronomy II[END_REF], 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; [START_REF] Claudi | Groundbased and Airborne Instrumentation for Astronomy II[END_REF]) 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 [START_REF] Vigan | Ground-based and Airborne Instrumentation for Astronomy V[END_REF], or for two different polarizations, in a so-called Dual Polarimetry Imaging mode [START_REF] Langlois | Ground-based and Airborne Instrumentation for Astronomy V[END_REF]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 [START_REF] Garufi | [END_REF] and shows much promise for the characterization of brown dwarfs (van Holstein et al. 2017) and exoplanets (van Holstein 2020, 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. , 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 [START_REF] Langlois | [END_REF]) 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 never accreted into the planets, as in the case of HR 4796A [START_REF] Perrin | [END_REF]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; [START_REF] Marois | [END_REF] 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 [START_REF] Tinbergen | Astronomical Polarimetry[END_REF]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;[START_REF] Ginski | [END_REF]de Boer 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 (van Holstein et al. 2020), 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 [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Combettes | [END_REF]Pustelnik et al. 2016). Such methods have been used over decades in astrophysics for the physical parameters estimation [START_REF] Titterington | [END_REF], 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 [START_REF] Thiébaut | [END_REF]. 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 using ADI data [START_REF] Pairet | Proceedings of SPARS'19[END_REF][START_REF] Flasseur | 2019 27th European Signal Processing Conference (EUSIPCO)[END_REF] 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 RHAP-SODIE (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 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 Optical devices HWP Optical devices Analysers

ψ 1 = 0 • ψ 2 = 90 • j = 1 j = 2 I u n I p n θn { I { Q, U × × d j,k,m { ν j,k, 0 • • 45 • • 22, 5 • • 77, 5 • • α k
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).

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) can be modeled by (see Eq. ( 17) of van Holstein et al. 2020):

S det j = M pol j T -θ der M der T θ der T(-α) M HWP T(α) M M4 T θ alt M UT T(θ par ) S (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

α k (HWP) ψ j (Analyzer) ν j,k,1 ν j,k,2 ν j,k,3 0 • 0 • (left : j = 1) 1 / 2 1 / 2 0 0 • 90 • (right : j = 2) 1 / 2 -1 / 2 0 45 • 0 • (left : j = 1) 1 / 2 -1 / 2 0 45 • 90 • (right : j = 2) 1 / 2 1 / 2 0 22.5 • 0 • (left : j = 1) 1 / 2 0 1 / 2 22.5 • 90 • (right : j = 2) 1 / 2 0 -1 / 2 77.5 • 0 • (left : j = 1) 1 / 2 0 -1 / 2 77.5 • 90 • (right : j = 2) 1 / 2 0 1 / 2
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).

detector. It then follows from Eq. ( 1) that the quantity measured by the detector is a simple linear combination of the input Stokes parameters:

I det j = 4 =1 ν j, (Θ) S (2)
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:

I det j,k = 3 =1 ν j,k, S (3)
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 (van Holstein et al. 2020):

I det j,k = 1 2 I + 1 2 cos 4 α k +2 ψ j Q+ 1 2 sin 4 α k +2 ψ j U , (4
) 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:

   I = I u + I p Q = I p cos(2 θ) U = I p sin(2 θ) (5)
and conversely by 1 :

     I p = Q 2 + U 2 θ = (1/2) arctan (U/Q) mod π I u = I -Q 2 + U 2 . (6)
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 . ( 7
)
For a perfect SPHERE/IRDIS-like instrument and not considering field rotation, combining Eqs. ( 4) and (5) yields:

I det j,k = 1 2 I u + I p cos 2 θ -2 α k -ψ j , (8)
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:

I det j,k,m = 3 =1 N n=1 ν j,k, H j,k,m,n S ,n , (9)
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 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}, 1 In this representation there is a ±180 • degeneracy for the linear polarization angle θ.

Notation Description d j,k,m

Measured data, Eq.( 12)

I det j,k,m
Model of the data, Eq. ( 9)

H j,k,m,n
Instrumental PSF, Eq. (10) Σ j,k,m

Variance of data, Eq. (15) ν j,k, Polarization effects, Eq. (3) S = (I, Q, U ) Stokes parameters, Eq. ( 5)

I u
Unpolarized intensity, Eq. ( 6)

I p
Polarized intensity, Eq. ( 6) θ angle of polarization, Eq. ( 6) X Parameters of interest, Eq. (13) X Estimated parameters, Eq. ( 13) X gt Ground truth, Eq. ( 28) f data Data fidelity term, Eq. ( 14)

W j,k,m
Weights, Eq. ( 15)

f ρ
Regularization term, Eq. ( 16)-( 20)

λ ρ
Regul. contribution, Eq. ( 16)-( 20)

µ ρ
Regularization threshold, Eq. ( 16)-( 20 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 2 . Following this analysis, our model of the spatial PSF is given by:

H j,k,m,n = N n =1 T j,k m,n A k n ,n , (10)
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:

A k = F -1 diag(p k ) F (11)
where F denotes a Fast Fourier Transform (FFT) operator of suitable size and diag(p k ) implements the frequencywise multiplication by pk = 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:

d j,k,m ≈ I det j,k,m (12)
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, 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 [START_REF] Titterington | [END_REF][START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF] 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):

X = arg min X∈C f (X) = f data (X) + ρ λ ρ f ρ (X) f prior (X) (13)
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 [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF] or equiv-alently by the following criterion:

f data (X) = j,k d j,k -I det j,k (X) 2 W j,k , (14)
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:

∀m ∈ 1, M , W j,k m,m = Σ -1 j,k,m
for valid data;

0 for invalid data, (15)
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 [START_REF] Titterington | [END_REF][START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF]. 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:

f I u (X) = n D n I u (X) 2 + µ 2 I u , (16)
f I p (X) = n D n I p (X) 2 + µ 2 I p , (17)
f I (X) = n D n I(X) 2 + µ 2 I , (18)
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:

D n u = u n1+1,n2 -u n1,n2 u n1,n2+1 -u n1,n2 . ( 19
)
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:

f Q+U (X) = n D n Q(X) 2 + D n U (X) 2 + µ 2 Q+U . (20)
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 [START_REF] Thiébaut | Astronomical Data Analysis II[END_REF]. See [START_REF] Denneulin | GRETSI LILLE 2019[END_REF]Denneulin 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 λ Unsupervised tuning of the hyperparameters with GSURE [START_REF] Eldar | [END_REF], 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 λ ρ /µ ρ :

λ ρ D n X ρ 2 + µ 2 ρ ≈ λ ρ µ ρ + λ ρ 2 µ ρ D n X ρ 2 . ( 21
)
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 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 f I u 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 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. µ Q+U = 10 -3.5 µ Q+U = 10 -2.5 µ Q+U = 10 -1.5

Hyperparameters λ I p λ I u λ I λ Q+U Set 1 = 0 ≥ 0 = 0 ≥ 0 Set 2 = 0 = 0 ≥ 0 ≥ 0 Table 3: Set of
(λ/µ) Q+U = 2 (λ/µ) Q+U = 3 (λ/µ) Q+U = 4

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:

C = (I u , I p , θ) ∈ R N ×3 ∀n ∈ 1, N , I u n ≥ 0, I p n ≥ 0 . ( 22
)
Expressed for the Stokes parameters X = (I, Q, U ), the positivity yields an epigraphical constraint:

C = (I, Q, U ) ∈ R N ×3 ∀n ∈ 1, N , I n ≥ Q 2 n + U 2 n . ( 23
)
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:

C = (I u , Q, U ) ∈ R N ×3 ∀n ∈ 1, N , I u n ≥ 0 . (24)

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 [START_REF] Combettes | [END_REF], 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 [START_REF] Thiébaut | Astronomical Data Analysis II[END_REF]) 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.

Data calibrations

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, thedefective 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 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. Intensities are in log-scale to enhance the faint diffraction patterns.

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.

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 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:

ε Q = n∈Ω Q n /I u n /N Ω ε U = n∈Ω U n /I u n /N Ω , (25)
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 det j,k,m = 3 =1 ν j,k, N n=1 (T j,k ) m,n S ,n , (26)
(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 [START_REF] Tinbergen | Astronomical Polarimetry[END_REF]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:

( Q, U ) ∈ arg min (Q,U )∈R N ×R N Q D. D. -AQ 2 + U D. D. -AU 2 + λ (Q+U ) D. D. f (Q+U ) D. D. (Q, U ) (27)
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:

MSE tot = N I u n=1 E I u -I u gt 2 + N I p n=1 E I p -I p gt 2 + N θ n=1 E angle e 2i( θ-θgt) 2 . ( 28
)
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. 7).

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.

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 stateof-the-art method. For each reconstruction, we present the map of the projected intensities, by using the standard azimuthal Stokes parameters Q φ and U φ estimated from: 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 K = ∆ tot t / ∆ t + K miss .

Q φ = Q n cos (2ψ n ) + U n sin (2ψ n ) U φ = U n cos (2ψ n ) -Q n sin (2ψ n ) , (29)
Target λ no-dec

I u λ no-dec Q+U µ no-dec I u µ no-dec Q+U λ dec I u λ dec Q+U µ dec I u
µ dec Q+U PSF TW Hydrae 10 4.0 10 4.0 10 -4.0 10 -4.0 10 4.0 10 4.0 10 -5.0 10 -5.0 Moffat IM Lupus 10 3.0 10 4.0 10 -0.8 10 -4.0 10 5.0 10 5.5 10 -3.0 10 -3.0 Airy MY Lupus 10 2.0 10 0.4 10 -0.7 10 -3.6 10 1.0 10 0.5 10 -3.0 10 -3.0 Moffat RY Lupus 10 4.0 10 5.5 10 -3.0 10 -3.0 10 5.0 10 3.0 10 -3.0 10 -3.0 Moffat T Chae 10 0.0 10 -1.0 10 -3.0 10 -3.0 10 0.5 10 -1.0 10 -2.0 10 -2.0 Moffat RXJ 1615 10 0.0 10 0.0 10 -4.0 10 1.0 10 1.0 10 -1.5 10 -3.0 10 -3.0 Moffat HD 106906 10 2.5 10 1.5 10 -4.0 10 -4.0 10 5.0 10 1.7 10 2.2 10 -1.8 Airy HD 61005 10 5.0 10 4.1 10 -0.7 10 -3.4 10 5.0 10 4.1 10 -0.7 10 -3.9 Moffat AU MIC 10 5.0 10 4.3 10 -4.2 10 -3.2 10 5.0 10 4.3 10 -1.2 10 -4.2 Moffat Table 5: Informations of the parameters used for the target reconstruction presented in this section: the name of the target and the values of the hyperparameters λ no-dec

I u
, λ 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.

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,and 16. Fig. 11 and Fig. 12 present the Double Difference and RHAPSODIE reconstructions of the protoplanery disks TW Hydrae [START_REF] Van Boekel | [END_REF], IM Lupus (Avenhaus et al. 2018) and MY Lupus (Avenhaus et al. 2018). Fig. 13 and 14 presents RHAPSODIE reconstructions of transition disks RY Lupus [START_REF] Langlois | [END_REF], T Chamaeleontis (Pohl et al. 2017), RXJ 1615(de Boer et al. 2016). Fig. 15 and Fig. 16 presents 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 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 images according to the required noise/angular resolution trade-off. For several of these disks (IM Lup, MY Lup on 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. 2020). The image clearly shows the 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 10906 (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 (HD61005 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 Hol-stein 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) 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 a percent in the polarization measurements as shown by (van Holstein 2020, in prep.) In addition, the effect of bad pixel interpolation could also have some impact when reaching 0.1% polarimetric accuracy.

demonstrates the benefits of advanced signal processing methods in this domain.

Acknowledgements

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

  ρ 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).

  hyperparameters used in the present work.

Fig. 2 :

 2 Fig. 2: Comparison of the non-linear 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.

Fig. 3 :

 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 log-scale to enhance the faint diffraction patterns.

Fig. 4 :Fig. 5 :Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :

 456789 Fig. 4: Visual comparison of the reconstructed polarized intensity I p with the state-of-the-art Double Difference and the RHAPSODIE methods without deconvolution.

  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.

  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 :

 13 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 :

 14 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 106 906, 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

  Fig. A.2: Schematic describing the process of the data simulation. 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%.

Table 2 :

 2 Notations.k ∈ 1, K and m ∈ 1, M . Table2summarizes the main notations used in this paper.

	)

  SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTH-EAS/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 co-funded by CNES. The authors are grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its financial support within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) 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.C-0359, 095.C-0273, 0102.C-0916, 096.C-0523, 097.C-0523, 096.C-0248.

Article number, page 2 of 23 Denneulin et al.: Reconstruction of high-contrast polarized circumstellar environments

if the pixel size of the polarimetric maps is not chosen to be equal to the angular size of the detector pixels

Article number, page 6 of 23 Denneulin et al.: Reconstruction of high-contrast polarized circumstellar environments

Article number, page 8 of 23 Denneulin et al.: Reconstruction of high-contrast polarized circumstellar environments

Article number, page 10 of 23 Denneulin et al.: Reconstruction of high-contrast polarized circumstellar environments

Article number, page 14 of 23 Denneulin et al.: Reconstruction of high-contrast polarized circumstellar environments

Article number, page 18 of 23 Denneulin et al.: Reconstruction of high-contrast polarized circumstellar environments

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 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 SNR, which would allow for more in-depth studies of dust properties. Finally, the method will enable more accurate measurements of the polarized intensity which is critical to construct scattering phase functions.

RHAPSODIE is the first regularized inverse approach implemented for high contrast polarimetric imaging. It 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 I u gt , I p gt 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 SNR. 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 udisk . 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 udisk component is mixed with the unpolarized I ustar stellar components and a point source companion (star close to the host star of different brightness). The unpolaized intensity is represented as I u = I ustar + I udisk . It is important to keep in mind that unmixing both disk and stellar unpolarized components is not possible from 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 Signal-to-Noise Ratio (SNR), given for all pixel n ∈ {1, . . . , N } by:

where σ 2 ro is the read-out noise variance. The difference between τ total and τ disk is that the last one does not take in account the unpolarized star residuals. The figure A.1 present the SNR 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 SNR and the total polarization ratio are the weakest especially in the case of small I udisk . Yet the SNR grows with the separation from the star center (like the stellar contribution or when the polarized contribution of the disk increases. 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. Are also introduced 10% of bad pixels choosen randomly. The weights related to each acquisition are simulated at the same time following (15). The datasets are composed of K αrot = 8 HWP cycles, with K αacc = 2 acquisitions per positions in each cycles, giving K α = 16 total images per position of Half-Wave Plate (HWP), for a total of K = 4K α = 64 images per dataset.

Several datasets are 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 K × 2N noise realization on each dataset, the random seed is reset to the same value. This allow the reproductibility 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 pre-processed. The bad pixels are interpolated; then the left and right part of the images are cut, recentered and rotated.