Revealing asymmetrical dust distribution in the inner regions of HD 141569

We obtained polarimetric differential imaging of a gas-rich debris disk around HD 141569A with SPHERE in the H-band to compare the scattering properties of the innermost ring at 44 au with former observations in total intensity with the same instrument. In polarimetric imaging, we observed that the intensity of the ring peaks in the south-east, mostly in the forward direction, whereas in total intensity imaging, the ring is detected only at the south. This noticeable characteristic suggests a non-uniform dust density in the ring. We implemented a density function varying azimuthally along the ring and generated synthetic images both in polarimetry and in total intensity, which are then compared to the actual data. We find that the dust density peaks in the south-west at an azimuthal angle of $220^{\circ} \sim 238^{\circ}$ with a rather broad width of $61^{\circ} \sim 127^{\circ}$. Although there are still uncertainties that remain in the determination of the anisotropic scattering factor, the implementation of an azimuthal density variation to fit the data proved to be robust. Upon elaborating on the origin of this dust density distribution, we conclude that it could be the result of a massive collision when we account for the effect of the high gas mass that is present in the system on the dynamics of grains. Using the outcome of this modelization, we further measured the polarized scattering phase function for the observed scattering angle between 33$^{\circ}$ and 147$^{\circ}$ as well as the spectral reflectance of the southern part of the ring between 0.98 $\mu$m and 2.1 $\mu$m. We tentatively derived the grain properties by comparing these quantities with MCFOST models and assuming Mie scattering. Our preliminary interpretation indicates a mixture of porous sub-micron sized astro-silicate and carbonaceous grains.


Introduction
The formation of planetary systems is a complex process involving several stages that are not yet fully understood. When the gas from the primordial disk is dissipated and planets are formed, typically after 5-10 Myrs, the remaining mass of the system is distributed in planetesimals, which, due to the gravitational interplay with planets, are usually arranged in the form of one or several rings orbiting as far as tens of astronomical units (au) from the star (depending on the stellar mass), akin to the asteroid and Kuiper belts in our own Solar System. The current paradigm implies that planetesimals in these rings undergo mutual destructive collisions, initializing a collisional cascade that produces fragments of various sizes, all the way down to small dust grains (Wyatt 2008). These small dust grains at the tail end of the collisional cascade can be detected in the form of a photometric excess in the infrared, but often in imaging as well. Sys-tems in which such second generation fragment production is underway are referred to as debris disks.
A significant fraction of these debris disks have been imaged, from the visible up to the millimetric domains, revealing a great variety of spatial structures, such as arms, spirals, radial gaps, warps, clumps, and two-sided asymmetries . Several scenarios have been explored for the possible origin of these structures, such as dynamical interactions with planets (e.g., Beust et al. 2009;Thebault et al. 2012;Nesvold & Kuchner 2015), interactions with leftover or second-generation gas (Takeuchi & Artymowicz 2001;Lyra & Kuchner 2013), or dynamical perturbations of companion stars (Thébault 2012). One scenario that has received a lot of attention in the past decade involves giant impacts between Ceres-to-planet-sized embryos, which could create transient bright spots or arms (Thebault & Kral 2018), but also more long-lived structures. The typical longlived outcome of a large-scale collision is an asymmetric struc-ture with a bright clump at the location of the initial breakup, through which all debris must pass, and a more extended and diffuse structure on the opposite side (Jackson et al. 2014;Kral et al. 2015). The survival time of these asymmetries can be of the order of 10 5 up to several 10 6 years.
Observations in direct imaging at near-infrared (NIR) or optical wavelengths are sensitive to the stellar light scattered by the dust. The way the light is scattered depends on the disk morphology (size and inclination, in particular), grains properties (composition, shape, and size), and surface density. As a result, direct imaging has the ability to provide constraints on these characteristics. The latest facility instruments, Spectro-Polarimetric Highcontrast Exoplanet REsearch (SPHERE, Beuzit et al. 2019), Gemini Planet Imager (GPI, Macintosh et al. 2014), and Subaru Coronagraphic Extreme Adaptive Optics project (SCExAO, Jovanovic et al. 2015) have now achieved an impressive record of disk studies which have been efficiently complemented by Atacama Large Millimeter/submillimeter Array (ALMA) results (e.g., MacGregor et al. 2018;Olofsson et al. 2019).
In the evolution of planetary systems, the disk around the Herbig Ae/Be star HD 141569A (∼ 5 Myr, Weinberger et al. 2000;Merín et al. 2004) occupies a unique place as it shares the characteristics of both protoplanetary and debris disks, and is regarded as a turning point object between these two categories. Initially classified as a transitional disk, it is now considered a gas-rich debris disk. Studying this object provides the opportunity to shed light on a very specific moment in the evolution of planetary systems. With favorable characteristics in proximity (110.63 +0.91 −0.88 pc, Gaia Collaboration 2018) and size (400 au), this unique system has been the subject of intensive studies, the objective being to reach the innermost regions where planets are expected to form.
The disk of HD 141569A was first resolved in the NIR coronagraphic observations with HST/NICMOS2 at 1.6 µm (Augereau et al. 1999a) and at 1.1 µm (Weinberger et al. 1999), allowing for two ring-like components to be resolved, located at roughly 250 au and 410 au from the star. Follow-up observations at visible wavelengths laid out a detailed view of the structure of these rings composed of arcs and spirals (Mouillet et al. 2001;Clampin et al. 2003), later confirmed in the NIR with an 8-m class ground-based telescope (Biller et al. 2015;Mazoyer et al. 2016).
Based on the spectral energy distribution (SED), the presence of material nested within the two outer rings, inwards of 100 au, was suggested by Augereau et al. (1999a) and then established owing to mid-IR thermal imaging (Fisher et al. 2000;Marsh et al. 2002;Thi et al. 2014). More recently, Konishi et al. (2016) and Perrot et al. (2016), with HST/STIS and SPHERE, respectively, identified substructures in the central area in scattered light. In particular, SPHERE provided a unique view of the inner system by revealing, for the first time, several rings the most prominent being located at about 45 au. This finding was confirmed with L band imaging by Mawet et al. (2017) and in polarimetry with GPI (Bruzzone et al. 2020). The innermost ring is of great interest because it features a strong north-south brightness asymmetry along the major axis of the disk, whereas we would rather expect an east-west asymmetry due to scattering according to the disk inclination (∼ 57 • ). For this reason, Perrot et al. (2016) suggested that the density of the innermost ring is not constant azimuthally, which has not, thus far, been accounted for in modeling works (Bruzzone et al. 2020). As for the gas content, the most detailed observations obtained at a high angular resolution with ALMA ( 13 CO J = 1 → 0 transition acquired at a spatial resolution of 0.76 × 0.56 ) allowed us to assess the CO gas distribution at similar distances than in scattered light (Di Folco et al. 2020). The corresponding map is asymmetrical with a significant excess of emission in the western part of the disk (from position angles of 220 • to 290 • ).
The purpose of the analysis presented in this paper is to take advantage of both the polarimetric and total intensity data obtained with SPHERE in order to solve for the density and scattering functions. Section 2 presents the observations and data reduction methods. The morphology of the disk and the modeling of the innermost ring both in polarized and total intensities are explained in Sects. 3 and 4, respectively. In Sects. 5 and 6, we present the methods for extracting the phase function of the ring, followed by the spectroscopy of the southern part of the ring in the total intensity data. In Sect. 7 we model the grain properties. In Sect. 8 we discuss the plausible physical origin of the ring.
The total intensity data has been acquired in May 2015 using the dual band imaging (DBI, Vigan et al. 2010) mode of IRDIS with filters H2 and H3, and with the Integral Field Spectrograph (IFS, Claudi et al. 2008) in Y J mode (0.95 − 1.35µm, R ∼ 54). More details of the observations can be found in Perrot et al. (2016). The second and third epochs were obtained in March 2016 and June 2017 with a different setup (IRDIFS-ext), which combines IRDIS in the K1 and K2 filters and IFS in Y JH mode (0.95 − 1.55 µm, R ∼ 33). All the data were processed using the data reduction and handling (DRH) pipeline (Pavlov et al. 2008;Mesa et al. 2015) at the SPHERE Data Centre (DC) 1 (Delorme et al. 2017). This preliminary data-reduction process included sky and dark subtractions, flat-field corrections, star-centering using calibration spots, distortion correction (Maire et al. 2016), bad-pixel removal, and wavelength calibration. The point-spread function (PSF) was obtained by moving the star outside of the coronagraphic mask and placing a neutral density (ND2) in the beam. The data are processed with Karhunen-Loève Image Projection (KLIP, Soummer et al. 2012) to perform advanced angular differential imaging (ADI, Marois et al. 2006).
The polarimetric data, on the other hand, were reduced with IRDAP, the IRDIS Data reduction for Accurate Polarimetry (van Holstein et al. 2020, version 1.3.0). IRDAP addresses the crosstalk and the Instrumental Polarisation (IP) effects of both the telescope and IRDIS. The pipeline first pre-processesed the Notes. The columns from left to right indicate the date of observations in UT, the filter, the field rotation, the number of polarimetric cycles, the individual exposure time (DIT) in seconds, the total number of exposures, the total exposure time in seconds, the DIMM seeing measured in arcseconds, the correlation time τ 0 in milliseconds, the variation of the flux during the sequence in %, and the true north (TN) offset in degrees.
data by applying dark subtraction, flat fielding, bad-pixel correction, and frame centering. It then computes the double-difference and double-sum images (Tinbergen 1996) to obtain the linear Stokes parameters (Q and U) and the corresponding total intensities (I Q and I U ), which are later corrected for IP effects and crosstalk. IRDAP provides the linearly polarized intensity, angle, and degree of linearly polarized light of the source, and computes azimuthal Stokes parameters Q φ and U φ (de Boer et al. 2020), with or without the star polarization removed. This is to say that IRDAP can identify the polarization signal of the star that may be present in the images possibly due to a spatially unresolved inner disk component. Such a stellar polarization signal creates a halo of polarized light that produces a butterfly pattern in the Q φ and U φ images, thereby distorting the images at close separations from the star. Throughout the paper, we use the star polarization subtracted images. Under the assumptions of single scattering and low inclinations (Canovas et al. 2015), Q φ thus represents the polarimetric signal, and U φ the noise.

