© ESO 2019

1 Introduction

Debris disks are the leftovers of stellar, and planetary, formation processes. As the primordial gas-rich, massive circumstellar disk evolves, the small dust grains grow and form planetesimals, which may become the building blocks of future planets. With a half-life time of a few million years (Hernández et al. 2007), this primordial disk rapidly transitions toward its debris disk phase, losing the vast majority of its gaseous content. In this debris disk phase, the planetesimals leftover from the planet formation process collisionally erode to produce the small dust grains that are observed.

Recent advances in high contrast imaging instruments, such as the Gemini Planet Imager (GPI, Perrin et al. 2015) or the Very Large Telescope/Spectro-Polarimetric High-contrast Exoplanet REsearch (VLT/SPHERE, Beuzit et al. 2019) provide us with new avenues to investigate debris disks at high angular resolution. For instance, Olofsson et al. (2016) and Milli et al. (2017) studied the scattering phase function over a wide range of scattering angles for HD 61005 and HR 4796 A, respectively, using the SPHERE instrument. Such studies help better constrain the nature of the small dust grains in young debris disks. Lee & Chiang (2016) present numerical simulations that can explain a variety of morphologies of debris disks. Their code, which can account for the effect of stellar radiation pressure, was used by Esposito et al. (2016) to model GPI observations of the debris disk around HD 61005. Overall, thanks to the exquisite spatial resolution provided by this new generation of instruments, we are able to perform in-depth studies of young and bright debris disks, trying to constrain the collisional activity responsible for the production of small dust grains.

Since Jura (1991) reported the detection of mid- and far-infrared (IR) excess emission with the InfraRed Astronomy Satellite, the debris disk around HR 4796 A has been intensely studied from the ground and from space with ever increasing image quality. HR 4796 A is young (8 ± 2 Myr, Stauffer et al. 1995) and nearby (71.9 ± 0.7 pc, Gaia Collaboration 2018). The A-type star is thought to be in a visual binary system, and hosts one of the brightest debris disks (fractional luminosity of about 5 × 10⁻³, Moór et al. 2006). The secondary is an M-type star at a separation of 7.7′′ (Jura et al. 1993). Scattered light and thermal emission observations have revealed a narrow ring of dust at about 77 au from the central A-type star (e.g., Jayawardhana et al. 1998; Wyatt et al. 1999; Augereau et al. 1999; Wahhaj et al. 2005; Schneider et al. 2009, 2018; Thalmann et al. 2011; Moerchen et al. 2011; Lagrange et al. 2012; Rodigas et al. 2015; Milli et al. 2015, 2017, 2019; Perrin et al. 2015; Kennedy et al. 2018). Several studies have reported that the disk displays a brightness asymmetry along its major axis, which is most likely related to the non-zero eccentricity of the disk. It has been postulated that planets could shepherd the parent planetesimals, inducing the observed azimuthal asymmetries via secular interactions (e.g., Wyatt et al. 1999). In this paper, we present SPHERE/ZIMPOL observations of HR 4796 A and we investigate the collisional activity and the production of small dust grains by constraining the azimuthal and radial distribution of the dust.

2 Observations and data reduction

The data were obtained with SPHERE/ZIMPOL (Schmid et al. 2018) in its polarimetric mode P2, corresponding to field stabilized observations. They are part of the SPHERE Guaranteed Time Observations¹. The target HR 4796 A was observed on the night of 24 May 2016 during a sequence alternating between two short unsaturated polarimetric cycles in fast polarimetry (called fast) and two deeper saturated polarimetric cycles in slow polarimetry (called slow). We repeated three times the pattern fast slow fast with the derotator position angle set to 0°, 30°, and 60° respectively. Two more cycles were actually observed, with position angles of 120° and 150°, but they are not used in this paper. For the last two cycles, columns of deeply saturated pixels fall right along the semi-major axis of the disk, and excluding them leads to an overall cleaner final reduction along the semi-major axis, but noisier in the regions close to the semi-minor axis. Given that the purpose of this paper is to study the radial profiles along the semi-major axes of the disk, we opted to exclude those additional cycles. This strategy enables us to get unsaturated frames that bracket the deep saturated exposures in order to calibrate the photometry. Furthermore, the different derotator position angles introduce an additional diversity parameter to smooth out any residual pattern. In this paper, we only focus on the polarized image of the disk without any absolute flux calibration. Therefore, only the deep unsaturated slow polarimetry images are used for the image discussed later. The absolute polarized flux of the disk and polarized fraction are treated in a separate paper (Milli et al. 2019). The extreme adaptive optics (SAXO; Fusco et al. 2006) yielded a point-spread function (PSF) with a full width at half maximum of 30– 40 mas, larger than the expected diffraction limit because of the low-wind effect (see also Milli et al. 2019).

The images were reduced with custom Python routines. No calibration frames (such as dark, bias or flat field) are applied to avoid introducing additional sources of noise that might degrade the polarized sensitivity. The ZIMPOL instrument is indeed designed to beat atmospheric and instrumental speckles by using the concept of polarimetric differential imaging with a masked CCD synchronized with a ferroelectric liquid crystal modulated at a rate of 27 Hz in Slow Polarimetry (Schmid et al. 2012). This unique design allows the signal coming from the two orthogonal polarization directions to be captured by the exact same pixels, self-calibrating any differential and pixel-to-pixel effect after applying the polarimetric subtraction. A half-wave plate (HWP) is introduced very early in the optical trainto calibrate the instrumental polarization. Each polarimetric cycle includes four positions of the HWP: 0°, 22.5°, 45° and 67.5°. We applied the double difference technique to obtain the two Stokes parameter Q and U out of these 4 HWP positions (Avenhaus et al. 2014). Additional instrumental polarization coming upstream from the HWP remains and are removed by subtracting from the Stokes Q and U a scaled version of the intensityimage I, as describedin Ginski et al. (2016) and Canovas et al. (2011). We then constructed local Stokes vectors Q_ϕ (shown in Fig. 1) and U_ϕ, containing the astrophysical signal and an estimate of the noise, respectively (Benisty et al. 2015; Olofsson et al. 2016, see also Canovas et al. 2015 for a discussion on the effect of optical depth).

3 Grain-size dependent dust distribution

Scattered light observations in the near-infrared are sensitive to the small-end of the grain size distribution (e.g., Mulders et al. 2013). Therefore, to properly characterize the dust production and the collisional activity, one needs an appropriate prescription for the behavior of the μm-sized dust grains. Besides gravity, the dominant effect that can affect the orbital parameters of small dust grains is stellar radiation pressure (Wyatt 2005 showed that Poynting-Robertson drag does not play a significant role in massive debris disks, but see Kennedy & Piette 2015 for a discussion about detecting dust in the inner regions using nulling interferometry). Consequently, in this paper, we try to constrain the radial and azimuthal distribution of the dust, taking into account the effect of radiation pressure.

