Open Access
Issue
A&A
Volume 653, September 2021
Article Number A138
Number of page(s) 22
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/202039618
Published online 24 September 2021

© L. Denneulin et al. 2021

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

With the 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 POLarimeter (ZIMPOL; Schmid et al. 2018) which can also be used to observe circumstellar environments in polarimetry, the IRDIS instrument is one of the three SPHERE instruments (Beuzit et al. 2019). SPHERE/IRDIS is able to acquire two simultaneous images at two different wavelengths, in a so-called Dual Band Imaging (DBI) mode (Vigan et al. 2014), or for two different polarizations, in a so-called Dual Polarimetry Imaging mode (Langlois et al. 2014; de Boer et al. 2020). Both circumstellar disks and self-luminous giant exoplanets or companion brown dwarfs can be characterized thanks to these new instruments in direct-imaging processes at these wavelengths.

The NIR polarimetric mode of SPHERE/IRDIS at the Very Large Telescope (VLT), which is described in Beuzit et al. (2019), de Boer et al. (2020), van Holstein et al. (2020), has been proven to be very successful in the detection of circumstellar disks in scattered light (Garufi, et al. 2017) and shows much promise for the characterization of brown dwarfs (van Holstein et al. 2017) and exoplanets (van Holstein et al. 2021, and in prep.) when they are surrounded by circumsubstellar disks.

There are three particular types of circumstellar disks that are the subject of recent studies: protoplanetary disks, transition disks, and debris disks. Observations of the protoplanetary and transition disks, in terms of their morphology and linked to hydrodynamical simulations, allow for the study of their formation scenario, as in studies of HD 142527 (Price et al. 2018), IM Lup, RU Lup (Avenhaus et al. 2018), and GSC 07396-759 (Sissa et al. 2018). Their observations are valuable because their shapes can serve as the signposts for the formation of one or several exoplanets. In fact, during their formation, the planets “clean” the dust off their orbits, creating gaps without dust, as in the case of RXJ 1615, MY Lup, PDS 66 (Avenhaus et al. 2018), and PDS 70 (Keppler et al. 2018, 2019; Haffert et al. 2019). The planet formation scenario can be explained with hydrodynamical simulations as the cases of HL Tau (Dipierro et al. 2015), HD 163296 (Pinte et al. 2019). Due to gravity, exoplanets can also create spiral arms, as in RY Lup (Langlois et al. 2018) and MWC 758 (Benisty et al. 2015) where the presence of exoplanets is predicted by hydrodynamical simulations. Debris disks are the oldest step in the evolution of circumstellar disks, when there is already one or several planets in the system and the gas is almost completely consumed. These disks are composed of dust and grain that never accreted onto the planets, as in the case of HR 4796A (Perrin et al. 2015; Milli et al. 2019).

These environments can be observed with SPHERE in the near-infrared and in the visible. However, such observations are difficult because of the high contrast between the light of the environment and the residual light from the host star. As a result, when acquiring images, the light of the environment is contaminated by the diffraction stains from the host star. Two methods can be used to disentangle the light of the disk from that of the star: angular differential imaging (ADI; Marois et al. 2006) and differential polarimetric imaging (DPI; van Holstein et al. 2020). The ADI technique makes use of the fact that the stellar residuals are fixed in the pupil plane and the object of interest artificially rotates. The resulting diversity makes it possible to disentangle the light of the object of interest from the residual light of the star. Still, such a method does not allow for a good reconstruction of the disk morphology as it is impacted by artifacts due to self-subtraction. Moreover, the method fails when the environment is nearly rotation-invariant. The DPI observations allow for access to the morphology of the disks, without the artifacts, by using the difference of polarization states between the light scattered by the environment and the light of the host star.

The 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 are equipped to deal with the noise statistics. This results in some limitations in sensitivity in the case of a low signal-to-noise ratio (S/N). Lastly, a deconvolution may be performed to get rid of the blurring by the instrumental point spread function. As this is done without accounting for the noise statistics after all the earlier processing, the results are not optimal given the available data. Still, these state-of-the-art methods have proven over the years to be sufficiently efficient to produce quality results. However, studies of circumstellar disks are often limited to analyses of the orientation (position angle and inclination) and morphology (rings, gaps, cavities, and spiral arms) of the disks (Muto et al. 2012; Quanz et al. 2013; Ginski et al. 2016; de Boer et al. 2016a). Quantitative polarimetric measurements of circumstellar disks and substellar companions are currently very challenging because existing data-reduction methods do not properly estimate the sources of the errors from both noise and detector calibration. They also require complete polarimetric cycles and do not account for the instrument convolution. For observations of circumstellar disks (van Holstein et al. 2020), calibrating the instrumental polarization effects with a sufficiently high accuracy has already yielded several improvements.

Over the last decades in image processing, it has been proven that the reconstruction of parameters of interest benefits from a global “inverse problems” approach, taking into account the noise statistics and all instrumental effects, rather than “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 estimation of the physical parameters (Titterington 1985), mostly in adaptive optics (Borde & Traub 2006) as well as 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. 2020) in polarimetric radio-interferometry using more sophisticated tools as non-smooth 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. However, 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 Sect. 2, we develop the physical model of the data obtained with the ESO/VLT SPHERE IRDIS instrument. This includes the polarimetric parametrization, the convolution by the PSF and the observing sequence. In Sect. 3, we describe Reconstruction of High-contrAst Polarized Sources and Deconvolution for near-Infrared Environments (RHAPSODIE1), the method we developed. Section 4 is dedicated to the calibration of the detector, the instrument, and the instrumental polarization. Finally, in Sect. 5, we present the results obtained with RHAPSODIE on both simulated and astrophysical data.

2. Modeling polarimetric data

The principal parameters of interest for studies of circumstellar environments in polarimetry are the intensity, Ip, 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 Iu 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, Iu.

The estimation of the parameters (Iu, Ip, θ) from polarimetric data is the objective of the present contribution. Stokes parameters are, however, more suitable in striving to account for the effects of the instrument on the observable polarization as the model of the data happens to be a simple linear combination of these parameters. Stokes parameters account for the total light, the linearly polarized light, and the circularly polarized light. Since circular polarization is mostly generated by magnetic interactions and double scattering, it is often negligible in the case of circumstellar environments and thus not measured by the SPHERE or GPI instruments. Ultimately, it is possible to reconstruct the parameters of interest Iu, Ip, and θ from a combination of the Stokes parameters.

2.1. Polarization effects