Morphology of the disk
This section describes the H−band polarimetric image of the disk around HD 141569A, and provides a comparison with the total intensity data by Perrot et al. (2016). Figure 1a displays the polarimetric Q φ image in a 5 ×5 field of view, which is smoothed with a five-pixel (61 mas compared to the diffractionlimited resolution of 42 mas) Gaussian kernel. While the previously known outer ring at 250 au is barely detected, we can still appreciate an east-west brightness asymmetry that is in apparent contradiction with the total intensity data, where this ring is more homogeneous. None of the fine structures identified in this ring by Perrot et al. (2016) are visible in Fig. 1a due to the poor signal to noise ratio (S/N). The even fainter outer most ring at 400 au is not visible in the DPI data for the same reason. In DPI, the transmission of the instrument is reduced as compared to the total intensity observations due to the use of the IRDIS polarizers and also because only a fraction of the light is actually polarized. The black and white circular marks in the lower left corner of Fig. 1a are caused by clusters of bad pixels. 2 The zoom-in Q φ and U φ images in Fig. 1b and c show the inner 2 ×2 region of the system, which are smoothed similarly as in Fig. 1a. In the Q φ image, we can recognize presumably the main components discovered in Perrot et al. (2016), namely the ringlets R1, R2, R3 as labelled in the zoom-in total intensity image (Fig. 1d), which are also visible in its corresponding S/N map in Fig. A.1. Starting from the shortest distances, R3 is located at around 45 au. In Fig. 1b, R3 features an interesting intensity distribution, with a brighter lobe in the southeast, and becomes almost undetected in the western side. R2 is a sort of hook-like structure apparently emanating from the northeast side of the inner ring R3 and extending upwards to the north. It can be confused with the northern tip of R3 but is actually at a larger distance. A portion of R2 is likely also visible in the southeast. R1 could be an arc, a portion of a ring, or a spiral structure extending in the south-east direction.
In the IFS images, the western side of R2 and the southern and the north-western part of R3 are clearly visible at all wavelengths. The reduced images with their corresponding S/N maps are presented in Figs. B.1 and B.2, respectively. To reinforce the identification of the rings, we took advantage of the three epochs and followed the procedure described in Gratton et al. (2019) in which the data are processed with both angular and spectral differential imaging (ASDI-PCA, Mesa et al. 2015), and their corresponding S/N maps are multiplied. The result shown in Fig. C.1 definitely confirms the detection of R2 and gives a stronger evidence of the north-south asymmetry in R3. However, we refrain from deriving any physical parameters from the outcome of this procedure given that the intensity is not well preserved.
In the IRDIS image as shown in Fig. 1d, a bright blob is seen in the northeast at a distance of 387 mas with a diameter of 90 mas. This blob is not detected in any of the IFS wavelength channels (Fig. B.1), nor in the IRDIS K band images (Perrot et al. 2016). If it were a dust feature, it would also be visible in the polarised intensity map, however, this is not the case. Moreover, it is aligned with the wind direction. Altogether, these arguments lead us to conclude that this source is an artifact superimposed over R3. In addition, the clump-like feature reported at the southern ansa by Perrot et al. (2016) is detected with both IRDIS (Fig. 1d) and IFS in total intensity (Fig. B.1). On the contrary, the absence of this feature in the Q φ image (Fig. 1b) could indicate an ADI-induced artifact (Milli et al. 2012).
While both DPI and ADI images reveal material at the very same location, the polarized and total intensities diverge strongly, owing to a combination of phase function effects, ADI self-subtraction effects, and, possibly, density effects. Focusing on the ring R3, the total intensity is predominantly larger in the southern part. The minor axis of this ring is strongly af-Article number, page 3 of 21 A&A proofs: manuscript no. main fected by self-subtraction as a result of a moderate field rotation (42.7 • ) and short angular separations (0.41 ). As for the case of HR 4796, any forward scattering peak expected toward the east would be damped by the self-subtraction (Milli et al. 2017). On the contrary, the polarized image reveals the disk in a larger range of position angles. The polarized intensity of the ring (see details in Sect. 4.3) is increasing from the north toward the east, and is reaching an apparent plateau between PA ∼115 • and ∼180 • . Then, the intensity drops again toward the west side while the disk becomes undetectable. Therefore, a natural as-sumption would be that the density of R3 is larger in the southeast, but this is at odds with respect to the total intensity image, which shows a nearly east-west symmetrical intensity in the southern part.
The contours displayed in Fig. 2a represent the elliptical fits performed on the ADI image by Perrot et al. (2016, Table A.2) to measure the geometrical parameters of R1, R2, and R3. When overlaid in the polarimetric Q φ image in Fig. 2b, the elliptical fits differ mostly in the northern part, where they are largely unconstrained in the former total intensity observations (Fig. 2a). Also, the ringlets do not look as sharp as in total intensity. This is likely due to the ADI technique itself, which is known to impact the local shape and photometry of extended objects (Milli et al. 2012). DPI is not affected by such a bias and, in this particular case, it is certainly more reliable in constraining the observed morphology.
To better constrain the geometry of the ring R3, we isolated the structure by imposing an aperture on the Q φ image encompassing only R3. We used a full elliptical mask as an aperture with semi-major and semi-minor axes of 0.35 and 0.18 for the inner ellipse and 0.48 and 0.27 for the outer ellipse as shown in Fig. 3 in white dashes. The same aperture is used in Sect. 4 in the context of the modeling. The spine of the ring is measured as for instance in Boccaletti et al. (2013) by applying a Gaussian fit to each direction (steps of 1 • ). We used a non-linear least squares algorithm to fit an non-centered ellipse, and obtained the following parameter values: i = 56.5 ± 0.7 • , PA = 356.8 ± 0.4 • , x off = 0.2 ± 0.2 au, y off = 0.6 ± 0.2 au. The inclination and PA are in close agreement with Perrot et al. (2016).
These results are broadly consistent with previous observations with GPI. Bruzzone et al. (2020) reported a diffuse ring between 0.4 and 0.9 in the H-band GPI polarized image, featuring a global east-west asymmetry as expected from the forward scattering of sub-micron sized grains. They also report the observed south-eastern brightness difference that is seen in Fig. 2b, particularly obvious after the subtraction of a symmetrical model (Fig. 8, Bruzzone et al. 2020). They interpreted this excess flux in the South as an arc-like over-density along the southern part, in reminiscence of the other arc-like structures observed at larger separation (Clampin et al. 2003). In their modeling approach, they did not consider the ringlets individually. As a result, they derived rather unconstrained values for the inclination (60 ± 10 • ) and PA (5 ± 10 • , dashed light yellow ellipse in Fig. 2b), which stand in contrast with the accuracy of our measurements to a certain extent.