3.1 Description of the code

The code used in this study is inspired by the work presented in Lee & Chiang (2016). We start with an analytical description of the belt of planetesimal-sized parent bodies that produce the observed dust. The ring is defined by 6 parameters: a reference radius r₀, eccentricity e, position angle on the sky ϕ (positive from north to east), inclination i, argument of periapsis ω, and width of the belt δ_r. The radial distribution of the parent bodies follows a normal distribution centered at r₀ with a standard deviation δ_r. All the dust grains, produced in a collisional cascade will originate from those parent bodies, whose sizes, numbers, or masses do not need to be explicitly defined in the code. However, the number of grains of different sizes released from those unseen parents bodies are defined by a size distribution between a minimum and a maximum size. Within this range of sizes, the grains are affected by radiation pressure from the star, and their dynamical evolution is size-dependent. One has to note that all parent bodies have a “forced eccentricity”, as they all share the same e and the same ω.

As input parameters, we consider a grain size distribution, which initially follows the “ideal” prescription of Dohnanyi (1969, a differential power-law of the form dn(s) ∝ s^−3.5ds, where s is the grain size). The distribution is divided in n_g intervals between the minimum and maximum grain sizes (s_min and s_max, respectively). The number of grains in each interval is computed as Eq. (2) of Dullemond & Dominik (2008) (1)

where Δlog(s) is the width of each bin in logarithmic space. Since the n_g grain sizes are logarithmically spaced, each Δlog(s) is the same, except the first and last ones which are half of that value (so that the grain size distribution is exactly between s_min and s_max). For each grain size, we then compute thedimensionless β ratio between radiation pressure and gravitational forces (Burns et al. 1979). For a given s, the value of β depends on the dust properties (optical constants and density) and the stellar properties (mass and luminosity), and is evaluated as (2)

where L_⋆ and M_⋆ are the stellar mass and luminosity, G the gravitational constant, ρ the dust density. The radiation pressure efficiency Q_pr(s) is equal to Q_ext(s, λ) − g_sca(s) × Q_sca(s, λ) averaged over the stellar spectrum, with Q_ext and Q_sca the extinction and scattering efficiencies (computed using the Mie theory), and g_sca the asymmetry parameter (the integral over 4π steradian of the phase function times the cosine of the scattering angle).

To decide where the collision releasing a dust grain takes place, we use a prior distribution on the mean anomaly. This “collisional distribution” can either be uniform or a normal distribution (centered either at the pericenter or the apocenter). The standard deviation when using the normal distribution is noted δ_ω. This implies that dust grains are not released uniformly in the disk, but that they can be released preferentially in a localized region of the disk, depending on the azimuth. The mean anomaly is then converted to the true anomaly ν, by solving the Kepler equation. The effect of radiation pressure on the orbital parameters of a single dust grain is parametrized as in Wyatt et al. (1999), Wyatt (2006), and Lee & Chiang (2016). Assuming that the parameters of the dust grain, upon its release, are a (drawn from the normal distribution of width δ_r centered at r₀), e, β, and ω, then its “updated” orbital parameters (a_n, e_n, and ω_n) are computed as (3)

For a given β, we make 3000 realizations of the prior collisional distribution. For each realization (giving a set of orbital parameters a_n, e_n, and ω_n), we check if the updated eccentricity e_n is larger or equal to zero and strictly smaller than unity to avoid hyperbolic orbits. Similarly to Löhne et al. (2017), the blow-out size does therefore depend on where the grains are launched from (e.g., pericenter or apocenter). If the orbit is bound, the code then populates it with 500 dust particles, uniformly distributed in mean anomaly. We finally then draw from a normal distribution of standard deviation h∕r = 0.04 to account for the vertical dispersion on the disk. The opening angle is set to 0.04 following Thébault (2009). The (x, y, z) positions of each particle are registered, and we then find the pixel of the image that is closest to the (x, y, z) values depending on the inclination and position angle of the disk (with the same pixel scale as the observations being modeled). We therefore produce number density maps for each value of β.

The modeling strategy described above does not take into account an important effect discussed in Strubbe & Chiang (2006) and Thébault & Wu (2008), which is that small grains produced inside the belt on high eccentricity (bound) orbits will spend most of their time in the collision-free outer regions where they cannot be collisionally destroyed. This will significantly enhance their collisional lifetimes, and thus their number density as compared to what would be obtained by simply spreading out the number density obtained with Eq. (1) over their whole orbit. To a first order, Strubbe & Chiang (2006) and Thébault & Wu (2008) found that a correcting “enhancement” factor should be applied to the high-β grain number density, which is roughly proportional to their total orbital period divided by the time spent within the birth ring. We take here the simplified expression for this correction factor given in Lee & Chiang (2016) and apply the following strategy. For each of 3000 particles at a given β that are not sent on hyperbolic orbits, we compute this correcting factor as , with α = 3∕2 (e being the eccentricity of the parent belt, Lee & Chiang 2016), and ν the true anomaly at the moment of the collision. This correction should naturally produce a surface brightness profile of r^−3.5. However, such a profile is an asymptotic behavior that is reached relatively far away from the parent planetesimal belt, but not right after the birth ring (Thebault et al. 2012). When computing the number density maps for each β values, the contribution of each particle is multiplied by this correction factor. One should also note that we do not take grain–grain collisions into account, which, as demonstrated in Löhne et al. (2017) can have an impact on the radial extent of the disk.

Once the 3000 × 500 particles have been launched, the scattering angle between the central star and the observer is computed for each pixel in the image. The code then computes one image per grain size bin by multiplying the number density of each pixel by S₁₂ × πs² × Q_sca∕(4πr²), where r is the distance to the star, Q_sca the scattering efficiency, and S₁₂ the polarized phase function (which can be computed using the Mie theory, or other means). By using the S₁₂ element, we are effectively computing the Q_ϕ image directly, andnot the Q and U images. The code can compute scattered light images replacing S₁₂ with the S₁₁ element of the Müller matrix. It can also compute thermal emission images by multiplying the number density map by 4πs² Q_absπB_ν(T_dust) (the dust temperature being evaluated by equating the energy received and emitted by a dust grain of size s at the distance r from the star, Eq. (5) of Wolf & Hillenbrand 2005). The final image is the collapse of all individual images, weighted by n(s).

3.2 Stellar parameters