The four Stokes parameters S = (I, Q, U, V) ∈ ℝ4 describe the state of polarized light: I is the total intensity accounting for the polarized and unpolarized light, Q and U are the intensities of the light linearly polarized along 2 directions at 45° to each other, and V is the intensity of the circularly polarized light. Under this formalism, polarization effects by an instrument like SPHERE/IRDIS (see Fig. 1) can be modeled by (see Eq. (17) of van Holstein et al. 2020):

(1)

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

where are the Stokes parameters on the detector after the left (j = 1) or right (j = 2) polarizer of the analyzer set while S represents the Stokes parameters at the entrance of the telescope. In the above equation, M denotes a Mueller matrix accounting for the polarization effects of a specific part of the instrument: for the left or right polarizer of the analyzer set, Mder for the optical derotator, MHWP is for the half-wave plate (HWP), MM4 for the fourth mirror of the telescope and MUT for the three mirrors (M1 to M3) constituting the telescope. The term T(θ) denotes a rotation matrix of the polarization axes by an angle θ: θder is the derotator angle, α is the HWP angle, θalt is the altitude angle, and θpar is the parallactic angle of the pointing of the alt-azimuthal telescope.

In the optical, only the total intensity out of the Stokes parameters can be measured by any existing detector. It then follows from Eq. (1) that the quantity measured by the detector is a simple linear combination of the input Stokes parameters:

(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 , rotating the angle α of the HWP introduces a modulation of the contribution of the Stokes parameters Q and U in , which can be exploited to disentangle the Stokes parameters I, Q and U. The Stokes parameter V characterizing the circularly polarized light cannot be measured with an instrument such as SPHERE/IRDIS (a modulation by a 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:

(3)

where νj, k,  = νj, k) with Θk the set of angles during the kth acquisition.

The instrumental polarization effects being reduced to those caused by the analyzers and the HWP and ignoring the rotation due to the altitude and parallactic angles, the detected intensity writes (van Holstein et al. 2020):

(4)

where ψj is the orientation angle of the left/right polarizer while αk is the HWP angle during the kth acquisition. Table 1 lists the values of the linear coefficients νj, k,  for a typical set of HWP angles.

Table 1.

Positions of HWP αk and orientations of the analyzer ψj and the corresponding values of the coefficients νj, k,  assuming no instrumental polarization and ignoring field rotation.

2.2. Parameters of interest

The model of the detected intensity given in Eq. (3) is linear with regard to the Stokes parameters S = (I, Q, U), which makes its formulation suited to inverse problem solving. However, to study circumstellar environments in polarimetry, the knowledge of the linearly polarized light Ip and the polarization angle θ is crucial. Both set of parameters are related as it follows:

(5)

and conversely by2:

(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 Iu, Ip, and θ:

(7)

For a perfect SPHERE/IRDIS-like instrument and not considering field rotation, combining Eqs. (4) and (5) yields:

(8)

which is the Malus law.

2.3. Accounting for the instrumental spatial PSF

In polarimetric imaging, each polarimetric parameter is a function of the two-dimensional (2D) FOV. We consider that the Stokes parameters are represented by images of N pixels each and denote as S, n the value of the nth pixel in the map of the th Stokes parameter.

Provided that polarization effects apply uniformly across the FOV and that the instrumental spatial point spread function (PSF) does not depend on the polarization of light, all Stokes parameters of a spatially incoherent source are independently and identically affected by the spatial PSF (e.g., Birdi et al. 2018; Smirnov 2011; Denneulin 2020). As the effects of the spatial PSF are linear, we can give the detected intensity for a given detector pixel as:

(9)

where is the intensity measured during the kth acquisition by the mth pixel of the sub-image corresponding to the jth polarizer of the analyzer set. The coefficients νj, k,  (accounting for the instrumental polarization) are defined in Eq. (3) and Hj, 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.

Table 2.

Notations.

It follows from our assumption on spatial and polarization effects applied independently that the discretized spatial PSF does not depend on the polarization index . Consequently, the spatial and polarization effects in Eq. (9) mutually commute.

Figure 1 and Eq. (1) provide a representation of the instrument from which we build a model of the spatial effects of the instrument. Accordingly, an image representing the spatial distribution of the light as the input of the instrument should undergo a succession of image transformations before reaching the detector. These transformations are either geometrical transformations (e.g., the rotation depending on the parallactic angle) or blurring transformations (e.g., by the telescope). Except in the neighborhood of the coronographic mask, the effects of the blur can be assumed to be 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 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 or spatial (de)magnification3. Following this analysis, our model of the spatial PSF is given by:

(10)

where Ak  :  ℝN → ℝN implements the shift-invariant blur of the input model maps while Tj, k  :  ℝN → ℝMj performs the geometrical transform of the blurred model maps for the jth polarizer of the analyzer set during the kth acquisition. N is the number of pixels in the model maps and Mj is the number of pixels of the detector (or sub-image) corresponding to the jth 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:

(11)

where F denotes a fast Fourier transform (FFT) operator of suitable size and implements the frequency-wise multiplication by , the discrete Fourier transform of pk, the shift-invariant PSF. We note that pk must be specified in the same reference frame as the FOV. If the shift-invariant PSF is calibrated from empirical images of, for instance, the host star, acquired by the detector, the inverse (or pseudo-inverse) of Tj, k must be applied to the empirical images.

2.4. Polarimetric data

During a sequence of observations with SPHERE/IRDIS in DPI mode, the HWP is rotated several times along a given cycle of angles α ∈ {0° ,45° ,22.5° ,67.5° }. In addition, the two polarizers of the analyzer set of SPHERE/IRDIS are imaged on two disjoint parts of the same detector. The resulting dataset consists in K frames composed of two (left and right) 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 sub-images, the available data are modeled by:

(12)

where is given in Eq. (9) for k ∈ ⟦1, K⟧ the index of the acquisition, j ∈ {1, 2} indicating the left and right polarizer of the analyzer, and m ∈ ⟦1, M⟧ the pixel index in the corresponding left and right 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 to be 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 of the actual data, the shot noise is the most important contribution and the number of electrons (photo-electrons 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(dj, k, m) is estimated in the calibration stage (see Sect. 4.1).

It is worth noting that since the proposed model is not valid in the neighborhood of the coronographic mask, the data pixels in this region have to be discarded. Moreover, the detector contains defective pixels (e.g., dead pixels, 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 Sect. 3.2.

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

3.1. Inverse problems approach

In polarimetric imaging, one is interested in recovering sampled maps of the polarimetric parameters, which can be the Stokes parameters (I, Q, U) or the intensities of unpolarized and linearly polarized light and the angle of the linear polarization (Iu, Ip, θ) or some mixture of these parameters. In order to remain as general as possible, we denote by X ∈ ℝN × L the set of parameters of interest to be recovered and by X ∈ ℝN, with  ∈ ⟦1, L⟧, the th parameter component which is a 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 as the ones that minimize a given objective function f(X) possibly under constraints expressed as X ∈ 𝒞 with 𝒞 the set of acceptable solutions. The objective function takes the form of the sum of a data-fidelity term fdata(X) and of regularization terms fρ(X):

(13)

where λρ ≥ 0 (∀ρ) are so-called hyperparameters introduced to tune the relative importance of the regularization terms. The data-fidelity term fdata(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 𝒞, e.g., to account for the requirement that intensities are non-negative quantities. These different terms and constraints are detailed in the following sub-sections.

3.2. Data fidelity

With knowledge of the appropriate statistics for the pre-processed data, the 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:

(14)

where denotes the Mahalanobis (1936) squared norm, dj, k = (dj, k, 1, …, dj, k, M) ∈ ℝM collects all the pixels (e.g., in lexicographic order) of the jth sub-image in the kth acquisition as defined in Eq. (12). Similarly, where the terms are given by the model in Eq. (9) applied to the Stokes parameters as a function of the parameters of interest X, S = S(X). For instance, if , then S(X) is obtained by Eq. (5). In the expression of the data fidelity term given by Eq. (14), Wj, 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 Sect. 2.4), we define the diagonal entries of the precision matrix as:

(15)

with Σj, k, m = Var(dj, k, m) the variance of a valid datum dj, k, m. Invalid data include dead pixels, pixels incorrectly modeled by our direct model because of saturation or of the coronograph, missing frames for a given HWP angle, and unusable frames due to overly strong atmospheric turbulence or improper coronograph centering.

We note that in Eq. (14), the Mahalanobis squared norms arise from our Gaussian approximation of the statistics, while the simple sum of these squared norms for each sub-image in each frame is justified by the fact that all frames and all sub-images are mutually independent.

3.3. 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 fdata(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).

3.3.1. Edge-preserving smoothness

We expect that the light distribution of circumstellar environments be mostly smooth with some sharp edges, hence, using an 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 Iu, to the polarized intensity Ip or to the total intensity I by the following regularization terms:

(16)

(17)

(18)

where we denote in boldface sampled maps of polarimetric parameters, for instance the image of the unpolarized intensity or Iu(X) to make explicit that it is uniquely determined by the sought parameters X. In the above expressions, μρ > 0 models the smoothing threshold and Dn : ℝN → ℝ2 is a linear operator which yields an approximation of the 2D spatial gradient of its argument around the nth pixel. This operator is implemented by means of finite differences; more specifically applying Dn to a sampled map u of a given parameter is expressed as:

(19)

where (n1, n2) denotes the row and column indices of the nth pixel in the map. At the edges of the support of the parameter maps, we simply assume flat boundary conditions and set the spatial gradient to zero there.

The regularizations in Eqs. (16)–(18) implement a hyperbolic version of a pseudo-norm of the spatial gradient of a given component of the light distribution which behaves as an L2-norm (i.e., quadratically) for gradients much smaller than μ and as an L1-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-L2 norm of the gradient, as with the last one in Eqs. (16)–(18), yields solutions in which strong changes tend to occur at the same locations in the sets of parameters. In order to encourage sharp edges to occur at the same places in the Stokes parameters Q and U, we also consider using the following regularization for these components:

(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 (Thiébaut 2002). We refer to Denneulin et al. (2019, 2020), Denneulin (2020) for a comparison of possible advanced regularizers.

3.3.2. Tuning of the hyperparameters

In the regularization function fprior(X), the terms defined in Eqs. (16)–(18) and Eq. (20) can be activated (or inhibited) by choosing the corresponding λρ > 0 (or λρ = 0). It is also required to tune the threshold level μρ > 0 in addition to the λρ multipliers. All these hyperparameters have an incidence on the recovered solution: the higher λρ the smoother the corresponding regularized component and a lowering of the threshold, μρ, allows us to capture sharper structures. A number of practical methods have been devised to automatically tune the hyperparameters: Stein’s Unbiased Risk Estimator (SURE: Stein 1981; Eldar 2008; Ramani et al. 2008; Deledalle et al. 2014), generalized cross-validation (GCV: Golub et al. 1979), the L-curve (Hansen & O’Leary 1993), and hierarchical Bayesian strategies (Molina 1994).

Unsupervised tuning of the hyperparameters with GSURE (Eldar 2008), in a prediction error formulation, has been considered in this context. Figure 2 shows the value of GSURE for several reconstructions of RXJ 1615 (Avenhaus et al. 2018). Since it is expected that there are relatively few sharp edges, ∥DnXρ∥≪μρ 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 λρρ:

(21)

thumbnail 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.

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 stud; in addition, for the results presented in this paper, we tuned the hyperparameters by hand by visually inspecting several reconstructions under different settings, as presented in Fig. 2.

The polarized parameter Ip contribute to I, Q, and U (see Eq. (5)). To avoid a contamination of the extracted polarized parameters Ip by the unpolarized component Iu, we regularize Iu and Ip separately. Hence, the regularization of the unpolarized component Iu should be done via defined in Eq. (16) rather than via fI defined in Eq. (18). For the polarized light, the joint regularization of Q and U by fQ + U defined in Eq. (20) is more effective than the regularization of Ip alone by defined in Eq. (17) which does not constrain the angle θ of the linear polarization. In Eq. (13), we therefore take the Set 1 or Set 2 of hyperparameters presented in Table 3. The latter combination is preferable as discussed previously.

Table 3.

Set of hyperparameters used in the present work.

3.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, Iu, and Ip should all be everywhere non-negative.

Since I = Iu + Ip, it is sufficient to require that Iu and Ip be nonnegative. Hence, for the set of parameters , the positivity constraint is expressed as:

(22)

As given for the Stokes parameters X = (I, Q, U), the positivity yields an epigraphical constraint:

(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 (for all pixels n), the positivity of the intensity Ip of the polarized light automatically holds if the parameters X = (Iu, Q, U) are considered. It is then sufficient to impose the positivity of the intensity Iu of the unpolarized light as expressed by the following feasible set:

(24)

3.5. Choice of the polarimetric parameters

Our method expresses the recovered parameters as the solution of a constrained optimization problem specified in Eq. (13). As explained in Sect. 3.3.2, the weights of imposed regularization is chosen via the values of the multipliers λρ.

The constraints can be implemented by the feasible 𝒞 for different choices of the parameters X. More specifically, X = (Iu, Ip. θ), X = (I, Q, U), or X = (Iu, Q, U) can be chosen. Whatever the choice for X, the relations given in Eqs. (5) and (6) can be used to estimate any parameter of interest given the recovered . These relations can also be used to compute the objective function f(X) which require the Stokes parameters needed by the direct model of the data Eq. (9) and various polarimetric component depending on the choice of the regularization.

With X = (I, Q, U), the positivity constraints take the form of an epigraphic constraint that is more difficult to enforce as being not separable in the parameters space. To solve the problem in Eq. (13) with such a constraint, an epigraphic projection is required, leading to the use of a 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 of X = (I, Q, U) may avoid some degeneracies because it ensures the convexity of the problem in Eq. (13). With X = (Iu, Ip, θ) or X = (Iu, 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 requires the provision of the bounds and a numerical function to compute the objective function and its gradient.

In the following, we compare the performances of RHAPSODIE for the polarimetric parameters X = (I, Q, U) and X = (Iu, 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 = (Iu, Q, U). Both parametrizations allow us to access the parameters of interest Ip and θ using Eq. (6). The best choice of polarimetric parameters is then used to process astrophysical datasets.

4. Data calibrations

4.1. Detector calibration

Before the application of a reconstruction method, the calibration of the data is essential to account for the noise and the artifacts linked to the measurement. It allows for the estimation and correction of any pollution induced by the sky background or the instrument as well as the detector behavior in terms of errors on pixel values.

We use an inverse method to calibrate the raw data from these effects: the quantity required for the calibration are jointly estimated from the likelihood of the calibration data direct model (Denneulin 2020). In this method, all the calibration data are expressed as a function of the different contributions (ı.e., flux, sky background, instrumental background, gain, noise, and quantum efficiency). All these quantities are then jointly estimated by the minimization of the 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 (dk)k ∈ ⟦1, K and their weights (Wk)k ∈ ⟦1, K.

4.2. Instrumental calibration

The instrumental calibration is a required step for the estimation, from dedicated data, of the instrumental PSF and of the star centers on each side of the detector. Since the star is placed behind the coronagraphic mask, a simultaneous PSF estimation is not possible. To estimate the PSF, we use a dedicated flux calibration (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 Ak. The PSF model does not include the spiders so that it may remain rotation-invariant (see Fig. 3).

thumbnail 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.

Table 4.

Information of the datasets used for the target reconstruction.

In the case of synthetic observations, the assumed 2D PSF is extracted from the real data RY Lup (e.g., similar to the top left image in Fig. 3). For astrophysical observations, we fit a 2D PSF model on pre-reduced PSF data for each observed target (e.g., Gaussian, Airy, and Moffat fits in Fig. 3).

During the coronographic observing sequence, the star point spread function peak is hidden by the coronagraphic mask and its position was determined using a special calibration (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 the resulting center estimations were averaged. In addition, we used the derotator position and the true north calibration from Maire et al. (2016) to extract the angle of rotation of the north axis. These instrumental calibration steps lead to the estimation of the transformation operators, Tj, k, which rotate and translate the maps of interest to make the centers and the north axes coincide with those in the data.

4.3. Polarization calibration

When the light is reflected by the optical devices in the instrument, some instrumental polarization is introduced, resulting in a loss of polarized intensity and cross-talk contamination between the Stokes parameters Q and U.

The classical method to compensate for this instrumental polarization is to employs the azimuthal Qϕ and Uϕ parameters to reduce the noise floor in the image (Avenhaus et al. 2014). Because this method is limited to 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 and then to perform the instrumental polarization correction. After computing the double difference, IRDAP uses a Mueller matrix model of the instrument carefully calibrated using real on sky data to correct for the polarized intensity Ip (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 Iu 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 correction factors εQ and εU are expressed as follows:

(25)

where NΩ is the number of pixels of the ring Ω. We then compute Qcor = QεQIu and Ucor = UεUIu in order to create an additional set of Qcor and Ucor images by subtracting the measured stellar polarization from the model-corrected Q and U images.

5. Applications on high-contrast polarimetric data

In this section, we compare two RHAPSODIE configurations. The first configuration is defined as “without deconvolution”, meaning that in the data fidelity term (14) does not include the convolution, leading to the simplification of the Eqs. (9) and (10):

(26)

(i.e., Ak is the identity). Such a configuration of RHAPSODIE aims to be comparable to 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 for such a case.

The second configuration is defined as “with deconvolution” and states for the full RHAPSODIE capabilities. Such a configuration of RHAPSODIE is aimed at showing the benefits of the global model compared to an a posteriori deconvolution with the aim of improving the angular resolution.

5.1. Application on synthetic data

The performance of the linear and non-linear methods are first evaluated on synthetic datasets, without (resp. with) the deconvolution displayed in Figs. 46 (resp. Figs. 79). These datasets are composed of unpolarized residual stellar flux, mixed with unpolarized and polarized disk flux. We produce several synthetic datasets following steps given in Appendix A, for different ratios of the polarization of the disk, called τdisk and defined in Eq. (A.1).

thumbnail Fig. 4.

Visual comparison of the reconstructed polarized intensity Ip with the state-of-the-art double difference and the RHAPSODIE methods without deconvolution.

thumbnail Fig. 5.

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

thumbnail Fig. 6.

Comparison of the MSE between the true map and the estimated map of the polarized intensity Ip and of the angle of polarization θ for the double difference and the linear and non-linear RHAPSODIE methods without deconvolution.

thumbnail Fig. 7.

Visual comparison of the reconstructed polarized intensity Ip with the state-of-the-art double difference and the RHAPSODIE methods with deconvolution.

thumbnail Fig. 8.

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

thumbnail Fig. 9.

Comparison of the MSE between the true map and the estimated map of the polarized intensity Ip and of the angle of polarization θ for the double difference and the linear and non-linear RHAPSODIE methods with deconvolution.

The results of the RHAPSODIE methods are compared to the results obtained with the classical double difference method (Tinbergen 2005; Avenhaus et al. 2014). The double difference is applied on recentered and rotated datasets with the bad pixels interpolated. For the comparison with deconvolution, the results of the Double Difference are deconvolved after the reduction. In order to provide fair comparisons with RHAPSODIE, we propose to use an inverse approach rather than applying a high-pass spatial filtering. The deconvolved double difference reconstructions are obtained by solving the following problem:

(27)

where and are the Stokes parameters reconstructed with the double difference and A represents the convolution by the PSF. This deconvolution method performs the deconvolution of Stokes parameters Q and U jointly in order to sharpen the polarized intensity image, but does not recover the polarization signal lost by averaging over close-by polarization signals with the opposite sign and may even introduce artificial structures that are not present in the original source.

For the reconstruction with RHAPSODIE, the hyperparameters of regularization, λρ and μρ, are chosen in order to minimize the total mean square error (MSE), i.e., the sum of the MSE between each estimated parameter , and , obtained from from Eq. (6), and the ground truth , and θgt. The total MSE is given by:

(28)

where, , , and Nθ are the number of pixels with a signal of interest in the respective Iu, Ip, and θ maps. For the deconvolution of the double difference results, and are chosen to minimize only the sum of the MSE on Iu and θ.

For the reconstructions without the deconvolution, Figs. 4 and 5 show that the reconstructions are less noisy with RHAPSODIE. The inner circle is always better reconstructed. Moreover, the non-linear reconstruction (i.e., minimization on Iu, Q and U) is better than the linear when τdisk grows. In both configurations of RHAPSODIE (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 reconstructions with the deconvolution, the RHAPSODIE methods have a better approach to the reconstructions of the polarized intensity, as shown on the reconstructions in Fig. 7, mostly for the inner ring, even if the background is less biased with the proposed deconvolved double difference (as shown on the error maps in Fig. 8). In fact, as shown on the MSE displayed in Fig. 9, for τdisk = 3%, RHAPSODIE delivers a more accurate polarized intensity estimation. The errors of the reconstructions are larger for the double difference than for the RHAPSODIE methods. The RHAPSODIE methods also allow us to achieve a better reconstruction of the angle of polarization, mostly with the non-linear 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 in Fig. 8. Even if according to the MSE in Fig. 9, for τdisk = 25%, the MSE is smaller for the linear RAPSODIE. in fact, the space between the outer ring and the thin ring is better reconstructed with such a configuration (see Figs. 7 and 8).

According to these results, the RHAPSODIE methods are better than the 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.

5.2. Astrophysical data

RHAPSODIE was applied to several IRDIS datasets dedicated to protoplanetary, transition, and debris disks to test the efficiency of our method and to compare it to the state-of-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:

(29)

with where (n1,n2) denotes the row and the column indices of the nth pixel in the map and those of the pixel center. We also present the maps of polarized intensities Ip and the angles of polarization θ estimated by the different methods.

First, the reconstructions of the target TW Hydrae with the double difference and RHAPSODIE without deconvolution are compared in Fig. 10, without and with the correction of the instrumental polarization. The images are scaled by the square of the separation to account for the drop-off of stellar illumination with distance.

thumbnail Fig. 10.

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

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 IRDAP method, these effects are compensated and the disk reconstructed is more accurate (Fig. 10c (iii)). As a result, the contamination of the Uϕ signal from 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 S/N in the area where the signal is low. On the other hand, it is slightly more sensitive to detector flat calibration accuracy. It is worth noticing that IRDAP has a dedicated correction of the detector response nonuniformity (flat field variation between the detector column) for the various amplifiers that is not implemented in RHAPSODIE. However, the flat calibrations of observations more recent than 2016 have improved and do allow for a better calibration and, as a consequence, they do not impact our method efficiency anymore.

The efficiency of the method we have developed is demonstrated on the Figs. 1117. Figures 11 and 12 present the double difference and RHAPSODIE reconstructions of the protoplanery disks TW Hydrae (van Boekel et al. 2017), IM Lupus (Avenhaus et al. 2018) and MY Lupus (Avenhaus et al. 2018). Figures 13, 14, and Fig. 17 present RHAPSODIE reconstructions of transition disks RY Lupus (Langlois et al. 2018) and T Chamaeleontis (Pohl et al. 2017). Figures 13 and 14 also present RHAPSODIE reconstructions of transition disk RXJ 1615 (de Boer et al. 2016a). Figures 15 and 16 present the RHAPSODIE reconstructions of the debris disks HD 106906 (Kalas et al. 2015; Lagrange et al. 2016), HD 61005 (Olofsson et al. 2016), and AU Mic (Boccaletti et al. 2018). These reconstructions are all normalized in contrast to the unpolarized stellar flux (estimated when registering the PSF 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 hyperparameters of regularization and the PSF models used for the reconstructions are listed in Table 5.

thumbnail Fig. 11.

Reconstructions of the polarized intensity Ip 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.

thumbnail 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.

thumbnail Fig. 13.

Reconstructions of the polarized intensity Ip 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.

thumbnail 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.

thumbnail Fig. 15.

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

thumbnail Fig. 16.

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

thumbnail Fig. 17.

Qϕ and Uϕ projections of the reconstructions of RY Lupus and T Chae. The maps are displayed in symmetrical logarithmic scale and normalized in contrast to the unpolarized stellar flux. The coronograph is masked in white.

Table 5.

Information related to the parameters used for the target reconstruction presented in this section.

The contrast between the selected disks polarized light and their host stars unpolarized light responsible for the level of stellar pollution are also very different (ranging from 3 × 10−2 to 8 × 10−6). The achievable contrast is more favorable for highly inclined disk, such as MY Lup in Fig. 11, because in such a case, the star likely shines partially through the disk, which is dimming the starlight and thus decreasing the contrast between the star and the disk.

In all cases, our method produces 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 methods confirms the benefits of our method, but it is more difficult to quantify the gain on real data without knowing the reality; and, thus, this approach should be strengthened through the use of numerical models for these objects. In addition, our method can produce deblurred results, which clearly enhances the angular resolution and is beneficial for interpreting the disk morphology and for studying its physical properties. In all cases, the deconvolution sharpens the polarized intensity image, helps to recover the polarization signal lost by averaging over 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 earlier in this paper, the hyperparameters could be further tuned to adjust the smoothing of the deconvolved images according to the required noise-and-angular resolution trade-off.

For several of these disks (IM Lup, MY Lup in Fig. 11), the outer edge of the disk and thus the lower disk surface have been detected in Avenhaus et al. (2018) and are confirmed by our method. The deconvolution of these datasets allows to further highlight these fainter features, to enhance or reveal the midplane gaps in the case of T Cha (in Fig. 13) which was not identified by Pohl et al. (2017). The reason for this is that without deconvolution, the PSF smears light from the disk upper and lower sides into the midplane gap.

Except for Au Mic, the debris disks presented in Fig. 15 are much fainter in contrast than the protoplanetary or the transition disks in polarized intensity. As a consequence, these datasets have required careful stellar polarization compensation. The HD 106906 debris disk is viewed close to edge-on in polarized light as reported in van Holstein et al. (2020), Esposito et al. (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.

We may remark that the RHAPSODIE deblurred reconstructions of RXJ 1615 (in Fig. 13) and HD 106906 (in Fig. 15) appear noisier than those without deconvolution. This is due to the hyperparameters that were purposely chosen so as to achieve the best angular resolution with the side effect of slight noise amplification. Thus, the noise is not negligible in the reconstruction. As shown in Fig. 2, increasing the regularization weight to smooth the background would lead to a loss of the thin structures.

We also analyzed two datasets (HD 61005 and Au Mic), which were taken under very bad conditions. For these two datasets, our method is also proven to be efficient in using incomplete polarization cycles to recover the disk polarized signal and to deconvolve this signal despite strong artifacts additionally produced by the rotating spiders (i.e., unmasked by the Lyot stop) when observing in field stabilized. The difference in the spider position during the polarimetric cycle results in an artificial polarimetric signal when using standard data reduction techniques. It is worth noticing that the strength of our method to deal with these artifacts comes from its ability to use weighted maps to account for them. For instance, the spiders are weighted by their variance during their rotation frame by frame, in addition to the inclusion of a static bad pixel ponderation. As a result, the contribution of the spiders to the polarimetric signal when using our method is decreased compared to the other methods. The noise created by these spiders remains, as seen in Fig. 15, and can generate artifacts in the deconvolution as seen for HD 61005. Counteracting the artifacts caused by the spiders can be efficiently done by performing DPI observation in pupil tracking mode as proposed in van Holstein et al. (2017). We have also validated the efficiency of RHAPSODIE in such observations which allows a further gain in the reduction of the instrumental artifacts from the telescope spiders.

Another advantage of our method is its ability to use incomplete polarimetric cycles, which are discarded by the classical methods. This capability of our method leads to an increase in S/N, which was quantified more precisely using our model of the data in Denneulin (2020). To benefit from this improvement, the instrumental polarization has to remain small because incomplete polarimetric cycles do not benefit from the instrumental polarization compensation performed by estimating both couples Q and −Q (and U and −U, respectively).

6. Conclusion

In this work, we develop a new method for extracting the polarimetric signal using an inverse problems approach that exploits a model of the measured signal formation process. The method includes a weighted data fidelity term, which takes into account the blur and the polarization due to the instrument, and effectively disentangles polarized signal of interest from stellar contamination. In order to avoid noise amplification in the minimization of the data fidelity term, an edge preserving smoothing penalization was added allowing to favor smooth estimates almost everywhere. The associated minimization problem is solved by standard optimization techniques. Our method enables to accurately measure the polarized intensity and angle of linear polarization of circumstellar disks by taking into account the noise propagation and the observed objects convolution. It has the capability to use incomplete polarimetry cycles (when the instrumental polarization is small), which enhances the sensitivity of these observations. It also takes proper account of bad pixels by using weighted maps instead of interpolating them. These bad pixels can cause systematic errors of several tenths of a percent in the polarization measurements as shown by van Holstein et al. (2021) In addition, the effect of bad pixel interpolation could also have some impact when reaching 0.1% polarimetric accuracy.

We validated the method on both simulated and archive data from SPHERE/IRDIS and compared its performance with state-of-the-art methods. We 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 S/N-and-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 S/N, 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 the scattering phase functions.

RHAPSODIE is the first regularized inverse approach implemented for high-contrast polarimetric imaging. It demonstrates the benefits of advanced signal processing methods in this domain.


1

The code of RHAPSODIE is available online at https://github.com/LaurenceDenneulin/Rhapsodie.jl

2

In this representation, there is a ±180° degeneracy for the linear polarization angle θ.

3

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

Acknowledgments

We acknowledge Rob Van Holstein for his help with the instrumental polarization correction based on careful instrumental calibrations he performed which are implemented in the IRDAP tool he developed. We also thank the anonymous referee for her/his careful reading of the manuscript as well as her/his insightful comments and suggestions. This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon (OSUL/CRAL). This work is supported by the French National Programms PNP and the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU 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.

References

  1. Adam, R., Ade, P. A., Aghanim, N., et al. 2016, A&A, 594, A8 [Google Scholar]
  2. Akiyama, K., Alberdi, A., Alef, W., et al. 2019, ApJ, 875, L4 [NASA ADS] [CrossRef] [Google Scholar]
  3. Avenhaus, H., Quanz, S. P., Schmid, H. M., et al. 2014, ApJ, 781, 87 [Google Scholar]
  4. Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44 [NASA ADS] [CrossRef] [Google Scholar]
  5. Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Berdeu, A., Soulez, F., Denis, L., Langlois, M., & Thiébaut, É. 2020, A&A, 635, A90 [CrossRef] [EDP Sciences] [Google Scholar]
  7. Beuzit, J.-L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Biraud, Y. 1969, A&A, 1, 124 [NASA ADS] [Google Scholar]
  9. Birdi, J., Repetti, A., & Wiaux, Y. 2018, MNRAS, 478, 4442 [CrossRef] [Google Scholar]
  10. Birdi, J., Repetti, A., & Wiaux, Y. 2020, MNRAS, 492, 3509 [Google Scholar]
  11. Boccaletti, A., Sezestre, E., Lagrange, A. M., et al. 2018, A&A, 614, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Borde, P. J., & Traub, W. A. 2006, ApJ, 638, 488 [NASA ADS] [CrossRef] [Google Scholar]
  13. Carrillo, R. E., McEwen, J. D., & Wiaux, Y. 2012, MNRAS, 426, 1223 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  14. Charbonnier, P., Blanc-Féraud, L., Aubert, G., & Barlaud, M. 1997, IEEE Trans. Image Process., 6, 298 [Google Scholar]
  15. Chierchia, G., Pustelnik, N., Pesquet-Popescu, B., & Pesquet, J.-C. 2014, IEEE Trans. Image Process., 23, 5531 [Google Scholar]
  16. Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, in Ground-based and Airborne Instrumentation for Astronomy II, Int. Soc. Opt. Photonics, 7014, 70143E [Google Scholar]
  17. Combettes, P., & Pesquet, J. C. 2011, Fixed-point algorithms for inverse problems in science and engineering, (Springer) 185 [Google Scholar]
  18. Combettes, P. L., & Wajs, V. R. 2005, Multiscale Model. Simul., 4, 1168 [CrossRef] [Google Scholar]
  19. de Boer, J., Salter, G., Benisty, M., et al. 2016, A&A, 595, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. de Boer, J., Langlois, M., van Holstein, R. G., et al. 2020, A&A, 633, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Deledalle, C.-A., Vaiter, S., Fadili, J., & Peyré, G. 2014, SIIMS, 7, 2448 [Google Scholar]
  22. Denneulin, L. 2020, PhD Thesis, Lyon, France [Google Scholar]
  23. Denneulin, L., Langlois, M., Pustelnik, N., & Thiébaut, E. 2019, in GRETSI LILLE 2019 [Google Scholar]
  24. Denneulin, L., Pustelnik, N., Langlois, M., Loris, I., & Thiébaut, E. 2020, iTwist, Nantes, France, Dec. 2–4, 2020 [Google Scholar]
  25. Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73 [NASA ADS] [CrossRef] [Google Scholar]
  26. Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in Ground-based and Airborne Instrumentation for Astronomy II, Int. Soc. Opt. Photonics, 7014, 70143L [Google Scholar]
  27. Eldar, Y. C. 2008, IEEE Trans. Signal Process., 57, 471 [Google Scholar]
  28. Esposito, T. M., Kalas, P., Fitzgerald, M. P., et al. 2020, AJ, 160, 24 [Google Scholar]
  29. Flasseur, O., Denis, L., Thiébaut, É., Olivier, T., & Fournier, C. 2019, in 2019 27th European Signal Processing Conference (EUSIPCO) (IEEE), 1 [Google Scholar]
  30. Garufi,, A., Benisty, M., Stolker, T., et al. 2017, The Messenger, 169, 32 [Google Scholar]
  31. Ginski, C., Stolker, T., Pinilla, P., et al. 2016, A&A, 595, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Golub, G. H., Heath, M., & Wahba, G. 1979, Technometrics, 21, 215 [CrossRef] [Google Scholar]
  33. Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
  34. Hansen, P. C., & O’Leary, D. P. 1993, SISC, 14, 1487 [Google Scholar]
  35. Högbom, J. 1974, A&AS, 15, 417 [Google Scholar]
  36. Kalas, P. G., Rajan, A., Wang, J. J., et al. 2015, ApJ, 814, 32 [NASA ADS] [CrossRef] [Google Scholar]
  37. Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Lagrange, A. M., Langlois, M., Gratton, R., et al. 2016, A&A, 586, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Langlois, M., Dohlen, K., Vigan, A., et al. 2014, in Ground-based and Airborne Instrumentation for Astronomy V, Int. Soc. Opt. Photonics, 9147, 91471R [Google Scholar]
  41. Langlois, M., Pohl, A., Lagrange, A.-M., et al. 2018, A& A, 614, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Lefkimmiatis, S., Ward, J. P., & Unser, M. 2013, IEEE Trans. Image Process., 22, 1873 [Google Scholar]
  43. Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proc. Natl. Acad. Sci., 111, 12661 [NASA ADS] [CrossRef] [Google Scholar]
  44. Mahalanobis, P. C. 1936, National Institute of Science of India [Google Scholar]
  45. Maire, A. L., Langlois, M., Dohlen, K., et al. 2016, in Ground-based and Airborne Instrumentation for Astronomy VI, Int. Soc. Opt. Photonics, 9908, 990834 [Google Scholar]
  46. Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
  47. Milli, J., Engler, N., Schmid, H. M., et al. 2019, A&A, 626, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Molina, R. 1994, IEEE Trans. Pattern Anal. Mach. Intell., 16, 1122 [Google Scholar]
  49. Muto, T., Grady, C., Hashimoto, J., et al. 2012, ApJ, 748, L22 [NASA ADS] [CrossRef] [Google Scholar]
  50. Nocedal, J., & Wright, S. J. 1999, Numerical Optimization, Springer Series in Operations Research (New York: Springer) [Google Scholar]
  51. Olofsson, J., Samland, M., Avenhaus, H., et al. 2016, A&A, 591, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Pairet, B., Jacques, L., & Cantalloube, F. 2019, Proceedings of SPARS’19, 1, 1 [Google Scholar]
  53. Perrin, M. D., Duchene, G., Millar-Blanchaer, M., et al. 2015, ApJ, 799, 182 [Google Scholar]
  54. Pinte, C., van der Plas, G., Menard, F., et al. 2019, Nat. Astron., 3, 1109 [Google Scholar]
  55. Pohl, A., Sissa, E., Langlois, M., et al. 2017, A&A, 605, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Price, D. J., Cuello, N., Pinte, C., et al. 2018, MNRAS, 477, 1270 [Google Scholar]
  57. Pustelnik, N., Benazza-Benhayia, A., Zheng, Y., & Pesquet, J.-C. 2016, Wiley Encyclopedia of Electrical and Electronics Engineering (Hoboken, NJ, USA: John Wiley & Sons, Inc.), 1 [Google Scholar]
  58. Quanz, S. P., Avenhaus, H., Buenzli, E., et al. 2013, ApJ, 766, L2 [NASA ADS] [CrossRef] [Google Scholar]
  59. Ramani, S., Blu, T., & Unser, M. 2008, IEEE Trans. Image Process., 17, 1540 [Google Scholar]
  60. Rudin, L. I., Osher, S., & Fatemi, E. 1992, Physica D: Nonlinear Phenomena, 60, 259 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  61. Schmid, H. M., Bazzon, A., Roelfsema, R., et al. 2018, A&A, 619, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Sissa, E., Olofsson, J., Vigan, A., et al. 2018, A&A, 613, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Smirnov, O. M. 2011, A&A, 527, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Starck, J.-L., Donoho, D. L., & Candès, E. J. 2003, A&A, 398, 785 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Stein, C. M. 1981, Ann. Stat., 1135 [Google Scholar]
  66. Tarantola, A. 2005, Inverse Problem Theory and Methods for Model Parameter Estimation (SIAM) [Google Scholar]
  67. Thiébaut, E. 2002, in Astronomical Data Analysis II, Int. Soc. Opt. Photonics, 4847, 174 [Google Scholar]
  68. Thiébaut, E., & Conan, J.-M. 1995, JOSA A, 12, 485 [Google Scholar]
  69. Tikhonov, A. N. 1963, Sov. Math. Dokl. [Google Scholar]
  70. Tinbergen, J. 2005, Astronomical Polarimetry (Cambridge University Press), google-Books-ID: SAS4JzAaMxkC [Google Scholar]
  71. Titterington, D. M. 1985, A&A, 144, 381 [NASA ADS] [Google Scholar]
  72. van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132 [NASA ADS] [CrossRef] [Google Scholar]
  73. van Holstein, R. G., Snik, F., Girard, J. H., et al. 2017, Proc. SPIE, 10400, 1040015 [Google Scholar]
  74. van Holstein, R. G., Girard, J. H., de Boer, J., et al. 2020, A&A, 633, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. van Holstein, R. G., Stolker, T., Jensen-Clem, R., et al. 2021, A&A, 647, A21 [EDP Sciences] [Google Scholar]
  76. Vigan, A., Langlois, M., Dohlen, K., et al. 2014, in Ground-based and Airborne Instrumentation for Astronomy V, Int. Soc. Opt. Photonics, 9147, 91474T [Google Scholar]

Appendix A: Synthetic dataset simulation

In order to evaluate the performance of the RHAPSODIE method, synthetic data have been created. These synthetic data are designed to reproduce astrophysical cases. First, the truth N = 128 × 128 maps , et θgt are created. Such a value of N pixels fits the main Region Of Interest (ROI) size.

The hardest disk structures to reconstruct are faint, small, and lightly polarized structures, and consequently have high contrast with the unpolarized stellar intensity and a low S/N. This is why the synthetic environment we generated a disk with three rings of equal brightness but with a different contrast with the unpolarized stellar intensity. This disk is partially polarized with a linearly polarization Ip and a polarization angle θ ∈ ] − π, π] and an unpolarized component Iu disk. The disk polarization ratio between both component is given by:

(A.1)

These synthetic images are then combined as maps of Stokes parameters I, Q, and U.

The Iu disk component is mixed with the unpolarized Iu star stellar components and a point source companion (star close to the host star of different brightness). The unpolarized intensity is represented as Iu = Iu star + Iu disk. It is important to keep in mind that unmixing both the disk and stellar unpolarized components is not possible from the DPI data without the diversity introduced by ADI. The τdisk value used to synthetise datasets is thus inaccessible in practice from observational polarimetric datasets. To assert the case difficulty, one can then use the total polarization, given for all pixel n ∈ {1, …, N} by:

(A.2)

and the S/N, given for all pixel n ∈ {1, …, N} by:

(A.3)

where is the read-out noise variance. The difference between τtotal and τdisk is that the last one does not take into account the unpolarized star residuals. Figure A.1 presents the S/N maps and the maps of total polarization ratio of the synthetic parameters generated for different τdisk. At the center, where the unpolarized star residual are the brightest, the S/N and the total polarization ratio are the weakest, especially in the case of small Iu disk. Yet the S/N grows with the separation from the star center (as with the stellar contribution or when the polarized contribution of the disk increases). Figure A.2 present the true simulated maps for τdisk = 10%.

thumbnail Fig. A.1.

Total ratios of polarization τtotal maps, with respect to the total intensities, and S/N maps for the different values of disk polarized fractions τdisk.

thumbnail Fig. A.2.

Schematic illustration describing the process of the data simulation. Starting from synthetic maps , , and θgt based on a disc model, we generate artificial calibrated and pre-processed datasets d and associated weights W, illustrated here for τdisk = 10%.

Finally, to generate synthetic calibrated data, the Stokes maps are combined following the expression of the data physical model (12), with K × M noise realizations from a fixed random seed. The values of 10% of the pixels, chosen at random, are replaced by zeros to mimic bad pixel behaviours. The weights related to each acquisition are simulated at the same time following Eq. (15). The datasets are composed of 8 HWP cycles, with two acquisitions per positions in each cycles, leading to 16 images per position of the half-wave plate (HWP). Thus, the total number of images per datasets is K = 64.

Several datasets were created with τdisk ∈ {3%, 7%, 10%, 15%, 25%}, corresponding to difficult cases for τdisk < 10% and less difficult brighter cases above this threshold. Before producing a K × 2N noise realization for each dataset, the random seed is reset to the same value. This allow the reproducibility of the results. In fact this realization is obtained by the multiplication of the standard deviation of the pixel to a Gaussian, centered and reduced gaussian realization. Since the random seed is the same for each dataset, the realization is the same for the given pixel, only the standard deviation changes.

In order to compare the results of the double difference to the results of the RHAPSODIE methods the dataset are pre-processed. The bad pixels are interpolated; then the left and right part of the images are cut, recentered, and rotated.

All Tables

Table 1.

Positions of HWP αk and orientations of the analyzer ψj and the corresponding values of the coefficients νj, k,  assuming no instrumental polarization and ignoring field rotation.

Table 2.

Notations.

Table 3.

Set of hyperparameters used in the present work.

Table 4.

Information of the datasets used for the target reconstruction.

Table 5.

Information related to the parameters used for the target reconstruction presented in this section.

All Figures

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

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail Fig. 4.

Visual comparison of the reconstructed polarized intensity Ip with the state-of-the-art double difference and the RHAPSODIE methods without deconvolution.

In the text
thumbnail Fig. 5.

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

In the text
thumbnail Fig. 6.

Comparison of the MSE between the true map and the estimated map of the polarized intensity Ip and of the angle of polarization θ for the double difference and the linear and non-linear RHAPSODIE methods without deconvolution.

In the text
thumbnail Fig. 7.

Visual comparison of the reconstructed polarized intensity Ip with the state-of-the-art double difference and the RHAPSODIE methods with deconvolution.

In the text
thumbnail Fig. 8.

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

In the text
thumbnail Fig. 9.

Comparison of the MSE between the true map and the estimated map of the polarized intensity Ip and of the angle of polarization θ for the double difference and the linear and non-linear RHAPSODIE methods with deconvolution.

In the text
thumbnail Fig. 10.

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

In the text
thumbnail Fig. 11.

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

In the text
thumbnail Fig. 12.

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

In the text
thumbnail Fig. 13.

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

In the text
thumbnail Fig. 14.

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

In the text
thumbnail Fig. 15.

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

In the text
thumbnail Fig. 16.

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

In the text
thumbnail Fig. 17.

Qϕ and Uϕ projections of the reconstructions of RY Lupus and T Chae. The maps are displayed in symmetrical logarithmic scale and normalized in contrast to the unpolarized stellar flux. The coronograph is masked in white.

In the text
thumbnail Fig. A.1.

Total ratios of polarization τtotal maps, with respect to the total intensities, and S/N maps for the different values of disk polarized fractions τdisk.

In the text
thumbnail Fig. A.2.

Schematic illustration describing the process of the data simulation. Starting from synthetic maps , , and θgt based on a disc model, we generate artificial calibrated and pre-processed datasets d and associated weights W, illustrated here for τdisk = 10%.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.