General assumptions
It has been suggested in Perrot et al. (2016) that the observed intensity distribution of R3 ( Fig. 1d) cannot be explained only by the light scattering properties of the dust of a homogeneous ring-like disk given the east-west inclination, thus arguing for an azimuthal variation of the dust density.
In the polarimetric image ( Fig. 1b), most of the disk intensity is located to the east, which confirms that the eastern side is the front side of the disk, if the scattering is predominantly forward. The lack of polarized and total intensity to the north further indicates probable depletion of dust. Thus, the DPI observation complements the ADI data in breaking the degeneracy due to the phase function of the dust, and confirms that the dust density of R3 is likely not azimuthally constant.
We used GRaTeR, a radiative transfer code (Augereau et al. 1999b), to generate synthetic scattered-light images of the ring R3, which we assumed to be optically thin (Thi et al. 2014). The scattering phase function (SPF), which determines the amount of light scattered at a particular scattering angle (θ) for total intensity observation is given by the Henyey Greenstein (HG, Henyey & Greenstein 1941) function ( f I,θ ). We consider Rayleigh approximation to account for the angular dependence of the linear polarization while modeling our DPI data. It has been experimentally demonstrated that the degree of linear polarization (DOLP) as a function of scattering angles for zodiacal and cometary grains has a bell shape (e.g., Leinert et al. 1976;Bertini et al. 2017;Frattin et al. 2019), which is well approximated by the Rayleigh scattering function. Since debris disks are expected to be reservoirs of cometary grains, the assumption of representing the polarized SPF with a model following a linear combination of HG and Rayleigh approximation was adopted in the analysis of HIP79977 (Engler et al. 2017) and HR 4796 ). Thus, for our data as well, we assume that the polarized SPF is expressed by: Previous works (Perrot et al. 2016;Bruzzone et al. 2020) have assumed that the density function varies radially (R) and vertically (Z). Here, we considered an additional component (A) to parameterize the azimuthal variation. The density function n is expressed as follows: where n 0 is the peak of the density function, r the radial dimension in the disk plane, z the vertical dimension perpendicular to the disk plane, and ϕ the azimuthal dimension in the disk plane.
In addition, we considered two approaches to explore the parameter space: grid-search and Monte-Carlo Markov Chain (MCMC, emcee package Foreman-Mackey et al. 2013) to derive the best model that could explain the observed data. This is to ensure that we obtain the best solution with well-constrained posterior probability distributions and also to compare the reliability of these widely used methods.

Azimuthal density distribution
In the effort to account for the origin of azimuthally variable dust density, the two leading scenarios, namely, the collisional breakup of a large planetesimal as well as the interaction with gas (further discussed in Sect. 8), make relatively similar predictions regarding the azimuthal geometry of the ring, at least from a qualitative point of view. Jackson et al. (2014) have indeed shown that in the aftermath of a planetesimal breakup, the dusty ring that forms peaks at the location of the breakup, remains relatively localized in an angular domain (whose width depends on the size of the parent body) around this peak, and is minimum on the opposite side. As for the gas-related scenario, we would expect the dust to follow, to a certain extent, the azimuthal profile of the gas, which for the inner regions of the HD 141569 system, has been shown to be concentrated in a particular range of angular separations and position angles around this peak (Di Folco et al. 2020).
In principle, the density distribution associated to the afore mentioned scenarios should be modelled in 2D, both in the radial and azimuthal directions, but the S/N of the polarimetric image does not allow us to reliably explore the radial density distribution beyond the region where the R3 ring peaks. Thus we took a conservative approach and considered a 1D prescription. As a consequence, and for the sake of simplicity, we established a hypothesis that assumes the dust density peaks at a given angular position in the ring and that the region of high-density has a given azimuthal width. We note that, in the simulations of Jackson et al. (2014), the dust density conveniently follows a Lorentzian function (see Fig.9a of Jackson et al. 2014). We thus implemented an azimuthal density gradient in the GRaTeR code such that the typical dust density function is multiplied by Article number, page 5 of 21 A&A proofs: manuscript no. main  Perrot et al. (2016). The same ellipses when overlapped in the IRDAP-reduced polarized intensity (b) indicate that the position of these structures were not well constrained in the total intensity data. The labels and arrows indicating R1, R2 and R3 in both images are exactly at the same position. The image (b) also shows the orientation of the best-fit model of GPI at the PA of 5 • . Both images are at the same spatial scale and represent the measured contrast. a Lorentzian expressed as: where d is the maximum density, ϕ 0 is the position angle of the collision site, and w is the width of the Lorentzian function. Throughout the paper, we fixed the position angle of the disk PA to 356 • and the disk aspect ratio to 0.02. While the later is chosen as a guess, Olofsson et al. (2020) showed that it can affect the measurement of phase function depending on the inclination of the disk. Therefore, to model the ring R3 we considered the following seven free parameters, which control the quantities R, Z and A: 1. The anisotropic scattering factor, g in the HG function. Forward scattering by the dust particles favours the value of g between 0 and 1 while the backward scattering renders g between 0 and -1. Here, g refers to g pol for polarimetric intensity and g total_int for total intensity, 2. The radius r 0 (au) of the disk where the dust density peaks, 3. The inclination i ( • ) of the disk, 4. The radial density profiles defined by two power law indexes α in and α out , 5. The azimuth of the peak in dust particle density, ϕ 0 ( • ), 6. The angular width of the Lorentzian function w ( • ).