To derive the stellar parameters, we first gathered optical and near-IR photometry using the VO SED Analyzer (VOSA) tool² (Bayo et al. 2008). The stellar effective temperature and surface gravity were retrieved using the VO-DUNES discovery tool³ (Eiroa et al. 2013), which explores VIZIER catalogs. We find a value of T_⋆ = 9700 K and a log (g) = 4.05. Fitting a Kurucz model (Castelli et al. 1997) to the optical and near-IR photometry, we find a stellar luminosity of L_⋆ = 25.75L_⊙. The distance of 71.9 pc combined with the dilution factor used to scale the; Kurucz model leads to a stellar radius of R_⋆ = 1.79R_⊙. The stellar mass is determined assuming that the star has finished contracting and we used the following relation between the stellar mass, radius and surface gravity, (4)

We obtain a mass of 1.31M_⊙.

3.3 Disk parameters

Milli et al. (2017) and Kennedy et al. (2018) presented comprehensive studies of the debris disk as seen in scattered light with VLT/SPHERE IRDIS and in thermal emission with ALMA, respectively. They constrain the main parameters of the disk, and because of our modeling strategy (see next sub-section), we used their results to set the inclination i to 76.6°. Kennedy et al. (2018) reported a value of 76.6° ± 0.2, Milli et al. (2017) a value of 76.45° ± 0.7, Schneider et al. (2018) a value of 75.9° ± 0.14, Schneider et al. (2009) a value of 75.88° ± 0.16, Thalmann et al. (2011) a value of 76.7° ± 0.5. Overall, our choice for the inclination is consistent with previous studies, especially given that it was derived from high signal-to-noise observations at high angular resolution. For the position angle ϕ, we used a maximum merit function as in Thalmann et al. (2011, see Olofsson et al. 2016 for details on how the elliptic mask is defined). We find that a value of ϕ = −151.6° provides a best fit to the observations, but this value will be re-evaluated during the modeling of the observations. Milli et al. (2017) found ϕ = −152.9°, but one should note that we use a different convention for the position angle as the one used in Milli et al. (2017). For ϕ = 0° we assume that the major axis is along the north-south axis, with the near side of the disk being toward the east, hence the 180° difference with the value of + 27.1° reported in Milli et al. (2017). One should also note that we do not constrain the direction in which the dust particles orbit around the star, and that the problem is symmetric.

3.4 Modeling strategy

Choosing the adequate scattering theory to compute the full (polarization and scattering) phase function when modeling debris disks observations still remains a challenge. The Mie theory (commonly used in the literature as it is computationally fast) seems to be insufficient to reproduce most of the spatially resolved observations (Lebreton et al. 2012; Rodigas et al. 2015; Olofsson et al. 2016; Milli et al. 2017). Other alternatives exist, such as the discrete dipole approximation (Purcell & Pennypacker 1973), but can be time consuming. Therefore, to alleviate this challenge, in this study we primarily focus on the radial profiles along the major axis of the disk. By doing so, and assuming that the disk is flat enough, we are probing the exact same scattering angles on both sides of the disk (close to 90°). Consequently, the exact value of the polarized phase function S₁₂ (integrated over all sizes) at this scattering angle does not matter when comparing the radial profiles along the northeast and southwest sides. The dependency of S₁₂(90°) as a function of the grain size remains, but given that observations in the optical are mostly sensitive to the small dust grains, which are dominating in number, we consider this effect to be of second order.

Nonetheless, modeling the radial cuts along the semi-major axis only does not allow us to really constrain the morphology of the disk. Preliminary results using solely radial cuts along the major axis led to clearly wrong results; the best fitting model would have a very large eccentricity (e ~ 0.2–0.3) and large reference radius (r₀ ~ 100–120 au) so that the major axis of the model would not be the same as in the observations but the radial cuts would intercept the disk at other azimuthal angles. Therefore, we considered two additional radial profiles at ϕ ± 10° from the major axis, to constrain the peak position of the disk, and obtain a more reliable determination of r₀, e, and ω. We chose not to consider additional cuts closer to the semi-minor axis as the signal-to-noise degrades quite significantly. Moreover, our study focuses on determining the radial density profiles, which are best constrained close to the semi-major axis. Figure 1 shows the location of those radial cuts.

The observations are de-rotated to align each of the different axes with the vertical direction. The radial profiles are then computed as the average of the polarized flux (from the Q_ϕ image) over a vertical slit centered on the central pixel, with a width of ± 5 pixels. For the uncertainties, we compute the standard deviation over the same vertical slit, from the U_ϕ image. For a given synthetic image, because of the finite number of particles that are used to generate the image, in some cases there can be some “shot noise”, with a given pixel being much brighter than its neighbors. This may lead to artificial local minima when trying to find the best fit model. To circumvent this issue, the model image is first smoothed by performing a median clipping over 3 × 3 neighboring pixels (the central pixel, which value is being estimated, is not included when computing the median). It is then convolved by a 2D Gaussian of standard deviation 1.22 × λ∕D, where D = 8 m the diameter of the telescope and λ = 0.735 μm. We then proceed similarly to the observations to extract the radial profiles, and scale them by finding the scaling factor that minimizes the difference between the northeast and southwest sides simultaneously. For the profile along the major axis (A_N and A_S on Fig. 1), the scaling factor is the same for both sides, as the polarized phase function should be the same. This allows us to test if the best fit model can successfully reproduce the brightness asymmetry between the north and south sides. For the radial cuts along B_N, B_S, C_N, and C_S, each scaling factor is determined independently as the polarized phase function is sampled at different angles (in principle, the phase function should be similar for the pairs C_N -B_S and B_N -C_S but we left the scaling factors unconstrained). Since the central region of the observations is contaminated by the instrumental PSF, we compute the goodness of fit between the ranges [−1.′′4, −0.′′95] and [0.′′ 95, 1.′′ 4] along the major axis, and [−1.′′3, −0.′′7] and [0.′′7, 1.′′ 3] for the other axis (these regions are highlighted by the wider white solid lines on Fig. 1).

The free parameters of the model are the reference radius r₀, the standard deviation of the radial distribution δ_r, the argument of periapsis ω, the standard deviation of the azimuthal collision probability at the pericenter δ_ω, the eccentricity e, and the position angle ϕ. The location of the star is fixed and is not allowed to vary in our modeling approach. As mentioned before the inclination is not a free parameter and is set at76.6°; the radial cuts are all close to the major axis of the disk, the worst direction to properly determine the inclination i. The grain size distribution exponent p is set to − 3.5. The grain size distribution is defined between s_min = 6μm (small enough that the grains in the first few bins are set on hyperbolic orbits) and s_max = 1.3 mm. The value of s_max is chosen so that the minimum value of β is 5 × 10⁻³ (optical constants of astrosilicates from Draine 2003, with a density of 3.5 g cm⁻³) and that the grain size distribution is sampled over a significant range of sizes (we set n_g = 200).