Modeling of the polarized intensity
To justify the need of introducing an azimuthal density variation in the modeling, we first performed a test by choosing g pol as the only free parameter and considering that the density is constant along the ring. Here, our objective is to show that whatever phase function we assume, it does not emulate the north-south asymmetry.
We simulated a series of 50 GRaTeR models from g pol =0.1 to g pol =0.99, and set the other four parameters to: r 0 = 45 au, i = 57 • , α in = 20, and α out = −20, as defined in Perrot et al. (2016). To make a comparison with the DPI image, the models were convolved with the PSF and scaled globally in intensity to the data, which were binned to 2 × 2-pixels. The best-fit model is determined with the reduced χ 2 minimization, which is measured within an aperture encompassing the ring R3 (Fig. 4a). We found the best g pol value to be 0.54 (reduced χ 2 =1.6). However, as shown in Fig. D.1a, the best-fit model with this method failed to reproduce the north-south asymmetry observed in the data (Fig. 4a), and produced excess intensity to the north as shown in the residual (Fig. D.1b).
Then, we implemented the modified density function in GRaTeR and created 144000 synthetic models. The values of the radius and the inclination remained fixed as before. The rest of the five free parameters (g pol , ϕ 0 , w, α in and α out ) were explored with a coarse grid-spacing. The second column of Table 2 shows the range of the inner and outer bounds of the free parameters, total number of values explored, and the best-fit values and error bars obtained. The data, the best-fit model from the gridsearch approach, and the residuals are shown in Fig. 4a, b, and c, respectively.
The best-fit model favors a high value of g pol , which is indicative of a strong forward scattering, and a density peaking at ϕ 0 = 219 ± 18 • , which is the southwest direction in the image. The width of the density function is rather large expanding on 145 ± 43 • (errors are computed according to Bhowmik et al. (2019)). Counter-intuitively, this first analysis suggests that the peak of density is on the southwest part of the disk, opposite to the observed intensity peak with respect to the major axis. Therefore, the large g pol value (very forward scattering dust) and an azimuthally dependent density profile are necessary to reproduce the intensity peak positioned on the southeast as observed in the polarimetric image.
The best values for the two parameters describing the slope of the radial density function, α in and α out , correspond to a sharp ring inward and a shallower profile outward. As a comparison, Perrot et al. (2016) found that the ring was very steep both inwards and outward when observed in total intensity. However, Perrot et al. (2016) did not systematically explore several values of α in and α out , but instead tested only three configurations with α in = −α out = [5, 10, 20]. In addition, the differences with Perrot et al. (2016) can also arise from the self-subtraction as long as the ADI method acts as a high-pass filter (Milli et al. 2012), with a tendency to remove the faint signal produced by the diffuse dust population extending outside of the main ring. In this particular case, DPI allows us to better constrain the actual radial width of this dust ring.
Given the outcome of the grid search in terms of the presence of a density enhancement at the opposite of the polarized intensity peak, we embarked on a more detailed exploration of the parameter space with MCMC to check the robustness of the solution. In this case, we considered seven free parameters with 24 walkers. We chose priors (shown in Table 2) that were exactly the same as the range of parameters in the grid search. The initial guess values of the free parameters were chosen to be close to the best-fit parameters provided by grid-search. Even though the length of the burn-in phase was 500 iterations, we kept running the MCMC with a total of 7738 iterations to roughly compute the equal number of models as generated in the grid search. Similar results were obtained with a mean acceptance fraction of 0.45, indicating a converging solution (Gelman & Rubin 1992). Figure 5 shows the posterior probability distributions of all the free parameters, which were computed using Python's corner module (Foreman-Mackey 2016). We derived uncertainties based on the 16 th (−1σ), 50 th (median value), and 84 th (+1σ) percentiles of the samples in the distributions, which are plotted in vertical lines in Figure 5. The median and the error bars for each parameter are shown in the third column of Table 2, which are found to be consistent with the grid-based approach, particularly with regard to the fact that the density is peaking at the southwest (ϕ 0 = 223.4 •+2.7 • −3.0 • ) and that the anisotropic scattering factor is very high (g pol = 0.92 +0.05 −0.06 ). Figures 4d and 4e show the model and the corresponding residual, which also look similar to the solution suggested by grid-search (Figs. 4b, and c). To depict the dust distribution in the ring inferred by the model, we provide a deprojected view of the ring density variation in Fig. 4f, which is convolved with the PSF.
We compared these geometric parameters with those obtained with GPI observations in Bruzzone et al. (2020). They found an inclination of 60 • ± 10 • , which is compatible with our findings. Radial extension is also consistent with a radius r 0 = 44 +8 −12 au and a low density slope outside the main ring α out to be −1.0 +0.5 −1.0 . This is much shallower than the value found in our analysis (α out = −4). However, the outer slope relies on the sensitivity of a specific observation and therefore depends on the observing conditions. Moreover, we also restrained the fit to the location of R3 to avoid mixing R2 and R3, which inevitably would result in a shallower slope.
We also performed a sanity check by running another MCMC with different initial conditions of the parameter space. The goal was to evaluate the possibility of a solution around the first intuitive interpretation, where both the intensity and the density peak lie in the southeast of the polarized intensity image. We started an MCMC chain of 70 walkers by considering the initial value of free parameters g pol , r 0 , i, ϕ 0 , w, α in , and α out to be 0.3, 48 au, 58 • , 100 • , 18 • , 15, and -4, respectively. We definitely observed a local minima around this position. The model corresponding to this local minima (reduced χ 2 =1.1) and its residual are shown in Fig. E.1. This model does not visually show an obvious difference with the one previously obtained (shown in Fig. 4d). However, the MCMC eventually clearly disfavors this local minima and converges toward the former best solution, where the density peak is in the southwest (same MCMC results were obtained as shown in last column of Table 2).

Modeling of the total intensity image
In this section, we use a similar modeling analysis on the total intensity data reduced using a KLIP-ADI algorithm with KLIP truncated at five modes. Given the lower quality of the data, which is further heavily self-subtracted by the post-processing algorithms, we chose only three free parameters (g total_int , ϕ 0 , w) out of the total seven parameters mentioned in Sect. 4.2, with an objective to check whether the total intensity data were compatible with a non-uniform azimuthal density. We started our first analysis by fixing the values of the rest of the parameters based on Perrot et al. (2016), such that r 0 =45 au, i=57 • , α in =20 and α out =-20. We generated 206400 synthetic forward models following the ADI processing (e.g., Bhowmik et al. 2019), which took the self-subtraction bias into account similarly to the total intensity data (Fig. 1d). A grid of models with fine sampling was generated with GRaTeR, while the reduced χ 2 minimization provided the best-fit values in Table 3. Figures 6a, b, and c show the total intensity data, the best-fit model, and the corresponding residual. We also used an MCMC analysis by choosing the best-fit values from the grid-search as the initial guess while maintaining the same prior as in DPI modeling in Table 2. Both the grid search and MCMC simulations (burn-in phase of 100 iterations) reached a similar conclusion: the density peak ϕ 0 is roughly located at around 230 • , which is consistent with the DPI modeling results. On the contrary, the width of the density function, w ∼ 40 • is much smaller than in DPI. This could be  Table 3 shows the median results and 1σ uncertainty, which although marginally consistent with former tests, still confirms the peak of density in the southwest in both DPI and ADI data. We also double checked our outcome by running MCMC with less aggressive ADI parameters (three modes of the KLIP-PCA instead of five). The resulting analysis is of similar quality than that involving five modes (χ 2 = 3.1 instead of 3.7), but it is still marginally compatible to the values in the last column in Table 3 (g = 0.14 ± 0.02, φ 0 = 228.5 ± 3.9, w = 69.8 ± 6.0). As a result, the modeling of the total intensity image seems quite sensitive to the input data, although the location of the density peak is still very consistent within the error bars.
Because of the high impact of the self-subtraction on the ADI dataset, we can assume that the former MCMC results on total intensity data are more prone to biases than in DPI in what concern the geometrical parameters. For this reason, we performed one more MCMC run where the only free parameter is the g total_int value, which we can expect to vary between DPI and ADI data sets as opposed to the others. The six geometrical parameters were kept fixed according to the DPI modeling: After 1000 iterations (24 walkers), we obtained the best value of g total_int as 0.17 +0.02 −0.02 (reduced χ 2 =3.9). This value is consistent at 3 σ with the value (0.23 +0.02 −0.02 ) obtained in the previous test, although still in the lower range of g total_int . Given the poor quality of the data (the northern side of the disk is not detected), we did not develop further this analysis.

Polarimetric phase function
In the model presented in Sect. 4, we assume that the azimuthal intensity distribution in the polarimetric image depends on the product of two terms, the polarized scattering phase function (being itself a product of the Henyey-Greenstein function and Rayleigh scattering), and the density function, as a consequence of the disk being optically thin. Given that in the previous sections we limit ourselves to analyzing the grain distribution of the disk, here we aim to interpret relevant information about the grain properties from the extracted polarized scattering phase function (pSPF) and a spectrum in the NIR wavelengths.
In this section, we provide a measurement of the pSPF assuming the density term is determined from the model fitting performed in the previous section. For that we need to demodulate the density term from the surface brightness (SB). We first extracted the SB of the ring R3 with an azimuthal sampling of 1 • by fitting a Gaussian function radially from the star within the same aperture as shown in Fig. 3. The SB as a function of azimuthal angle is shown in Fig. 8. The error bars are the standard deviations measured in the U φ image within the same aperture at every azimuthal angle. The SB of the best-fit model is also over-plotted to compare with the observed data. The SB of the ring peaks at 172 • and there is a steep decrease beyond this PA in accordance with the absence of flux in the western part of the ring. We notice in the same figure that the model does not exactly match the data and it deviates significantly from an azimuthal angle of ≈ 40 • to 140 • . The largest deviation of 71% is reached at an azimuthal angle of 109 • .
To estimate the contribution of the density term to the surface brightness, we built a synthetic total intensity image in which the photometric effect of the scattering is artificially removed, with an anisotropic scattering parameter g pol =0 (isotropic scattering). The other model parameters were kept the same as the ones found with the MCMC as given in Table 2. This GRaTeR model is convolved to the observational PSF, and scaled to the Q φ image. Hereafter, we refer to the density model as long as it is based solely on the density function. In the optically thin regime, the SB of the density model provides the dust density distribution of the ring. Therefore, in this framework, the pSPF is obtained by taking the ratio of the SB of the ring by that of the density model. The scattering phase angle is calculated as in Milli et al. (2017, Eq. 2), which here corresponds to values ranging from 33 • to 147 • in the ring R3 (which is 90 • ± i). The resulting pSPF is plotted in Fig. 9, separately for the northern and southern sides of the ring, and normalized to its maximum value which is reached at a scattering angle of 43.8 • . Table 3: Free parameters constrained with the grid-search and MCMC methods while modeling the total intensity data. The second column assumes fixed parameters from (Perrot et al. 2016): r 0 =45 au, i=57 • , α in =20, α out =-20, and the third column assumes the values obtained for the MCMC best-fit DPI model (Table 2) Fixed parameters from Perrot et al. (2016) Fixed parameters from DPI model While the process of measuring the pSPF from the data is dependent on the model due to the hypothesis regarding the density function, it should be noted that this is in disagreement with the numerical relation based on the Henyey-Greenstein and Rayleigh functions for g pol = 0.92 (χ 2 ∼ 30), as determined by the best-fit model. Therefore, we tried to estimate the value of g pol that generates an analytical SPF conforming to our measured pSPF. We found that a better matched pSPF can be obtained for g pol = 0.65±0.08, although it is not fully satisfactory with regard to χ 2 ∼ 12. The model and its corresponding residual is shown in Fig. F.1. This value happens to be in agreement, within 3σ error bars, with g pol = 0.90 ± 0.09, as obtained by the grid-search (Table 2) method.
Several debris disks have been modeled with two anisotropic scattering factors (Milli et al. 2017;Bhowmik et al. 2019), essentially one for the forward (g 1 ) and the other for the backward (g 2 ) direction. We also explored this scenario and found a higher reduced χ 2 value (χ 2 ∼ 30) compared to that of the best-fit model with a single anisotropic scattering factor with g pol = 0.65 ± 0.08. Using the latter instead of the MCMC value from Sect. 4 (g pol = 0.95), we generated the surface brightness curve while keeping the same density function. This resulted in a much better match of the actual SB, particularly between 40 • and 140 • , as can be seen in Fig. 8 (pink dotted line).

Spectroscopy
We performed a spectral analysis of R3 using the IFS and IRDIS data obtained in March 2016, as these cover the largest spectral range (YJH-K1K2) available and the S/N is higher in the IFS than for the May 2015 IFS data. Due to the large background noise in the IRDIS K2 band image, we retained the K1 band data only. In this section, we extract a total intensity spectrum of the southern part of the disk (Fig. 6a) that we correct for the self-subtraction. We then constrain the grain properties by fitting the spectrum with synthetic models.
To account for the self-subtraction biases in the ADI reduced data, we used the best-fit GRaTeR model obtained from the MCMC analysis (Table 3), following the procedure detailed in Bhowmik et al. (2019). The model is convolved with a cropped PSF (0.4 radius for IRDIS and 0.3 radius for IFS) such that 99.99% of PSF's total flux is retained. In our analysis, we considered only the southern part of the disk inside the elliptical contour shown in Fig. 3 and normalized the best-fit KLIP model to the data.
We measured the average azimuthal intensity of R3 at each wavelength on deprojected images assuming i = 56.3 • and PA = 356 • , then converted it into contrast per arcsec −2 by normalizing it to the PSF flux, along with using the respective pixel scales (12.25±0.01 mas/pix for IRDIS and 7.46±0.02 mas/pix for IFS Maire et al. 2016). The outcome is the average spectral reflectance as a function of wavelength displayed in Fig. 10.
Globally, the spectrum has a mild negative slope (blue color) of (−7.97 ± 0.53) × 10 −5 arcsec −2 µm −1 across the Y JHK1 spectral range. The error bars shown in Fig. 10 are the dispersion in the measurements of the de-projected residual images (similar to Fig. 6c) at each wavelength within the same southernelliptical aperture used to extract the photometric measurements. The larger error bars at 1.4 ± 0.05 µm are because of the telluric water absorption feature as a result of which the disk is very faint. At around 1.5 µm, the spectral reflectance features a broad absorption pattern which might be attributed to OH bonding resonance in materials. However, the reliability of this feature is questionable given the error bars, calling for better S/N images in the future.