For the polarized phase function, we use the analytical Henyey-Greenstein expression, as (5)

where θ is the scattering angle and g the anisotropic scattering factor (− 1 ≤ g ≤ 1). Similar approaches are commonly used for modeling polarimetric images of debris disks (e.g., Engler et al. 2017, who also included a term for the diffraction peak caused by grains larger than the wavelength). Nonetheless, in the preliminary versions of Milli et al. (2019) the authors found that the polarized phase function at scattering angles close to 90° is overall well reproduced using g ~ 0.3⁴ (the final best fit being obtained using a combination of two Henyey-Greenstein functions, to best match the brightness close to the semi-minor axis, which we do not try to reproduce here). We therefore use this value of g = 0.3 throughout the paper, to alleviate the number of free parameters. Discussing the polarized phase function and the implications it has on the dust properties is out of the scope of this paper and we refer the reader to Milli et al. (2019). The choice of the analytical form of the polarized phase function also alleviates some of the free parameters regarding the dust properties (such as the porosity as well, see Arnold et al. 2019). Indeed, now that the phase function is parametrized for all sizes, the absolute value of β associated with each grain size matters less; what matters is the global shape of the β(s) curve as it will determine the radial distribution of the dust grains, which we are modeling. With our approach we can therefore ignore most of the dust properties⁵ and focus on the impact that the radiation pressure has on the spatial distribution of the small dust grains in the disk, regardless of their truesizes and properties.

To determine the best-fitting parameters and estimate their uncertainties, we then used an affine invariant ensemble sample Monte-Carlo Markov chain (emcee package, Foreman-Mackey et al. 2013). The initial conditions are set to be close to the disk parameters reported in previous studies (for r₀, e, ω, and ϕ). We used 30 walkers and ran a short burn-in phase (length of 400 models) and then run chains of 1000 models for each walker. The uniform priors are also reported in Table 1. At the end, the mean acceptance fraction is of 0.27 (indicative of convergence and stability, Gelman & Rubin 1992), and the maximum auto-correlation length among all the free parameters is 68 steps.

Fig. 1 Reduced ZIMPOL image of HR 4796 A. The super-imposed lines A to C show the locations where we measure the radial profiles (the width of the lines does not correspond to the width of the slits used to measure the radial profiles).

3.5 Results

The radial profiles are displayed in Fig. 2 and the observations and best-fit model are shown in Fig. 3, with the same linear stretch (the difference in surface brightness closer to the semi-minor axis being due to the choice of the polarized phase function, see Milli et al. 2019). The location of the pericenter of the disk is marked by a blue circle on the right panel of Fig. 3. While the best fit model slightly under-estimates the polarized intensity at separations larger than 1′′ along the B_N and C_N axis (still within 3σ), and the peak positions along the B_S and C_S axis are not a perfect match to the data (even though the signal-to-noise is lower in those regions), overall, the radial profiles along the major axis (A_N and A_S) are well reproduced, both for the peak positions and the slopes, up to 1.′′4. Figure A.1 shows a top view of the weighted cross section of the best fit model.

The most probable parameters are summarized in Table 1 and the probability density functions are shown in Fig. 4. The uncertainties for the MCMC results are estimated from the 0.16 and 0.84 quartiles using the corner package (Foreman-Mackey 2016). The projected posterior distributions are shown in Fig. 4. We find that the pericenter should be located on the front side of the disk, with . The reference radius of the disk is of au, and the standard deviation of the radial distribution of the parent belt is au, while the standard deviation of the collisional distribution is δ_ω = 63.9°. We find that the eccentricity is , and the position angle is − 152.1°, close to the value of − 151.6° we previously found to define the location of the major axis.

Fig. 2 Radial profiles of disk along three cuts highlighted in Fig. 1 (error bars are 1σ), with best-fitmodel over-plotted in red. The gray shaded areas show where the goodness of fit is not estimated.

Fig. 3 Observations and best-fit model to the SPHERE/ZIMPOL observations of HR 4796 A, with the same linear stretch (left and right, respectively). For the model, the blue circle marks the location of the pericenter.

Fig. 4 Projected probability density distributions along with the determined uncertainties for the different parameters in the modeling as well as density plots. The contours correspond to the [0.12, 0.39, 0.68, 0.86] density percentiles.

4 Discussion

The debris disk around HR 4796 A has been resolved at high angular resolution on several occasions, with different instruments (e.g., Lagrange et al. 2012; Rodigas et al. 2015; Milli et al. 2017; Schneider et al. 2018; Kennedy et al. 2018). All these observations showed that the disk appears as a narrow ring. With the new ZIMPOL observations presented here, we also find the disk to be very narrow. To reproduce the observations, our modeling results suggest that the parent planetesimal belt follows a normal distribution with a standard deviation that can be as narrow as 3.6 au. One should note that the radial extent of the parent belt is on the order of the vertical of the disk, suggesting it is shaped as a thin torus. From the results of Rodigas et al. (2014), Milli et al. (2017) concluded that the observed width of the disk around HR 4796 A could be explained by a planet lighter than Saturn, inward of the ring, shepherding the debris disk (or even smaller if the planet is migrating, e.g., Perez et al. 2019). Given the comparable angular resolution provided by the SPHERE/IRDIS and ZIMPOL instruments, we do not revise the values reported in Milli et al. (2017) to explain the narrowness and eccentricity of the parent planetesimal belt. But the ZIMPOL observations provide new insights into the azimuthal distribution of the dust grains, as well as on the production of small dust grains in this young debris disk.

4.1 The pericenter glow of HR 4796 A

The pericenter glow effect was originally proposed by Wyatt et al. (1999) who modeled Keck II observations of HR 4796 A at 18.2 μm (also presented in Telesco et al. 2000). The asymmetry observed in thermal emission could be explained by the fact that the dust grains in the direction of the forced pericenter of the disk are closer to the star, and hence warmer as they receive more stellar light. The pericenter glow can also be observed in scattered light observations, the dust grains at the pericenter receiving, and therefore scattering, more light compared to the apocenter. Wyatt et al. (1999) reported an argument of periapsis of 26° for an eccentricity of 0.02, meaning that the pericenter is located close to the projected semi-major axis of the disk. Nonetheless, from their Fig. 7 showing an ensemble of possible solutions, a forced eccentricity of 0.07 would lead to an argument of periapsis of about ~75° which, despite different reference system, is compatible with our results. However, Moerchen et al. (2011) modeled Michelle and T-ReCS mid-IR observations, and found that the pericenter should be located along the projected semi-major axis of the disk (perpendicular to the line of sight), but with an even larger eccentricity of 0.06 compared to Wyatt et al. (1999, 0.02 when the pericenter is located near the projected semi-major axis). They also found a family of solutions, with increasing eccentricity when the pericenter is moved closer to the projected semi-minor axis. However, since then it has been shown that this is not compatible with several studies that consistently found the pericenter to be located closer to the projected semi-minor axis of the disk, without increasing e too much (e.g., Thalmann et al. 2011; Rodigas et al. 2014; Milli et al. 2017, and this work). In Fig. A.2, we show the radial cuts along the semi-major axis of the disk, for a model in which particles are released uniformly in mean anomaly (all the other parameters being the same as the best model obtained before). One can note that in this model the brightness asymmetry is not well reproduced, justifying why our best fit model requires more particles to be released near the pericenter.