Modeling the grain properties
The SPFs of debris disks such as HR 4796 Olofsson et al. 2020) and HD 15115 ) were compared with theoretical models of SPF generated using Mie theory in order to put constraints on grain properties. Following this approach, we considered models simulated using the radiative transfer code MCFOST (Pinte et al. 2006). The grains in these models are porous in nature and composed of a mixture of astro-silicates (Draine & Lee 1984), carbonaceous grains (Li & Greenberg 1997), and water-ices (Li & Greenberg 1998) that can partially occupy the porous fraction of the mixture. We created 780 synthetic pSPF models using the following parameters: the minimum dust grain size, s min , varying between 0.1-100 µm, the fraction of vacuum removed by the ice, p H 2 0 , ranging between 1% and 29%, the porosity without ice, P, between 0% -80% (Augereau et al. 1999b), and the silicates to carbonaceous grains volume fraction, q sior , which can take values of 1 (silicates only), 2 (same amount of silicates and carbonaceous grains), or ∞ (carbonaceous only). Beside, we fixed the slope of the power size distribution κ to -3.5 as predicted by Article number, page 9 of 21  Fig. 4: DPI data is displayed as a reference in the top row, together with the best-fit models (left) and the corresponding residuals (right) resulting from the grid-search (second row), and the MCMC (third row) analysis methods. Table 2 provides the parameters of the model for these two cases. The ellipses in dashed white lines show the aperture within which the reduced χ 2 minimization has been performed. The bottom row shows the deprojected density function corresponding to the MCMC best-fit model. Dohnanyi (1969) and maximum grain size s max = 1 mm. The total porosity of the grains, p wH 2 0 , is a derived quantity defined as p wH 2 O = P (1 − p H 2 O /100) (Augereau et al. 1999b). All the models are normalized to the maximum value of the pSPF before performing the reduced χ 2 analysis. Similarly, we used the 780 MCFOST models to fit the reflectance spectrum, by integrating over the range of accessible scattering angles (33 • to 147 • ), sampled and normalized on the 40 spectral channels corresponding to the IFS and IRDIS data, followed by the reduced χ 2 analysis.
Motivated by the goal to derive global grain properties, we aimed to find a single model fitting both the pSPF and the reflectance spectrum, based on the average χ 2 . The best fit displayed in Fig. 9 and Fig. 10 (violet line, reduced average χ 2 = 7.34) corresponds to s min = 0.1 µm, p H 2 0 = 1%, P = 40%, and q sior = 1 ( Table 4). The probability distribution (Fig. 11) shows that the minimum grain size is ranging between 0.1 -0.2 µm. This is apparently inconsistent with the larger value of 4 µm derived by Bruzzone et al. (2020). However, we recall here that Bruzzone et al. (2020) directly fits MCFOST synthetic images considering an azimuthally uniform ring and limited the composition to a spherical astro-silicates type of grains, so these results cannot be directly compared. In addition, the analysis favors P = 0.4 and p H 2 O = 0.01, which translates to a total porosity of 40%, so again it is discrepant with respect to a value of 0% by Bruzzone et al. (2020). Finally, while the maximum probability is reached for q sior = 1, there is a 25% probability that the ring has twice the amount of astro-silicates than carbonaceous grains.
Also, as seen in Fig. 10, the best model fails to fit the data point corresponding to the K1 filter at 2.11µm. If we ignore the K1 data point, the average χ 2 values tend to reduce noticeably to about 4. Making the spectrum bluer would have required a steeper grain size distribution which is unusual, hence, it is not considered here. Moreover, fitting the pSPF and the reflectance spectrum separately provides a similar solution for the grain properties for the former, but a much larger value of s min = 3.2 µm for the latter. Therefore, the reflectance spectrum seems less relevant to derive the grain properties, in this particular case. Overall, the exercise is showing the limits in deriving grain properties with both a limited range in scattering phase angles and in wavelengths, together with a strong assumption on the grains sphericity. A more detailed analysis taking into account all the available observables is beyond the scope of this work but this will be key in setting meaningful constraints on the grain properties in future studies.

Grain scattering and composition
The estimation of the scattering asymmetry factor g pol shows a strong dispersion. data yields g pol = 0.92 +0.05 −0.06 or g pol = 0.90 ± 0.09. Even though the residuals of the fit are acceptable, these are extremely high values as far as grain properties are concerned. This would imply very large grains with potentially negative signal in Q φ and non-null signal in U φ , at some scattering phase angles, which is incompatible with the data. We note that the error bars on the g pol parameter are particularly small and potentially unrealistic with the MCMC approach, although the mean values are consistent.
When relying on the measured SB of the same polarimetric data, we find a relatively smaller estimate for the g pol parameter (g pol = 0.65), which appears more consistent with the other constraints and typical values already measured in such systems Milli et al. 2019). The SB is indeed inconsistent with the g pol = 0.92 model at scattering phase angles between 40 and 140 • . However, in this region of the ring, the SB is essentially determined by the azimuthal density function, as can be seen in Fig. 8. In terms of image residuals, a model with g pol = 0.65 provides a larger reduced χ 2 (0.90 instead of 0.53) but still acceptable visually. We also note that the proximity of the rings R2 and R3 in the eastern part can cause biases in the image fitting procedure, while the SB extraction should be less sensitive to that. From this analysis we conclude that degeneracies between the density and the scattering parameters prevent us A&A proofs: manuscript no. main from deriving an accurate value of g pol , thus it is more reliable to adopt a lower limit g pol > 0.65.
MCMC model fitting on total intensity images provides an even lower value of g total_int = 0.14 − 0.23. While the scattering asymmetry factor can differ from polarized to total intensity data as a result of grain properties, and as already witnessed in other debris disks, the present case appears rather extreme. However, even if we generated forward models, we cannot exclude that the self-subtraction inherent to ADI, by suppressing the disk flux along the minor axis (which is oriented along the forwardbackward direction, where the observation should be the most sensitive to scattering effects) would lead to an underestimation of g total_int , as it was the case for HR 4796 (Milli et al. 2017). Reconciling the derived values of the anisotropic scattering factor would certainly require higher S/N polarimetric images and less biased total intensity images, which is beyond the scope of this paper.