To further check the validity of our model, we computed mock observations at 18.1 and 24.5μm, to be compared with the observations presented in Moerchen et al. (2011). These mock observations are shown in Fig. 6; we used the same code as before, to compute images at the two wavelengths only considering thermal emission from the dust grains (using the same parameters as the best fit model to the optical observations), convolved them with 2D Gaussian with full width at half maximum of 0.′′ 52 and 0.′′ 72, and added noise to the images to mimic the real observations. While we did not aim at fitting those observations, visual inspection suggests that we are obtaining very comparable results, the brightness asymmetry being more pronounced at 18.1μm than at 24.5μm (where it is only marginally detected).

We additionally visually compared the ALMA observations presented in Kennedy et al. (2018) with a model computed using the same code at the same wavelength as the ALMA observations (thermal emission only). Figure 5 shows the Briggs-weighted Band 7 image (left) and the thermal emission for our best fit model, convolved by a 2D Gaussian similar to the beam of the observations. Our best fit model can reproduce the overall shape of the disk, the width of the ring, and the brightness distribution over all azimuthal angles. ALMA 880μm observations being sensitive to larger grains, which are not subjected to strong radiation pressure, bound dust grains close to the cut-off size do not contribute significantly to the thermal emission at those wavelengths and the model traces the distribution of the parent planetesimal belt. For those large grains, even though they are preferentially released near the pericenter, their eccentricities are rather small (close to the eccentricity of the parent belt), and therefore when populating the orbits, dust particles are distributed almost uniformly in the azimuthal direction. The point being that, in our model, the larger the grains the less significant the value of δ_ω becomes; ALMA observations can hardly constrain if dust grains are released in a preferential location of the disk (e.g., the pericenter in our case). Overall, this brings confidence to our best-fit model of the dust production and distribution in the disk around HR 4796 A, as it can reproduce (at least to the first order) observations from optical to millimeter wavelengths. In that model the pericenter glow effect plays a very minor role to explain the brightness asymmetry.

Fig. 5 Left: ALMA 880μm Briggs-weighted image of HR 4796 A. Right: best fit model at the same wavelength, convolved with a similar beam (displayed in the bottom left corner of both images). No noise was added to the model.

Fig. 6 Mock observations with the Michelle instrument (18.1μm, left panel) and T-ReCS (24.5μm, right panel).

4.2 Dust production around HR 4796 A

The azimuthal number density of an eccentric ring should naturally peak at the location of the apocenter (Pan et al. 2016), and not at the location of the pericenter as is the case for HR 4796 A according to our best-fit model. As discussed in Pearce & Wyatt (2014), if a planet is interacting with a debris disk less massive than the planet, the number density should be higher at the apocenter. Such a model would therefore not be applicable to the disk around HR 4796 A (around which very stringent upper limits on the presence of companions have been placed, Milli et al. 2017). On the other side, Pearce & Wyatt (2015) also investigated the case of the interactions between an eccentric planet and a debris disk of comparable mass. They conclude that the end-result of those interactions usually is a double-ringed debris disk, which is not observed for the disk around HR 4796 A.

A possible mechanism that could help explain our results would be the violent collision of large bodies at 76.4 au from the star (e.g., Kral et al. 2015). The numerical simulations presented in Jackson et al. (2014) indicate that all the bodies released from the original collision point will have to pass through the same location at each orbit, which would locally increase the collision probabilities. This collision point would therefore become the main production site for any secondary dust (or gas) in the system. Jackson et al. (2014) studied the cases of eccentric progenitors, and concluded that the resulting brightness asymmetry highly depends on where the collision took place. If the collision happened near the apocenter of the progenitors, then the brightness asymmetry is “constructive” due to the increased density at the collision point and the fact that particles spend more time near the apocenter (Fig. 11-D of Jackson et al. 2014). On the other hand, if the collision took place close to the pericenter of the progenitors, then there is a competition⁶ between the over-densities due to the “pinch-point” at the pericenter, and the more time spent by particles at the apocenter (Fig. 11-C of Jackson et al. 2014). Nonetheless, the apocenter being more spread out than the pericenter (especially for the bound grains with the highest β value), the latter one may still appear brighter than the former. Overall, the observations of the disk around HR 4796 A could be explained if the disk is the outcome of a unique, massive collision of initially eccentric progenitors (and the collision should have taken place close to the pericenter of the progenitors). Because the velocity dispersion is larger at the pericenter compared to the apocenter, a collision near the former could, in principle, generate a larger amount of small dust grains. Nonetheless, given how bright the disk is (f_disk = L_disk∕L_⋆ ~ 5 × 10⁻³), the collision would have to have been an extremely rare event. If this is indeed the case, the implications for the formation of large planetesimals at distances larger than 70 au are strong (see Kenyon & Bromley 2005, as well as the discussion about the collisional status of the disk in Kennedy et al. 2018). One would need to form at least two very massive oligarchs in the outermost regions, on a timescale of a few Myr. Indeed, the mass of the body whose breakup is able to produce a disk of debris of fractional luminosity 5 × 10⁻³ is extremely large. Taking the Eqs. (2)–(5) of Wyatt et al. (2007) linking f_disk to the mass of a collisional cascade producing it, we find that, even in the very optimistic hypothesis that all the mass is contained in ≤1m bodies, one needs at least a few Earth masses of material for the disk to be as bright as f_disk = 5 × 10⁻³ at 75 au from its central star (Augereau et al. 1999 had already found a similar estimate for the amount of ≤ 1m bodies). As a consequence, the catastrophic breakup scenario requires the breakup of planetary object, probably at least in the super-Earth range. One may furthermore wonder if such a collision should not also have released significant amount of gas (as postulated for β Pictoris for instance, Dent et al. 2014). Despite sensitive observations, no gas has been detected in the outermost regions of HR 4796 A (Kennedy et al. 2018), but Iglesias et al. (2018) reported possible “falling evaporating bodies”, by detecting variable extra absorption lines in optical spectroscopic observations.

4.3 The effect of HR 4796 B on the debris disk

To estimate if the M-type star HR 4796 B can reasonably have an effect on the disk around HR 4796 A (projected separation of 568.3 au), through radiation pressure, stellar winds, or gravitational interactions (e.g., Thébault et al. 2010; Cuello et al. 2019), we first check the separation between both stars. We used their Gaia DR2 measurements (Gaia Collaboration 2018), which are reported in Table 2. The radial velocity of HR 4796 A was taken from Iglesias et al. (2018), and we estimated the one of HR 4796 B from UVES observations (program IDs 082.C-0218 and 089.C-0207). The first noteworthy difference is the parallaxes of both stars, which translate to distances of 71.9 and 70.9 pc for HR 4796 A and B, respectively. However, the measurements for HR 4796 A have large uncertainties due to its brightness and the star is flagged for possible astrometric errors. Therefore, we here assume that the B star has the same parallax as the A star and to evaluate if the former can have an impact on the disk around the A star, we checked if the system is bound. To that end, we estimated the escape velocity of B with respect to A, as , where G is the gravitational constant, M_⋆,A is the mass of the HR 4796 A (1.3M_⊙), and r the separation between the two stars. With the positions and velocities of both stars, we find an escape velocity of 2.01 km s⁻¹, compared to a relative velocity of 2.47 km s⁻¹ (estimated from the proper motions and radial velocities of both stars). As we assumed the same parallax for both stars, this is the most favorable case with the smallest three-dimensional separation. Therefore, with the available astrometric measurements, it seems that HR 4796 A and B are not bound to each other, most likely minimizing the possible impact that the B component can have on the debris disk, but this may have to be revisited with a more reliable astrometry for HR 4796 A in the near future.

4.4 An asymmetric halo of small dust grains

Schneider et al. (2018) presented deep HST observations of the disk around HR 4796 and revealed an extremely extended halo outside of the birth ring. They detect this “exo-ring” material up to ~ 10′′ along the northeast side, while the southwest side appears more compact. As a matter of fact, there seems to be an “arc-like” feature along the southwest side, which could be due to either interactions with the ISM gas, or with HR 4796 B. Based on the discussion above, the latter scenario is probably unlikely, if HR 4796 A and B are indeed not gravitationally bound. The authors mention the possibility that the dust grains in the halo may be unbound from the central star, but those dust grains should be evacuated from the system extremely quickly (see e.g., Thébault & Wu 2008), and therefore should hardly be detected. On the other hand, bound grains, with large β ratios can have their apocenter at 10–20 times the separation of the birth ring. Therefore, in the following, we attempt to reproduce the general shape of the exo-ring (extended emission along the northeast with an arc-like shape along the southwest), based on our current best-fit model, and we consider bound dust grains with high β ratios. This subsection is meant to be speculative, given the complexity of the problem and we simply aim at providing two scenarios to explain the HST observations.

First of all, our current best fit model, out-of-the-box cannot reproduce the HST observations. We find that the pericenter is located along the north side, and therefore, the disk is naturally more extended in the direction toward the apocenter, that is, along the south side, while the HST observations show otherwise. We therefore consider here two alternative explanations: first, the dust grains set on high eccentricity orbits are interacting with the local interstellar medium, and second, a slow precession of the pericenter of the debris disk is causing the apparent asymmetry.

4.4.1 Description of the code

When considering additional forces, such as gas drag, the analytical prescriptions given in Eq. (3) are no longer useful to estimate the orbital parameters of the dust grains. Therefore, we opted to use a simple Runge–Kutta integrator at the fourth order to estimate the positions and velocities of the particles. The forces considered here are the gravitational potential from the central star (weighted by 1 − β when considering the effect of radiation pressure), and the gas drag on each dust grain. We followed the work of Marzari & Thébault (2011; see also Pástor 2017) for the implementation, and the acceleration felt by a dust grain is estimated as (6)

where v is the velocity of the dust grain, v_H the velocity of the ISM gas, n_H and m_H are the density and mass of the hydrogen atoms in the ISM, and C_D is a drag coefficient (set to 2.5 as in Marzari & Thébault 2011). The cross section of a dust grain is given by A_s = πs².

To produce an image, we initially release 10 000 particles (uniformly distributed in mean anomaly) within the birth ring (described by its semi-major axis, eccentricity, and argument of pericenter, which are the same as the best fit model) and follow their trajectories in time. For the sake of simplicity, in those simulations there is no preferential collisions near the pericenter. At the initialization of the simulation, we first check if each particle is indeed bound to the system, by estimating its initial velocity, and comparing it to the escape velocity of the system. Only particles that are bound are kept in the simulation. At each time step (usually four years) we save the (x, y, z) positions of each of the particles, and then project those values onto the sky plane, depending on the inclination, position angle, and opening angle of the disk. For all the simulations presented in the following, we assume the Henyey-Greenstein approximation for the phase function, with g = 0 (isotropic scattering), and compute the surface brightness by multiplying the phase function by Q_scaπs²∕(4πr²). Also, we only consider a single grain size of s = 7.15μm, as this is the smallest size for which all the particles remain bound to the star, and have a high β ratio of ~ 0.46.

4.4.2 Interaction with local ISM gas

The simulations presented in this section last for 100 000 yr (25 000 steps of 4 yr), and the free parameters that we investigate are: the density of hydrogen atoms in the ISM n_H, the direction and velocity of the ISM gas (all encompassed in the v_H vector), and the direction of rotation of the dust grains in the disk. To simplify the definition of the problem, we assume that the v vector only represents the orbital velocity of the particles, while the v_H vector contains information about both the proper motion of the star and the direction of motion of the ISM gas. As mentioned before, a full exploration of the parameter space is out of the scope of this paper, and we simply aim at providing a qualitative assessment on how to explain the HST observations.