Considering planetesimal disruption as a possible origin for the ring
Large planetesimal disruptions in the recent past are a natural way of producing density asymmetries along the ring. Such disruptions may happen through direct collisions between planetesimals (Jackson et al. 2014;Lawler et al. 2015), or, if a massive planet is present in the vicinity, planetesimals can also be tidally disrupted if they pass within the tidal disruption radius of the planet (Cataldi et al. 2018). Collisions tend to dominate disruption rates in cases where the parent planetesimals are small, while tidal disruption dominates in cases where the planetesimals are large (Janson et al. 2020). Jackson et al. (2014) have shown that the collisional breakup of large planetesimals "naturally" produces strongly asymmetric rings, which are made of the small particles produced by the collisional cascade created in the initial breakup's aftermath. This disk is asymmetric because all fragments have eccentric orbits imparted by the velocity kick they get in the aftermath of the breakup, with a pericenter that is necessarily located at the location of the breakup. As a result, this creates a ring that is narrow and bright on the side of the breakup and more extended, diffuse, and faint on the opposite side. The survival time of these asymmetries is set by the progressive spread due to collisions among the fragments (Kral et al. 2015) or by orbital precession due to other planetary bod- Fig. 9: Polarized scattering phase function of the northern (pink) and southern (blue) part of the ring as a function of scattering angles. The solid black line shows the analytical SPF with g = 0.92 and the shaded area represents the analytical best-fit SPF with g = 0.65±0.08. The best fit to the spectrum with Mie grains is shown as a violet line. All plots are normalized to the maximum value of pSPF at the scattering angle of 43.8 • . Fig. 10: Reflectance spectrum of the ring R3. The violet line represent the spectra corresponding to MCFOST's best Mie model that is normalized to the reflectance spectrum.
ies (Jackson et al. 2014) and can typically be as large as a few million years at 50 au. In this collisional-origin scenario, one crucial parameter is the initial velocity kick imparted to the collisional fragments, which constrains the eccentricities of the fragment orbits and thus determines the ring's extent and shape. This velocity kick can be parameterized as a departure, σ k , with respect to the local circular Keplerian velocity ( k ). To a first order, σ k is set by the mass M 0 of the shattered large body and the location of the breakup. For example, bodies with masses of Ceres and Pluto colliding at 50 au around a Sun-like star would result in σ k / k = 0.05 and 0.1, respectively. The azimuthal extension of the high-density region around the breakup location is set by the value of σ k / k : the higher this ratio, the smaller the width, w, of the azimuthal Lorentzian distribution becomes and the more asymmetric the disk (see Fig. 9a of Jackson et al. (2014)). While Jackson et al. (2014) give no explicit analytical relation between w and M 0 , we can use their Fig.9a to derive an approximate relation between w and σ k /v k , for which we get: We can now use the fact that σ k /v k ∝ M 1/3 0 to estimate the mass M 0 of the object, whose breakup produces a Lorentzian density distribution having a width, w, that is compatible with the one obtained in our best-fit. We remain careful and consider the value obtained in our intensity data-fit, w = 57 o , and the value obtained in our polarimetric fit, w = 119 o , as bracketing values for w. In this case, we get from Eq.4 that 0.015 ≤ (σ v /v k ) ≤ 0.03. For a typical asteroidal density of ∼ 3g/cm 3 , this corresponds to an object of mass M 0 with 7.5 × 10 −6 M ⊕ ≤ M 0 ≤ 6 × 10 −5 M ⊕ and a radius comprised between ∼ 150 and ∼300 km.
This size domain is well below the smallest object (of Ceres size) that was explored in the simulations of Jackson et al. (2014). As the brightness of the breakup-produced ring is directly linked to the mass of the initially shattered object, this raises the question of whether or not a 150-to-300 km object can account for the brightness of the inner ring, which is around 10 −2 that of the star. This value is obtained by integrating the flux of the ring within the aperture as per the procedure explained in Sect. 6. If we consider the most optimistic case (in terms of the total geometrical cross-section), where all the mass of the initial object is entirely transformed into a cloud of micron-sized grains, we obtain that this dusty cloud intercepts a fraction of stellar light (supposing a perfect albedo=1 case) comprised between 2 × 10 −6 (for a 150 km parent body) and 1.5 × 10 −5 (for a 300 km one). This is at least three orders of magnitude less than the observed luminosity contrast. This observed luminosity contrast of ∼ 0.01, which should be of the order of the ring's optical depth, would correspond to a lunar-sized progenitor of at least radius ∼ 1500 km.
The collisional breakup scenario, as proposed in Jackson et al. (2014), thus seems highly unlikely, essentially because there is an inconsistency between the high luminosity of the disk, which would be the signature of a very large collisional progenitor, and the wide azimuthal extension of the density profile, which would on the contrary be the signature of a smaller progenitor.
We note that the alternative breakup scenario by tidal disruption in the vicinity of a giant planet encounters the same problem as the collisional breakup hypothesis: the velocity dispersion of the produced fragments has to be at least comparable to the escape velocity of the shattered object. This means that a large object (which is needed to generate the bright disc of debris that is observed) would also necessarily produce debris on high-e orbits, leading to a ring that is too peaked in azimuth when compared to observations. A way to circumvent this problem would be to assume that the ring is the result of a succession of more than a thousand collisional breakups (or tidal disruptions) of ∼ 150 − 300 km progenitors happening in close succession at exactly the same spatial location. While this scenario cannot be ruled out beforehand, it seems, however, highly contingent and non-generic. The presence of gas in HD 141569 has always been considered an important factor for explaining its particular morphology. The most recent measurements obtained with ALMA indicate large amounts of CO gas distributed at the same distances as the ringlets detected with SPHERE (Di Folco et al. 2020). Such a large quantity of gas could have a dynamical effect on the smallest grains in the system. The fact that the CO gas distribution is almost azimuthally homogeneous means that gas should naturally induce a more axisymmetric dust distribution than what is predicted in Jackson et al. (2014) in a gas-free case after the breakup of a massive object. Therefore, the azimuthally peaked distribution of dust observed with SPHERE could have been produced by a more massive object than what is estimated in Sect. 8.2, and gas would broaden the extent of the smallest dust grains azimuthally. This scenario could help solving the main inconsistency in the disruption scenario in which the progenitor is too small to account for the actual luminosity of the ring. Here, we assess the validity of this hypothesis, based on ordersof-magnitude estimates of the gas mass in the disk and of its dynamical effect on the dust.
Based on ALMA data, Di Folco et al. (2020) estimated the mass of CO gas, respectively of dust, to be M CO ∼ 0.1M ⊕ , and M dust ∼ 0.03 − 0.52 M ⊕ , leading to a CO-to-dust mass ratio in the range 0.19-3.3. Under the assumption that the gas is of secondary origin (and not primordial as posited in Di Folco et al. 2020), the hydrogen content can be much lower than in a protoplanetary disk ), but there may be a substantial amount of carbon and oxygen produced by the photodissociation of CO (Kral et al. 2019;Marino et al. 2020). The CO mass in the system is high enough, so that it is self-shielding (Visser et al. 2009) and carbon may also start shielding CO from photodissociating. Then, CO can viscously spread, which would explain why it is so broad (extending over 200 au) in a secondary origin context (Kral et al. 2019). The total carbon mass then may be lower or greater than the CO mass depending on unconstrained parameters, such as the viscosity α of the gas disk. However, assuming a high α value considering that the magnetorotational instability may be very active in these disks (see Kral & Latter 2016 and the first observational confirmations in Marino et al. 2020), the carbon and oxygen masses should supersede the CO mass (Kral et al. 2019;Marino et al. 2020), leading to a gasto-dust ratio larger than 1. Therefore, the gas-to-dust ratio lower limit is 0.19, but it is likely to be greater than 1.
Given this rather high gas-to-dust ratio for a potentially secondary gas disk (Kral et al. 2017;Matrà et al. 2017), we go on to explore the effect of gas on bound and unbound dust grains in this system. The gas surface density in CO at 44 au in HD 141569 is Σ CO ∼ 1.5 × 10 −4 g/cm 2 (Di Folco et al. 2020), or 8 × 10 −5 g/cm 2 using the more optically thin 13 CO line (and assuming a typical isotopic ratio of 77 between 12 CO and 13 CO) giving a lower limit to the total gas surface density. Now, we can compute the time for a bound grain to couple to gas. Using Richert et al. (2018) and Burns et al. (1979), the dimensionless stopping time (in orbital units) reads as: where s, ρ and Σ g are the grain size, grain density, and gas surface density, respectively 3 . Therefore, grains that are 50 microns should see their orbits strongly affected by gas in fewer than 100 orbits (Takeuchi & Artymowicz 2001), which is shorter than the typical collisional time for a disk of optical depth ∼ 0.01. This means that any azimuthal clustering of dust grains 3 We note that if primordial the total gas mass and surface density could be several orders of magnitude higher than for a secondary gas production.
should broaden on a timescale of only a dozen orbits for grains close to the blow-out size (which is quite high in this system with s blow ∼ 8 microns 4 ). Given the SPHERE images are representative of light scattered by dust, we also need to investigate a similar quantity for unbound grains which can be smaller than the blow-out size. These potentially unbound micron to sub-micron grains would be affected by gas drag and following the same procedure as in Bhowmik et al. (2019) (assuming T=35 K from Di Folco et al. 2020 and M = 2 M from Merín et al. 2004), we find: where ∆R is the width over which a grain feels the gas. This stopping time for unbound grains can be even shorter than T s . Grains below the blow-out size could strongly feel the gas and may circularize rapidly (and become bound), thus accumulating instead of leaving the system. This conclusion is similar to what was found for HD 32297 in Bhowmik et al. (2019).
These results show that small grains produced in a disruptive-like scenario, such as the one described in the previous section, may quickly move onto different orbits, which would broaden the expected density peak, both azimuthally and radially. Therefore, we can conclude that the disrupted mass we found in the previous section is, in fact, a lower limit.
When the gas mass is taken into account, the disruptive scenario cannot be ruled out any longer because a larger disrupted body could indeed both explain the disk brightness and its azimuthal extent. It also means that fewer disrupting events may actually be needed to explain the current disk brightness if several disrupted bodies are to be involved. However, dedicated simulations with, for example, an improved version of LIDT-DD (Kral et al. 2013(Kral et al. , 2015Thebault & Kral 2018) that includes gas drag would be needed to explore this avenue further and give a stronger conclusion.

Potential of photoelectric instability to explain the observations
The photoelectric instability (Lyra & Kuchner 2013) has been identified as a possible mechanism that could potentially alter disk structures and create asymmetric arc-like features, such as those observed in this study (Richert et al. 2018). The amount of CO gas derived by Di Folco et al. (2020) corresponds to the regime of "high-level of gas" described in Richert et al. (2018), for which the instability is the strongest and shows no anticorrelation between gas and dust. Interestingly, Di Folco et al.
(2020) demonstrated that the CO is not entirely axisymmetric, but is instead enhanced between 220 • and 290 • , thus peaking in the west-southwest part of the ring; hence, this finding is in qualitative agreement with the dust distribution inferred from scattered light observations. Both gas and dust seem slightly correlated.
The photoelectric linear instability may be stabilized for gasto-dust ratio < 1. As explained previously, we expect the gas-todust ratio to exceed 1 based on only CO and its daughter species carbon and oxygen. This instability may then explain the presence of the three rings (located within 100 au) observed with SPHERE as the gas observed with ALMA extends up to more than 200 au (White et al. 2016;Miley et al. 2018;Di Folco et al. 2020). This is one tempting possibility, although it goes beyond the scope of this paper to explore it further by running dedicated simulations of this instability for HD 141569.

Explaining the slight asymmetry observed in the gas
Finally, a relevant question considers whether the asymmetric structure of the gas revealed with ALMA at separations of ∼50 au (see previous subsection) can produce a pressure trap that could be responsible for the dust distribution. First, we note that the radial location of the gas asymmetry (Di Folco et al. 2020) seems to be at a slightly larger value than the ring studied in this paper, but there is some overlap and a higher resolution would be needed to pinpoint the gas peak location. Also, the excess of gas emission (CO clump) is 6 σ above the noise in the data by Di Folco et al. (2020), but given the line is optically thick, it does not give any information on the surface density compared to the rest of the gas disk, but we do note that an increased temperature at the ring location would indeed create more emission in CO. We cannot rule out a pressure trap (a more optically thin line would be needed) but we can suggest another possibility.
Due to the high quantity of dust in the ringlets and due to the photoelectric effect, there could be a lot of electrons at these specific locations that would heat the CO gas and increase its brightness, hence creating the illusion of a CO clump in the ALMA data. In this case, the significant dust asymmetry observed could thus be responsible for the CO clump.

Conclusions
We obtained SPHERE polarimetric data of the gas rich debris disk around HD 141569, as well as new total intensity data complementing the study presented in Perrot et al. (2016). We definitively confirmed the presence of several concentric rings around the star at few tens of au.
In this paper, we focus on constraining the dust distribution in the innermost detected ring located at ∼45 au, based on the combined information from polarized and total intensity images. We posited that the dust density must be non-uniform azimuthally to account for the observed intensity distributions. Assuming that the ring azimuthal density can be described with a Lorentzian function, we found that the dust is peaking in the southwest, while the intensity peak is located in the southeast, as a result of the light being preferentially scattered forward. This main finding is robust whatever the approach taken to model the polarimetric and total intensity images.
The values derived from total intensity and polarimetry are consistent providing overall a position angle of 229.1 • ± 9.0 • , as determined with MCMC. Interestingly, this location is consistent with the enhancement of CO measured with ALMA by Di Folco et al. (2020), possibly pointing to a link between the distributions of gas and dust. On the contrary, the azimuthal width of the density function is rather broad but differs by a factor of 2 between polarimetric (118.9 •+8.2 • −7.3 • ) and total intensity data (57.9 •+3.2 • −2.9 • ). In addition, the polarimetric image fitting procedure converges to a very high value of the anisotropic scattering factor (g pol = 0.92), while we obtained a moderate value for total intensity (g total_int = 0.14 − 0.23).
Based on the modeling of the images we extracted the polarized scattering phase function (pSPF) by demodulating the azimuthal density variation from the surface brightness. We concluded that a lower value of the anisotropic scattering factor (g = 0.65) provides a better match to the pSPF as compared to the image fitting. However, we note that the degeneracies between scattering and density can stand as an issue impeding a thorough measurement of this parameter. We also derived the spectral reflectance of the ring in the wavelength range of 0.98 − 2.1 µm. Modeling the pSPF and spectral reflectance with MCFOST and assuming Mie scattering provides a composition of porous sub-micron sized mixture of carbonaceous and astrosilicate grains. We argue that broader wavelength coverage and higher S/N data would be needed to derive more meaningful conclusions about the grain properties.
Finally, we discuss the implications of the dust distribution in terms of the dynamical history of the system. While massive collisions of planetesimals can produce a ring-like pattern with dust enhancement at the place of the collision, this scenario fails to account for the azimuthal dust distribution in the ring, since it suggests a small size for the progenitor, incompatible with the total luminosity of the ring. However, the gas mass is high enough (> 0.1 M ⊕ ) that gas can broaden any azimuthal clustering on a few dynamical timescales as grains move onto different orbits very quickly before being destroyed through collision. In this paper, we find that the disruption scenario can no longer be ruled out when gas dynamics is taken into account because a larger disrupting body can, in this case, account for both the high brightness of the disk and its azimuthal extent. We note that dedicated simulations including gas drag are required to further explore this scenario and reinforce our conclusion.