Figure 7 shows two simulations as well as the HST observations published in Schneider et al. (2018; left panel), where the only difference is the rotation direction of the dust grains in the disk (clockwise and counter-clockwise for the middle and right panels, respectively). The images have been convolved by a 2D Gaussian with a standard deviation of 1 pixel, and are shown with a log stretch to highlight the faint outer regions (as a consequence, the birth ring appears broader than in ZIMPOL simulations which are shown with a linear scale). For those simulations, we tried to fix as many free parameters as possible. We therefore assumed that the radial velocity of the ISM gas matches the one of the central star (Iglesias et al. 2018 detected several absorption lines at different velocities, but we cannot know the relative distances of those clouds, only that they are between the star and the observer). Furthermore, given that the arc-like shape is in the east-west direction, we assumed that the ISM gas has a null velocity in the north-south direction, and that the δ component of the v_H vector is equal to the proper motion of the star (see Table 2). This leaves us with α, n_H, and the direction of rotation of the dust grains in the disk as free parameters. Similarly to Pástor (2017, studying the disk around HD 61005), we find that the ISM density has to be quite significant to produce the arc-like feature in the southwest direction (but overall, the density and the amplitude of the velocity are degenerate parameters). For the simulations presented in Fig. 7, we set n_H = 125 cm⁻³ and α = −100 mas.yr⁻¹. This means that in those simulations, the star is moving 55 mas yr⁻¹ toward the west (following its proper motion in right ascension) and that the ISM gas is moving 45 mas yr⁻¹ toward the east. Overall, both simulations can produce the arc-shape structure as seen with HST, roughly at the same separation (~ 5′′), but when the dust grains are moving in the clockwise direction onto the sky plane, the northeast side of the disk is much fainter, and the south side of the birth ring appears brighter than the north side. This suggests that if the extended halo is indeed the consequence of interactions with the ISM gas, then, the dust grains would most likely be orbiting counter-clockwise around HR 4796 A. The main caveat in this scenario is the rather large density of the ISM gas required, compared to the surroundings of the solar system, similarly to what was found for HD 61005 by Pástor (2017). However, the volume density of the cold, dense, interstellar medium can vary between 5 and 100 cm⁻³ (e.g., Nguyen et al. 2019 and references therein). Therefore, our results, while on the higher end of the range, remain overall compatible with studies of the interstellar medium.

Fig. 7 HST observations, as published in Schneider et al. (2018) (left panel), and mock HST observations calculated from the N-body simulations of small, bound, dust grains around HR 4796 A, when considering interactions with the local ISM (clockwise and counter-clockwise rotations on the middle and right panels, respectively). The images are in log-scale.

4.4.3 Precession of the pericenter

The second scenario we investigate to explain the HST observations is based on the fact that an eccentric disk is more radially extended toward the apocenter. Several studies of the disk around HR 4796 A consistently found that the pericenter of the disk is located close to the (projected) semi-minor axis of the disk, on the north side. But, if the disk is precessing (see also Löhne et al. 2017), in the past the pericenter might have been located on the south side, resulting in a more extended disk along the north side. To test this hypothesis, we run similar N-body simulations, for a single grain size, without any additional forces (only gravitational and radiation pressure forces). The two free parameters of the simulations are the total duration t = t_sim and the precessing rate Ω of the pericenter (i.e., (ω_end − ω_start)∕t_sim). In this case, we assume that the pericenter of the disk is rotating counter-clockwise, from the southwest toward the northwest (just past the semi-minor axis of the disk). Therefore, we only consider dust grains that are also rotating counter-clockwise around the star.

We chose ω_start = −320°, ω_end = −254°, and a total duration t_sim = 200 000 yr, divided in n_sim = 500 different steps. At t = 0, the pericenter is located at ω_start, and the N-body simulation last over t_sim, with a timestep of four years and we save the position of all the particles for the last 20 iterations (hence the last 80 yr). The ith N-body simulation (out of the 500) will last slightly less time (t_sim − i × t_sim∕n_sim) and the location of the pericenter will have moved slightly (ω_start + i × (ω_end − ω_start)∕n_sim), and we still save the last 80 yr of the simulation. At the end, we simply collapse all the n_sim images together (see Thébault 2012 for a similar approach).

Figure 8 shows the result of this simulation (on an aggressive log stretch to reveal the outermost regions, also rendering the birth ring thicker than in other figures). One can see that indeed the southwest side is slightly dimmer than the northeast side, less extended in the radial direction, and that arc-like structures start to develop in the southwest side. The concentric ellipses close to the birth ring are due to the fact that we save the positions of the particles over the last 80 yr of the simulations, but the time step between each N-body simulations is larger than 80 yr (400 yr). It also appears, that we cannot reproduce the scale of the HST image in this simulation (the arc shape is located at about 5′′ in the HSTobservations, while the disk still appears bright at even 10′′). This would most likely suggest that we are observing grains with a β value slightly smaller than 0.46.

This approach is extremely simplistic, and heavily depends on the initial conditions and the duration of the simulation (i.e, when we decide to stop the simulation). The main caveat is that we assume that the dust particles can actually survive over t_sim without being destroyed. Ifwe assume β = 0.45, e = 0.07, a = 76.4 au, and draw uniformly 1000 values for the true anomaly between [0, 2π), we can estimate the orbital parameters of the dust grains following Eq. (3). We can then estimate the orbital period of those grains as . We find a mean value of 18 500 yr with a standard deviation of 15 800 yr (the distribution is strongly peaked at ~5 500 yr). This means, that over t_sim = 200 000 yr, the dust grains would, on average, pass trough the birth ring 10−40 times (this criterion, combined with the CPU time for the simulation, drove the choices for the precession rate and the total duration of this simple exercise). Depending on the optical depth of the birth ring, those grains may get destroyed before that, which may work in our favor. Indeed, to produce an arc-like feature at large distances from the star, one needs to break the symmetry. If the disk has been precessing for a long time, and none of the particles are destroyed, then the result would most likely be symmetric. However, if after several orbits, the oldest dust grains are destroyed in the birth ring, this can generate a possible asymmetry in the surface brightness of the disk, as we “lose memory” from past events.

Overall, our simulation does not properly match the observations, and remains speculative, but it seems to be going in the right direction to try and explain the HST observations presented in Schneider et al. (2018). Those observed images could well be result of a combination of both effects, precession of the disk and interactions with local interstellar medium. Furthermore, the scenario that we invoked earlier to explain the morphology of the disk (a collision at 76.4 au, following the work of Jackson et al. 2014), is not incompatible with precession of the disk. Jackson et al. (2014) mention that even thoughthe original collision point remains static on orbital time-scales, it should still be able to precess, due for instance to the presence of other massive bodies in the vicinity of the birth ring. Further investigation of the lifetime of the dust grains, precession rate, and the reason for the precession of the disk is out of the scope of this paper, the intent here being simply to propose a scenario that can reconcile all the available observations of the disk.

Fig. 8 Mock HST observations calculated from the N-body simulations of small, bound, dust grains around HR 4796 A, when considering that the pericenter is precessing over time. The image is in log-scale.

5 Conclusion

We presented high angular resolution SPHERE/ZIMPOL observations that reveal the bright debris disk around HR 4796 A with exquisite angular resolution (see also Milli et al. 2019 for a discussion about the polarized phase function). At optical wavelengths, the northeast side of the disk is brighter than the southwest side, which we aimed at reproducing in this paper. We modeled the radial profiles along (and close to) the semi-major axis of the disk with a code that includes a simple prescription of the effect of radiation pressure. With this code, which is faster (but also much simpler) than models including a more accurate treatment of collisional and dynamical evolution of debris disks (e.g., Kral et al. 2013, 2015; Löhne et al. 2017), we are able to reproduce the observed profiles on both sides of the disk. We can reproduce the outer edge of the disk without invoking the presence of an outward planet truncating the disk (similarly to Thebault et al. 2012) and we find that the underlying planetesimal belt can be as narrow as a few au. As previously stated in the literature, this could be related to the presence of an unseen planet, inward of the planetesimal belt, shepherding the debris disk. We find, similarly to other studies in the past years, that the pericenter of the disk is located close to the projected semi-minor axis of the disk. We show that with such a configuration, the pericenter glow, that had been postulated to explain marginally resolved observations, has in fact very little impact on the azimuthal brightness distribution of the disk. To reproduce the observed brightness asymmetry, we find that small dust grains must be preferentially released, as a result of collisions between larger bodies, close to the pericenter of the disk. Finally, our best-fit model can self-consistently reproduce most of the available observations of the disk, from optical to millimeter wavelengths. The only dataset that remains challenging to explain are the recently published HST observations that reveal an extended halo at large separations from the star. After concluding that HR 4796 B may not be bound with HR 4796 A, thus minimizing its possible effect on the disk, we propose two possible scenarios that could help explaining the HST dataset; interactions with the local interstellar medium and precession of the disk. Even though the results remain speculative, those hypotheses (or a combination of those) could help reproduce the very faint halo of small dust grains. But further investigation is required to confirm our findings, for instance by performing more refined dynamical simulations (e.g., accounting for grain-grain collisions during the evolution of an eccentric disk, Löhne et al. 2017)

Acknowledgements

We thank Glenn Schneider for kindly providing the HST observations. We thank the anonymous referee for providing useful comments that helped improving the paper, especially the description of the code. This research has made use of the SIMBAD database (operated at CDS, Strasbourg, France) and of the Spanish Virtual Observatory (http:// svo.cab.inta-csic.es, supported from the Spanish MICINN/MINECO through grants AyA2008-02156, AyA2011-24 052). This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013), as well as the TOPCAT software (Taylor 2005). SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF – Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016). We also acknowledge financial support from the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU in France. This work has also been supported by a grant from the French Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56). The project is supported by CNRS, by the Agence Nationale de la Recherche (ANR-14-CE33-0018). It has also been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). M.R.M., H.M.S., and S.D. are pleased to acknowledge this financial support of the SNSF. Finally, this work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice) and Observatoire de Paris/LESIA (Paris). We thank P. Delorme and E. Lagadec (SPHERE Data Centre) for their efficient help during the data reduction process. J.O., A.B., J.C.B., D.I., M.M., M.R.S., and C.Z. acknowledge support from the ICM (Iniciativa Científica Milenio) via the Nucleo Milenio de Formación planetaria grant. J.O. acknowledges support from the Universidad de Valparaíso and from Fondecyt (grant 1180395). A.B. acknowledges support from Fondecyt (grant 1190748). M.M. acknowledges financial support from the Chinese Academy of Sciences (CAS) through a CAS-CONICYT Postdoctoral Fellowship administered by the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. FMe acknowledges funding from ANR of France under contract number ANR-16-CE31-0013. J.C.B. acknowledges support from Proyecto FONDECYT postdoctorado 2018 nro. 3180716. G.M.K. is supported by the Royal Society as a Royal Society University Research Fellow. A.Z. acknowledges support from the CONICYT + PAI/Convocatoria nacional subvención a la instalación en la academia, convocatoria 2017 + Folio PAI77170087.

Appendix A Miscellaneous

Figure A.1 shows the de-projected, weighted cross section, of the best fit model to the ZIMPOL observations, also displaying how we define the argument of periapsis and the position angle of the disk. Figure A.2 shows the cut along the semi-major axis (profiles A_N and A_S) for a model in which the particles are not released from a preferential location in the disk. All the parameters of the model are the same as the best fit model to the ZIMPOL observations otherwise. This demonstrates that for our best fit solution, pericenter glow alone cannot fully reproduce the brightness asymmetry observed along the major axis of the disk. The other radial cuts are not shown and the main point of the figure is to discuss the contribution of pericenter glow to the brightness distribution.

Fig. A.1 Top view of the weighted cross section of the best fit model, indicating the position angle and argument of periapsis.

Fig. A.2 Cut along the semi-major axis of the disk when particles are released uniformly in mean anomaly in the disk (i.e., no preferential collision at the pericenter).

References

All Tables

All Figures

	Fig. 1 Reduced ZIMPOL image of HR 4796 A. The super-imposed lines A to C show the locations where we measure the radial profiles (the width of the lines does not correspond to the width of the slits used to measure the radial profiles).
In the text

	Fig. 2 Radial profiles of disk along three cuts highlighted in Fig. 1 (error bars are 1σ), with best-fitmodel over-plotted in red. The gray shaded areas show where the goodness of fit is not estimated.
In the text

	Fig. 3 Observations and best-fit model to the SPHERE/ZIMPOL observations of HR 4796 A, with the same linear stretch (left and right, respectively). For the model, the blue circle marks the location of the pericenter.
In the text

	Fig. 4 Projected probability density distributions along with the determined uncertainties for the different parameters in the modeling as well as density plots. The contours correspond to the [0.12, 0.39, 0.68, 0.86] density percentiles.
In the text

	Fig. 5 Left: ALMA 880μm Briggs-weighted image of HR 4796 A. Right: best fit model at the same wavelength, convolved with a similar beam (displayed in the bottom left corner of both images). No noise was added to the model.
In the text

	Fig. 6 Mock observations with the Michelle instrument (18.1μm, left panel) and T-ReCS (24.5μm, right panel).
In the text

	Fig. 7 HST observations, as published in Schneider et al. (2018) (left panel), and mock HST observations calculated from the N-body simulations of small, bound, dust grains around HR 4796 A, when considering interactions with the local ISM (clockwise and counter-clockwise rotations on the middle and right panels, respectively). The images are in log-scale.
In the text

	Fig. 8 Mock HST observations calculated from the N-body simulations of small, bound, dust grains around HR 4796 A, when considering that the pericenter is precessing over time. The image is in log-scale.
In the text

	Fig. A.1 Top view of the weighted cross section of the best fit model, indicating the position angle and argument of periapsis.
In the text

	Fig. A.2 Cut along the semi-major axis of the disk when particles are released uniformly in mean anomaly in the disk (i.e., no preferential collision at the pericenter).
In the text