© ESO 2021

1 Introduction

Debris disks are the natural by-product of the planet formation process. They are second-generation dusty circumstellar disks created by collisions of planetesimals that have formed in previously existing planet-forming disks (Kral et al. 2018; Hughes et al. 2018). This process releases a large amount of micron-sized dust grains whose presence can be inferred either through the infrared excess over the photospheric level of their host star or from spatially resolved scattered (linearly polarized) light observations in the optical or near-infrared. The average lifetime of those small dust grains, however, is much shorter than the typical age of the star, and thus they have to be replenished continuously from the larger bodies. The most efficient process for removing these micron-sized dust grains from the system is usually the pressure exerted by the radiation field of the central star (e.g., Krivov 2010). For low-mass stars, however, the irradiation field is generally much weaker, while stellar winds can be strong (see, e.g., Reidemeister et al. 2011). Thus, stellar winds are most likely the dominant, yet poorly constrained, mechanism driving the rapid removal of particles. The issue has not yet been sufficiently addressed because debris disks around low-mass stars are rarely detected.

Debris disks are found around about 20% of A-type stars (Su et al. 2006; Eiroa et al. 2013; Chen et al. 2014; Matthews et al. 2014) and are mostly detected through the excess emission in the far-infrared. Similar to the searches targeted at solar-type and more massive stars, several large surveys have been conducted to search for cold debris disks around M dwarfs in the mid- and far-infrared (Gautier et al. 2007; Avenhaus et al. 2012; Kennedy et al. 2018) as well as in the submillimeter domain (Lestrade et al. 2006, 2009). These and other surveys yielded numerous debris disk candidates around young M stars (for details see, e.g., Luppe et al. 2020, and references therein). However, up to now, only five disks around young M stars have been confirmed based on more than one independent observation, such as spatially resolved imaging (in scattered light and/or the submillimeter to millimeter domain): AU Mic (Kalas et al. 2004; Boccaletti et al. 2015), TWA 7 (Choquet et al. 2016; Olofsson et al. 2018; Bayo et al. 2019; Esposito et al. 2020; Ren et al. 2021), TWA 25 (Choquet et al. 2016), GJ 581 (Lestrade et al. 2012), Fomalhaut C (LP 876-10; Cronin-Coltsmann et al. 2021), and GSC 07396-00759 (Sissa et al. 2018). Given the low number of spatially resolved detections, we still know very little about debris disks around low-mass stars, and yet several of these disks show very interesting and peculiar features, such as possible spiral-like structures, the gas around TWA 7 (Olofsson et al. 2018; Matrà et al. 2019), and the fast-moving arc-like structures around AU Mic (Boccaletti et al. 2018). Debris disks around low-mass stars are also very interesting targets, especially with respect to the potential connection between their occurrence and the possible presence of planets (Raymond et al. 2011). Although only about 2% of the M stars have been found to host giant planets (Johnson et al. 2007), rocky planets appear to be more frequently detected around low-mass stars, as pointed out by Dressing & Charbonneau (2015), who estimated a cumulative planet occurrence rate of 2.5 ± 0.2 planets per M dwarf with radii 1–4 R_⊕ and periods shorter than 200 days.

However, the reason for the paucity of debris disk detections around low-mass M-type stars remains unclear (see, e.g., Luppe et al. 2020)since stars of all spectral types appear to have a similar detection frequency as protoplanetary disks in the earlier stages of their evolution (see, e.g., Andrews & Williams 2005). In fact, studies of the Lambda Orionis star-forming region carried out on a wide set of spectroscopically confirmed members of the central star cluster Collinder 69 (~ 5–12 Myr) by Bayo et al. (2012), or the ALMA 887 μm survey of the disk population in the nearby 2 Myr old Chamaeleon I star-forming region conducted by Pascucci et al. (2016), indicate that disks around young, low-mass stars with M_⋆ ≲ 0.6M_⊙ are in fact morefrequent than those around higher-mass hosts. Therefore, it is more likely that current observations may simply not be sensitive enough because dust experiences significantly less heating around low-luminosity M dwarfs than around more massive stars. Consequently, these colder disks are significantly more difficult to detect at the typical wavelengths (e.g., 24–160 μm) that were used to build statistics of debris disks, making the excess emission hard to detect in the spectral energy distribution (SED) of these stars (Morey & Lestrade 2014). This might increase the chances of missing potential targets when compiling a sample for an observation or survey, and better suited alternatives are needed.

One alternative, as mentioned, is high angular resolution imaging observations in total intensity or polarized light, at optical or near-infrared wavelengths. Pioneered by Hubble Space Telescope (HST) observations (e.g., Graham et al. 2007; Stark et al. 2014), the availability and the recent advances in high contrast imaging instruments, such as the Gemini Planet Imager (GPI; Macintosh et al. 2006) or the Very Large Telescope/Spectro-Polarimetric High-contrast Exoplanet REsearch (VLT/SPHERE; Beuzit et al. 2019), provide direct access to the absorption and scattering properties of the grains and have opened new avenues to detect, resolve, and investigate debris disks at high angular resolution. Studying how stellar light is scattered off of the dust grains, either using the color of the disk between different bands (e.g., Debes et al. 2008; Rodigas et al. 2015) or through studying the phase function over a wide range of scattering angles – that is, the angle between the star, the dust grain, and the observer (e.g., Olofsson et al. 2016; Milli et al. 2017, 2019; Ren et al. 2019) – allows better constraints to be put on the properties of the dust, such as their typical grain sizes as well as their porosity, shape, and composition.

The target of our study is GSC 07396-00759, a young, nearby (71.43 ± 0.26 pc; Gaia Collaboration 2018) M1-type star. The star, classified as a weak-line T-Tauri star (Kastner et al. 2011), is probably a member of the β Pictoris Moving Group (β Pic MG, ≈ 18 ± 2 Myr; Miret-Roig et al. 2020), which is known to harbor numerous debris disk host stars, such as AU Mic and CP-72 2713 (Moór et al. 2020). Total intensity observations of GSC 07396-00759 obtained by Sissa et al. (2018) have revealed an extended (r₀ = 70 au) and nearly edge-on disk (i = 83°), probably containing sub-micron-sized grains, a possible indicator of a strong interaction of the stellar radiation field with the disk. By comparing their best fitting model with the observed SED of GSC 07396-00759 up to 22 μm (WISE/W4, <4.49 mJy; Cutri et al. 2014), the longest wavelength available at this point, Sissa et al. (2018) estimated an upper limit for the dust mass in the disk around GSC 07396-00759 of M_dust ~ 0.33 M_⊕, corresponding to an upper limit for the fractional luminosity of L_disk∕L_⋆ ≤ 4 × 10⁻³ due to the lack of far-infrared photometry data points. The disk furthermore displays a small brightness asymmetry, swept-back wings (warps), and ripples in the spine of the disk on both sides of the disk. The available proper motion and radial velocity data suggest that GSC 07396-00759 is also likely associated with V4046 Sgr AB, forming a loosely bound hierarchical multiple system (a_proj. ≈ 0.06 pc; Torres et al. 2006; Kastner et al. 2011; Sissa et al. 2018). V4046 Sgr AB is itself a close binary with accretion signatures (Stempels & Gahm 2004) and a gas-rich circumbinary disk (Rosenfeld et al. 2013; Rapson et al. 2015).

In this paper, we aim to compare newly obtained near-infrared polarimetric observations at high angular resolution of GSC 07396-00759 with radiative transfer modeling to study the morphology of the disk and to probe dust properties under the influence of radiation pressure and stellar winds, very rare constraints for low-mass stars and thus far only determined for a few debris disks, such as AU Mic (Augereau & Beust 2006) or ε Eri (Reidemeister et al. 2011). In Sect. 2, we describe the observations obtained with the SPHERE InfraRed Dual-band Imager and Spectrograph (IRDIS) instrument at near-infrared wavelengths. The model used to analyze the observationsis presented in Sect. 3, followed by the results reported in Sect. 4 and their discussion in Sect. 5. We conclude with a short summary and our conclusions in Sect. 6.

2 Near-infrared polarimetric imaging

2.1 Observations and data reduction

The observation of GSC 07396-00759 took place on June 22, 2018 (UTC; Program ID: 1100.C-0481(R), P.I.: J.L. Beuzit) as part of the SPHERE Guaranteed Time Observations. The data were obtained with the SPHERE/IRDIS (pixel scale of 12.25 mas, ~ 11′′ × 11′′ field of view; Dohlen et al. 2008) using the field-stabilized, dual-beam polarimetric imaging (DPI) mode (Langlois et al. 2014; de Boer et al. 2020; van Holstein et al. 2020), and employing the H-band (BB_H) filter with a central wavelength λ_c = 1625 nm and a width Δλ = 290 nm. To further increase the contrast the Apodized Lyot Coronagraph (mask diameter: 185 mas, Carbillet et al. 2011; Guerri et al. 2011)was used to mask the central star (Hmag = 8.76, Rmag = 12.01, Cutri et al. 2003; Zacharias et al. 2017). This allows for an efficient suppression of the stellar light, which is assumed to be only marginally linearly polarized, while keeping the scattered (i.e., polarized) light from the circumstellar disk relatively unaffected.

Each polarimetric observation consists of a set of four linear-polarization components, called Stokes Q⁺, Q⁻, U⁺, and U⁻, that are obtained by subtracting the two beams with orthogonal polarization states recorded simultaneously on the detector and tuning their polarization direction with a half-wave plate (HWP) with positions of 0°, 45°, 22.5°, and 67.5°, respectively. We took one exposure for each of the Stokes Q⁺, Q⁻, U⁺, and U⁻ components, each with a detector integration time (DIT) of 64 s. The polarization cycle of Q⁺, Q⁻, U⁺, and U⁻ was then repeated twelve (N_pol) times, adding up to a total of 48 exposures with a total integration time (t_exp) of 51.2 min. An overview of our observations is presented in Table 1. We also list the observing conditions at the time of observation, such as the average coherence time (τ₀), the average seeing, estimated from the Differential Image Motion Monitor (DIMM), the airmass (Am), and the Strehl ratio in the H-band, measured from the flux image.

Since the star is obscured by the coronagraphic mask in the science images, a flux and center calibration frame were taken in addition to the science observations. The flux calibration frames were taken with the central star moved away from the coronagraphic mask. The center calibration frames were taken after the star was aligned behind the coronagraphic mask and with the deformable mirror system (SAXO; Fusco et al. 2006) used to introduce a waffle pattern to create equidistant calibration spots outside of the coronagraphic mask. This allowed us to accurately determine the stellar position behind the coronagraph from the center image while we compute the Strehl ratio in the H-band for GSC 07396-00759 directly from the flux image with the star off-centered from the coronagraph. The H-band Strehl ratios estimated by the adaptive optics system during the observation are recorded in separate FITS files¹ (called GEN-SPARTA data). From these data we find Strehl ratios for GSC 07396-00759 ranging from 0.6 to 0.77, and in good agreement with our measured Strehl ratio of 0.67 from the flux image.

The data were reduced using the IRDAP² (IRDIS Data reduction for Accurate Polarimetry, van Holstein et al. 2020; de Boer et al. 2020) pipeline. IRDAP is a dedicated pipeline for the reduction of polarimetric data obtained with IRDIS, capable of differentiating and correcting instrumental and stellar polarization. For a detailed description of the reduction procedure and a discussion of the applied corrections we refer the reader to van Holstein et al. (2020) and de Boer et al. (2020). The reduction can be summarized as follows.

After applying the standard calibration routines, including sky-frame subtraction, flat-fielding and bad-pixel correction, the images are split into two individual frames representing the left and right sides of the IRDIS detector, corresponding to the parallel and perpendicular polarized beams, respectively. Then, the precise position of the central star is measured using the star center calibration frames on both image sides separately, and the right side of the image is shifted to a common center and subtracted from the left side. The pipeline then applies the double-difference method (see, e.g., Tinbergen 1996) to obtain the linear Stokes parameters Q and U corrected for instrumental polarization created downstream of the HWP and the corresponding total-intensity images computed with the double sum, respectively. However, the double difference does not remove instrumental polarization caused by the telescope and instrument mirrors upstream from the HWP, which is assumed to be proportional to the total intensity image, as shown in Canovas et al. (2011), nor does it remove the most important cross talk contributions (de Boer et al. 2020). IRDAP uses a Mueller matrix model to determine the polarimetric response function for the polarimetric imager (de Boer et al. 2020)and to correct for these instrumental polarization effects. This model describes the complete optical path of SPHERE/IRDIS (i.e., the telescope and instrument) and has been fully validated with measurements using SPHERE’s internal source and observations of unpolarized standard stars (van Holstein et al. 2020). The images of Stokes Q and U incident on the telescope are computed by setting up a system of equations describing every measurement of Q and U and solving it – for every pixel individually – using linear least-squares. The Q and U images thus created, however, may still contain some stellar polarization that is constrained by measuring the flux in the Q and U images around regions that should be virtually devoid of polarized signal from the disk. The final product of the reduction pipeline are the images of the azimuthal Stokes parameter Q_ϕ and U_ϕ (for definition see de Boer et al. 2020), where Q_ϕ > 0 is equivalent to a azimuthal polarization component (with respect to the position of the star), Q_ϕ < 0 to a radial component, and ± U_ϕ signal is equivalent to polarization angles oriented at ± 45° with respect to the azimuthal component, respectively.

2.2 Observational results

Figure 1 shows the final collapsed total-intensity, Q_ϕ and U_ϕ images of our IRDIS DPI H-band observations of GSC 07396-00759. As can be seen, the total-intensity image (left panel) shows numerous point sources, and in the field of view (FoV; cropped to ~ 2.2′′ × 2.2′′) of our IRDIS observation alone, we identified more than ten different sources. In fact, Sissa et al. (2018) detected a total of 109 point sources in their dual-band imaging observations taken with IRDIS as part of the SHINE survey (Vigan et al. 2021), in the night of June 15, 2017. However, by using observations of GSC 07396-00759 taken in the International Deep Planet Survey (IDPS; Vigan et al. 2012; Galicher et al. 2016) with the NIRC2 instrument at Keck, and with a time difference of about 10 yr, they were able to identify 70 objects as background objects based on common proper motion. Another 32 sources could be rejected as companions based on their location in the color-magnitude diagram (CMD). For the seven remaining objects the results could not be determined unambiguously; however, all but one have separations above 5′′ (a_proj. > 400 AU), and the colors of those objects are similar to that of the identified background objects (Sissa et al. 2018). It is thus highly likely that these are unrelated background objects too.

The middle-panel of Fig. 1 shows the reduced Q_ϕ image. Noticeable is the almost edge-on diskextending about 1.3′′ to the southeast (SE) and northwest (NW). Similar to the findings by Sissa et al. (2018), we detect an asymmetric brightness distribution,with the SE side of the disk appearing slightly brighter than the NW side. The disk also appears to be warped on the NW side, indicating a more complicated structure of the disk. The image however shows one of the major advantages of polarimetricobservations. That is, since stellar light is usually not polarized, the image is not contaminated by background stars, which makes it easier to detect circumstellar matter. In fact, except for one source or artifact located about 0.5′′ north of the central star, the disk appears to be free of contamination from other sources visible in the total intensity image. We masked this source separately during our subsequent analysis to minimize its contribution.

We also note that on average the U_ϕ image (Fig. 1, right panel) is mostly devoid of signal. This is consistent with what we expect in case of single scattering in a centrally illuminated, optically thin disk. We, thus, assume that the U_ϕ image to the first order contains only noise and measure these uncertainties from the standard deviation in concentric annuli with a width of 1 pixel in the U_ϕ image to create a noise map. This allows us to estimate the goodness of fit for our models as well as to create the signal-to-noise ratio (S/N) map shown in Fig. 2, computed as the ratio of Q_ϕ image and the estimated noise map. To minimize the influence of the higher noise close to the central star on the modeling, we introduce a circular numerical mask with r = 0.3′′ as inner boundary as well as an ellipse as outer limit for the modeling, which are also indicated in Fig. 2. Values inside the circular mask and outside the ellipse are not considered in the modeling.

Fig. 1 SPHERE/IRDIS polarimetric imaging observations of GSC 07396-00759. Shown from left to right are the total intensity Itot (log stretch) and the Qϕ and Uϕ images (both in linear stretch). The Qϕ and Uϕ images were convolved with a Gaussian point spread function (σ = 2 pixel) to increase the visibility in this plot. The coronagraphic mask is indicated by the green circular region in each panel. North is to the top and east to the left in each panel. It should be noted that the disk is not visible in the total intensity image (left panel). The structure appearing to extend from the NW to the SE is a point spread function artifact.

Fig. 2 Signal-to-noise map estimated from the Qϕ and noise images (without convolution) of our IRDIS DPI observations. The plot is shown with a linear stretch between [− 3σ, 3σ]. The inner and outer model boundaries are indicated by the circle and the ellipse in the image. North is to the top and east to the left.

2.3 Stellar properties

Along with the observations of GSC 07396-00759, and a thorough analysis for possible companions, Sissa et al. (2018) collected and estimated various stellar parameters, and we refer the reader to their paper for a detailed description of these parameters. Nevertheless, we revised the most important parameter, the distance to GSC 07396-00759, using new astrometric data from the Gaia satellite mission data release 2 (Gaia DR2, Gaia Collaboration 2018). In Table 2, we summarize the resulting stellar parameters, which we estimated as follows.

Employing the Virtual Observatory SED Analyzer tool³ (VOSA; Bayo et al. 2008), we first collected photometric data ranging from the ultraviolet (UV) at ~ 0.297μm (XMM-OM/UVW1, 0.38 mJy, Page et al. 2012) to the mid-infrared wavelengths, with the longest wavelength being 22 μm at this point. We then applied a BT-Settl-CIFIST model grid of theoretical spectra (Baraffe et al. 2015) to derive the stellar parameters from the constructed SED. The data indicate a small UV excess that we account for by reducing the weight of the UV contribution to the fit. Furthermore, allowing the interstellar extinction (A_v) to vary between 0 and 0.4, we find an effective temperature of T_eff = 3800 ± 100 K and a stellar luminosity of L_⋆ = 0.132L_⊙ for GSC 07396-00759. From the position in the Hertzsprung–Russell diagram (HRD) combined with theoretical isochrones and mass tracks (Baraffe et al. 2015) we estimate an upper limit for the age of ≲ 20 Myr and a stellar mass of ~ 0.62M_⊙ (but with a flat distribution between 0.6 and 0.68 M_⊙). From the dilution factor of the SED fit, combined with the distance of 71.43 pc, we estimate a stellar radius of R_⋆ ~ 0.71R_⊙.

When compared with the previous estimates by Sissa et al. (2018), we find a slightly higher effective temperature, but a significantly smaller stellar radius for GSC 07396-00759. This, however, is not entirely surprising ifwe take into account that i) the luminosity, and thus the stellar radius, depends on the distance to the star, and ii) in particular young stars that are members of young moving groups, such as the β Pic MG, appear over-luminous, and hence inflated, when compared to older (field) stars (Malo et al. 2014b). This, combined with the choice of the model grid could explain the discrepancy in the estimated stellar properties between our results and those obtained by Sissa et al. (2018).

Additionally, we analyze available light curves of GSC 07396-00759 observed with the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015), and during the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017), respectively, to measure the stellar rotation period. GSC 07396-00759 was observed with TESS (600–1000 nm), covering an observation period of about 27 days, between 2019-06-19 and 2019-07-17, whereas the ASAS-SN light curve was obtained in the V -band over a time period of 2.5 yr between 2016-03-10 and 2018-09-22. Using Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration 2018), we estimate a stellar rotation period P_rot = 11.63 ± 0.02 d from the Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) of the TESS light curve (see Fig. A.1), and P_rot = 12.06 ± 0.02 d from the ASAS-SN light curve (see Fig. A.2), respectively. The uncertainty associated with each period was estimated by resampling the light curve using the bootstrap method (with replacement), and corresponds to the 95% confidence interval of our sample estimates. We note, however, that these are only statistical uncertainties and they do not represent the totality of the error budget of the light curves. For example, the errors do not reflect properly that the TESS light curve only covers about two periods, or that the period estimates are probably affected by two additional sources (Δmag = 2) that fall on the same TESS pixel as GSC 07396-00759, given the pixel scale of 21′′ per pixel. Nevertheless, both period estimates are in good agreement with the period of P_rot = 12.05 ± 0.5 d reported by Messina et al. (2017). From the period estimated from the TESS light curve, together with the evaluated radius (R_⋆ ~ 0.71R_⊙) and a v sini = 3.0 ± 1.5 km s⁻¹ (Malo et al. 2014a) we estimate a stellar inclination of i ~ 80° ± 20° (see, e.g., Eq. (1); Justesen & Albrecht 2020), consistent with a stellar rotation that is likely co-planar with that of the disk, as proposed by Sissa et al. (2018). Very few disks have a measured stellar inclination to which the disk inclination can be compared (see, e.g., Greaves et al. 2014). Our results therefore might be helpful in future studies of the relation between the stellar inclination and the inclination of the disk.

As mentioned above, the model fit of the SED of GSC 07396-00759 suggests a small UV excess. We also found that GSC 07396-00759 has various detections in the X-ray band, reported in the fourth generation of serendipitous source catalogs (4XMM; Traulsen et al. 2020). The detected release of energy in the X-ray and UV portion of the stellar spectrum suggests that GSC 07396-00759 is subject to strong coronal emission, and possibly intense and frequent flares, common for young and low-mass stars, such as GSC 07396-00759. These stellar flares are caused by the reconnection of magnetic field loops on the surface of the star (Mondrik et al. 2019). How frequently these stellar flares occur also depends on how magnetically active the star is. For example, AU Mic, which is similar in luminosity and age to GSC 07396-00759 (both are members of the β Pic MG), is known to frequently present X-ray and extreme-UV flares at a rate of about 0.9 flares per hour, during which the extreme-UV and X-ray luminosities increase by a typical factor 10 (see, e.g., Augereau & Beust 2006, and references therein). In fact, recent studies of the flare activity of young (≤100 Myr) K and M stars in the Upper Sco region using Kepler K2 data found that early- and late-type M stars might have 10 000 to 80 000 times as many high-energy flares, so-called super-flares with E ≥ 5 × 10³⁴ erg, than solar-like stars (Guenther et al. 2019). Such an enormous activity will certainly affect the direct environment around the stars, especially potential disks and/or planets that orbit the star.

Although we do not detect flares in the available light curves of GSC 07396-00759, such an increased flare rate might also be the case for GSC 07396-00759, because not only is the star presumably young (≈ 20 Myr), but it is likely also magnetically active, as indicated by the detection of the Hα and Hβ lines in emission (Sissa et al. 2018), and the moderate rotation period of P_rot ≈ 12 days. Thus, even if GSC 07396-00759 is currently in a quiescent state, the present coronal activity of the star and its effect on the disk may not be negligible, given the observed variations in the X-ray flux and the apparent UV excess detected in GSC 07396-00759.

3 Analysis

Besides gravity, the orbital parameters of small micron-sized particles, observed in the visible and near-infrared using scattered and linear-polarized light, are mainly affected by radiation pressure and stellar winds. In particular the latter can dominate the effects on circumstellar grains in debris disks around late-type, low-mass stars as first pointed out by Plavchan et al. (2005). Hence, in this study we try to constrain the dust properties, as well as the distribution of dust grains in the debris disk around GSC 07396-00759, taking into account not only radiation pressure but also the effect from the stellar wind of the central star.

3.1 Code description

The code used in this study is based on the same code as was described in Olofsson et al. (2019). However, since we have made some changes to the code, in the following, we give a short summary of the code description and the modifications we made to also accountfor the effects of stellar wind pressure on the dust grains.

The observed dust is released from a belt of planetesimal-sized parent bodies (their numbers, sizes and masses do not need to be specified), which is defined by six parameters: a reference radius r₀, which is the semimajor axis of the orbit, eccentricity e, inclination i, argument of periapsis ω, the position angle on the sky φ, measured positive from north to east, and the width of the ring δ_r. All the dust grains originate from those parent bodies that, a priori, all share the same e and ω, and have a radial distribution that follows a normal distribution centered at a radius r₀ with a standard deviation of δ_r. The dust grain-size distribution, used in the code, follows a differential power law, (1)

where s is the grain size, and p = 3.5 is the slope index of a particle distribution derived from an idealized collisional cascade, following Dohnanyi (1969). This distribution is divided in n_g intervals (equidistant in logarithmic space between s_min and s_max, respectively), and the number density of grains in each bin is computed using Eq. (2) of Dullemond & Dominik (2008).

For each grain size the code then computes the dimensionless parameter β, which, in previous versions of the code, was defined as the ratio between the radiation pressure and the gravitational forces (Burns et al. 1979). If this ratio exceeds 0.5 for an initially circular orbit, the dust particle gets pushed into a hyperbolic orbit and leaves the system (Krivov et al. 2006). Hence, for a given grain composition and porosity, this β parameter cutoff can be used as a very close approximation of the effective radius below which dust particles will leave the system (Arnold et al. 2019). Following the investigation of AU Mic (Augereau & Beust 2006; Sezestre et al. 2017) we modified the code to also include effect of the stellar wind pressure force on the dust grains. The net pressure force acting on a grain is then defined by β = β_RP + β_SW (Sezestre et al. 2017), where β_RP is the ratio between the radiation pressure forces and the gravitational forces, and β_SW is the ratio between the wind pressure and gravitational forces. The individual contributions of the two pressure forces to β can be estimated via (2)

where Ṁ_⋆ is the stellar mass-loss rate, V_SW is the stellar wind speed, C_D the dimensionless free molecular drag coefficient, which we take equal to 2 (see, e.g., Augereau & Beust 2006, and references therein), G the gravitational constant, M_⋆ the mass of the star, and ρ the volumetric mass density of the dust, respectively, and (3)

where L_⋆ the stellar luminosity, Q_RP the dimensionless radiation pressure efficiency, which along the wavelength also depends on the dust grain size and composition, F_λ the stellar flux, and λ the wavelength. As can be seen from the equations, in general and for sufficiently large dust grains β varies with s⁻¹. It should be noted however that for smaller grain sizes, the relationship between β and s becomes more complex and increasingly dependent on the grain composition, the stellar mass-loss rate Ṁ_⋆, and stellar wind speed V_SW (Sezestre et al. 2017).

For a given β the code draws 3000 samples from an uninformative uniform prior distribution for the mean anomaly, to decide where the dust grain is located upon its release. For each realization, the code then calculates the “updated” orbital parameter (a_n, e_n, and ω_n) using Eq. (2) in Wyatt et al. (1999); Wyatt (2006); Lee & Chiang (2016). In case the updated eccentricity e_n is larger or equal to zero and strictly smaller than unity (to avoid hyperbolic orbits), the resulting orbits are populated with 300 dust particles, uniformly distributed in mean anomaly. The vertical distribution of the disk is accounted for by drawing from a normal distribution with a standard deviation h = 0.04 × r, following Thébault (2009). This allows us to account for column density effects, as explained in further detail in Olofsson et al. (2020). The (x, y, z) positions of each particle are registered, and depending on inclination and position angle of the disk, the corresponding closest pixel is determined, thus producing number density maps for each value of β.

Furthermore, and as discussed in Olofsson et al. (2019), when computing the number density maps for each β value, the contributionof each particle is also multiplied by a correction factor (Strubbe & Chiang 2006; Lee & Chiang 2016), which is roughly proportionalto their total orbital period divided by the time spent within the birth ring. This correcting “enhancement” factor to the high-β grain number density allows us to at least to the first order account for the fact 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, hence enhancing significantly their collisional lifetimes, and thus their number density (Strubbe & Chiang 2006; Thébault & Wu 2008).

Once each of the 3000 × 300 particle has been launched the code computes the scattering angle between the central star and the observer for each particle in the image. The polarized light images per grain size bin are then computed by multiplying the estimated number density of each pixel by S₁₂ × πs² × Q_sca∕(4πr²), where r is the distance from the particle to the star, Q_sca the scattering efficiency, and S₁₂ the polarized phase function. The final image is the collapse of all individual images for each grain size bin, weighted by the grain size distribution n(s).

3.2 Modeling strategy

3.2.1 Caveats and modeling approach

Most likely due to a degeneracy between the minimum grain size and the porosity of the dust particles or the high inclination of the disk around GSC 07396-00759, our primary modeling attempts were not converging on a unique solution (see Sect. 4 for a detailed discussion). To simplify the problem, we therefore used a different, two-step approach in which we first focused on the morphology of the disk and then on the dust properties. The procedure can be summarized as follows.

We first try to alleviate the influence of some of the dust properties by replacing the polarized phase function in our disk model with an analytical form, the Henyey–Greenstein (HG) approximation (Henyey & Greenstein 1941) described in Sect. 3.2.2. This form of the phase function is parameterized for all grain sizes. Thus, the absolute value of β for each grain size becomes less relevant, and the spatial distribution of the dust particles is mostly determined by the global shape of the β function (β(s) ∝ s⁻¹). By sampling the entire range of β values below 0.5, we can therefore estimate the spatial distribution of the dust particles in the disk, without having to pay too much attention to their actual sizes and properties.

Outside of the birth ring there is also an overabundance of small dust grains, mostly due to the effect of radiation pressure, as discussed for example in Thébault & Wu (2008). While the grain size distribution is expected to follow the “classical" distribution with an exponent p = 3.5 (see above) in the birth ring, outside of the birth ring the smallest dust grains have an increased collisional lifetime as they are set on highly eccentric orbits and, thus, spend most of their time near their apoapsis (i.e., they survive longer and hence contribute even more to the flux observed in scattered light). In fact, the deviation of the size distribution from -3.5 can actually be quite significant (see, e.g., Fig. 3 of Strubbe & Chiang 2006), resulting in the small dust grains dominating outside of the birth ring. Therefore, the phase function we retrieve is the one for the dominant/most representative grain size, which we assume to be close to the blow out size of the dust particle.

With the radial distribution and the phase function established by our disk model, we can then, in the second step, try to conversely infer other properties of the dust particles. We achieve this by comparing the inferred analytical polarized phase function with a model-dependent polarized phase function computed by other means such as effective medium theory (EMT), the discrete dipole approximation (DDA; Purcell & Pennypacker 1973; Draine 1988), or a distribution of hollow spheres (DHS; see Min et al. 2005). Depending on the grain composition and porosity we can then try to estimate individual parameters, such as the effective radius (s_min), below which particles will not stay in the system, or the stellar mass-loss rate (Ṁ_⋆).

3.2.2 Modeling the disk

The free parameters of our disk model are the reference radius of the radial distribution r₀, its standard deviation δ_r, the inclination i, and the position angle ϕ. The high inclination of the disk around GSC 07396-00759, however, makes it hard to reasonably constrain the eccentricity. We therefore assume a circular disk and fix the eccentricity to e = 0, making the argument of periapsis ω irrelevant for the modeling. For the stellar parameters we use a bolometric luminosity of 0.13 L_⋆, a stellar mass of 0.62 M_⊙, and a distance of 71.43 pc.

We sampled the grain size distribution so that the entire range of β values from 0.5 down to 10⁻⁵ is probed, that is to say, the dust particles remain gravitationally bound over the entire grain size distribution between s_min and s_max. To determine the polarized phase function S₁₂ as a function of the scattering angle, we use a modification of the analytical HG approximation (Henyey & Greenstein 1941), which is given by (4)

where θ is the scattering angle and g the anisotropic scattering factor (−1 ≤ g ≤ 1), which governs the scattering efficiency as a function of the scattering angle θ. The first factor in Eq. (4 table 4) is the HG function that describes the scattered flux produced by the photons hitting and interacting with the dust particles, while the second term in Eq. (4) takes into account the angle dependence of the linear polarization produced by the particle scattering, using the Rayleigh scattering function as a simple approximation. Thus, for g = 0 (isotropic scattering) the maximum of scattered polarized flux occurs at θ = 90°, while for g > 0 the maximum is shifted to smaller scattering angles resulting in an asymmetry in the amount of polarized light received from the front and backsides of the disk (Engler et al. 2017).

The best solution to our model with five free parameters (r₀, δ_r, i, ϕ, and g) is determined using an affine invariant ensemble sample Monte Carlo Markov chain (MCMC; emcee package, Foreman-Mackey et al. 2013). For the initial conditions of the model we use the disk parameters (r₀, i, and ϕ) reported by Sissa et al. (2018). We then draw random samples from uniform priors we also report in Table 3, using an MCMC composed of 30 walkers, a burn-in phase of 1000 models, followed by the actual modeling using chains of 10 000 models for each walker. To speed up the modeling process, we exclude all image data points located outside of an elliptical mask defined by a semimajor and semiminor axis of 2.3′′ and 0.6′′, and a position angle of 150°, respectively. We also place a numerical mask with a radius of 0.3′′ on the center to reduce the influence of the increased noise around the central star on the modeling (see Fig. 2). The goodness of fit for each model is estimated as the sum of the squared residuals divided by the uncertainties. The results of the disk modeling are presented in Sect. 4.

Fig. 3 Observation and model images for GSC 07396-00759. Shown from left to right are the measured Qϕ, best-fit model, and residual image obtained from our disk model using a linear scaling for each panel. The regions outside of the ellipse and within the circle (dashed green lines) are excluded from the χ2 calculation. The coronagraphic mask is indicated by the (shaded) circular region in each panel. North is to the top and east to left in each panel.

3.2.3 Polarized phase function

In the second step of our analysis we use the parameterized polarized phase function obtained in the previous step to try to infer further dust properties by comparing it to a model-dependent phase function. In this study, we use the OpacityTool, a dedicated Fortran package⁴ of the DIANA project (see Woitke et al. 2016; Toon & Ackerman 1981), to compute dust opacities for the purpose of disk modeling. The package uses a DHS (see Min et al. 2005) to compute absorption and scattering properties, as well as six different elements of the scattering matrix, including S₁₁, the scattering phase function, and S₁₂, the polarized phase function, respectively. We used the same grain size distribution as in the previous step to calculate polarized phase functions.

Preliminary tests, however, suggested that a size distribution between s_min and s_max (s_max ≫ s_min) could not match the phase function derived from the modeling of the disk, given the limitations of the data and the observed degeneraciesin the model. Therefore, we chose to only compare our results to models with one representative grain size to find the typical grain size that can best explain the inferred polarized phase function. For each grain size we calculate polarized phase functions varying the porosities from 0.01 to 0.99 in steps of 0.01, and the grain size (s) between 0.01 μm (s_min) and 1 mm (s_max), respectively. The OpacityTool calculates the opacities and other properties using a mixture of 85% amorphous laboratory silicates (Dorschner et al. 1995, Mg_0.7Fe_0.3SiO₃) with 15% amorphous carbon (Zubko et al. 1996, BE-sample). The effective refractory index of the porous material is calculated by applying the Bruggeman (1935) mixing rule. The maximum hollow volume ratio is set to f_max = 0.8. This filling factor f_max, computationally, represents the maximum volume fraction occupied by the central void in the hollow sphere, while in practice, this parameter represents the more general amount of deviation from a perfect homogeneous sphere (Min et al. 2016).

Using the grain size and porosity as the only free parameters, we then compute a grid of models to estimate whether we can reproduce the parameterized phase function. As a consequence of the high inclination of the disk and the employed inner mask not all scattering angles between 0 and 180 degrees can be sampled. To account for this, we use the best-fit parameters determined for the disk and the applied masks to estimate range of valid scattering angle. The result of this exercise is shown in Fig. B.1, yielding observable scattering angle from 9.14° to 170.86° (θ_min,max = 90° ± 80.86°). We also note at this point that due to the method employed, the uncertainties in our disk model may be underestimated. We therefore utilize the 3σ error bounds on the value of g, rather than the 1σ values.

4 Results

4.1 Geometric properties of the disk

We first used the model described in the previous section to reproduce the observed disk around GSC 07396-00759, and to determine the most probable values for each of the five free parameters (r₀, δ_r, i, ϕ, and g). To assess the convergence and stability of the MCMC solution, we estimated the maximum autocorrelation length among all parameters (i.e., the average autocorrelation time) after each iteration. The fitting is considered converged when the number of iterations is larger than 100 times the average autocorrelation time and its changes, between subsequent iterations, are less than 1%. In Appendix C.1, we show the evolution of the autocorrelation time as a function of the iteration step over the course of the modeling. At the end of the modeling, the average autocorrelation time was 79 steps and the mean acceptance fraction (Gelman & Rubin 1992) for our best fitting model was 0.47. The best-fit model, along with the residuals and the observations, are presented with the same linear scaling in the center, right and left panels of Fig. 3, respectively.The most probable parameters and their uncertainties are summarized in Table 3, while the probability density functions, plotted with corner package (Foreman-Mackey 2016), are shown in Fig. 4.

We found that our observations are best described by an extended disk with a dust distribution centered at a radius r₀ ≈ 107 ± 2 au, with a standard deviation of the radial distribution of δ_r ≈ 27 ± 1 au, which is highly inclined at an inclination of i ≈ 84.3° ± 0.3, and a position angle ϕ ≈ 148.7° ± 0.1. For the analytical polarized phase function we found a coefficient g = 0.60 ± 0.03.

The uncertainties for the MCMC result parameters are estimated from the 0.16 and 0.84 quartiles using the corner package, and are also shown along with the projected posterior distributions in Fig. 4. These uncertainties, however, are most likely underestimated and should be taken with caution. This might be because we computed the goodness of the fit using a noise map derived from the standard deviation in concentric annuli in the U_ϕ image, a strategy commonly employed in direct imaging studies. We therefore may be underestimating the true uncertainties, resulting in larger χ² values. This in turn might also have led to a narrower probability distribution as the Monte Carlo algorithm samples a smaller range of parameter values.

Fig. 4 Projected posterior distribution and density plots of the free parameters used in the disk modeling. The plot additionally shows the 50%, 16%, and 84% quartiles (vertical dashed lines), representing the distribution median and the 1 σ uncertainties (lower and upper bound), respectively. The 1D histograms represent the probability distributions of each parameter marginalized over the other free parameter. The contour maps represent the central 68.3%, 95.5%, and 99.73% of the 2D probability distributions of different combinations of parameters, marginalized over one another.

4.2 Dust properties

Next we applied a grid of models, computed with the OpacityTool, as described in Sect. 3.2.3, to estimate whether we can reproduce the inferred polarized phase function (S_{12, HG}) considering a DHS. The resulting χ² map is shown in Fig. 5. Also shown are the contours that enclose the regions where the χ² is smaller than the lower 5% (dashed line) and 15% (solid line) quartiles. We found that the analytical form of the polarized phase function (S_{12, HG}) is best reproduced by small micron-sized dust grains with s ~ 0.34μm. Due to a degeneracy between the grain size and the dust porosity, however, we cannot readily constrain the porosity for this grain size (see discussion below).

In Fig. 6, we present the parameterized polarized phase function (S_{12, HG}) as a function of the scattering angle θ from our disk modeling for (black line and gray area), along with the best-fit model of the polarized phase function (S_{12, DHS}, red line). Also shown is the 1 σ range of models (red area, enclosing values between the 16% and 84% quartiles of the likelihood distribution), estimated from the phase functions for s = 0.32–0.36 μm, and the full range of modeled porosities (0.01–0.99). As can be seen, and despite this being the best solution to our model approach, there are notable differences between both curves, but in particular for scattering angles of less than 90°. The χ² map, furthermore, shows “low” χ² values not only around 0.34 μm, but also for particles of 1 μm with porosities around 0.8 as well as for larger particles (s > 100μm) with porosities of less than 0.6. The corresponding range of phase functions, determined similar to the best-fit results, are shown as well in Fig. 6. These phase functions clearly provide a worse fit to our approximation (shown in black) then the best solution using a DHS model (shown in red).

The presence of multiple solutions, on the other hand, could indicate additional hidden degeneracies beyond the one on the minimum size and porosity mentioned, restricted only to the second step of our alternative two step approach. For instance, we use a parameterized approximation (S_{12, HG}) of the real phase function (S₁₂), which effectively limits the examinable parameter space. We then fit the obtained curve with the phase function computed for a DHS model.However, it is questionable whether we can reasonably reproduce the analytical HG approximation of the phase function with the phase function from a DHS model.

To address this question we created a map similar to the χ² map shown in Fig. 5, where each point of the grid represents the best g value fittingthe DHS model for a given set of grain size and porosity. The resulting map of g values (hereafter called ) is shown in Fig. 7. Also shown are the regions where the results overlap with the range of g values from the HG approximation. As can be seen, for the smaller grains (0.3–1 μm) the value estimates are mostly consistent with the range of g values using the HG approximation, while the g values for the larger particles appear systematically higher. The data also indicate that there are possibly isolated solutions for particles with sizes between 5 μm and 10 μm. We also summarize these results in the lower part of Table 3.

Fig. 5 2D χ2 map constructed from the computed χ2 at each point (model) of the parameter grid. The free parameters of the model are the porosity and the grain size. The contours enclose the regions where the χ2 is smaller than the 5% (dashed line) and 15% (solid line) quartiles. To increase the contrast, the upper limit of the χ2 color map isset to the 50% quartiles.

Fig. 6 Polarized phase function (S12, HG) as a functionof scattering angle θ. Shown are the parameterized polarized phase functions from our best-fit disk model for g = 0.60 ± 0.09 (black line and gray area) and the range (red area, enclosing values between the 16% and 84% quartiles) of the best-fit model (s = 0.34μm, porosity = 0.01–0.99) of this phase function, considering a DHS. Additionally shown are the polarized phase functions resulting from the “low”-χ2 regions indicated in Fig. 5. The phase functions were sampled around grain sizes and porosities of about 1 μm and 0.80 (blue area) and grains larger than 100 μm and porosities below 0.6 (green area), respectively. The range of possible scattering angles (θ = 9.14° to θ = 170.86°) observable with the used mask is indicated by the blue shaded area.

Fig. 7 2D map of the anisotropic scattering factor computed at each point of the grain-size, porosity parameter grid with the DHS model. The contours for (blue lines) and (red lines) indicate the range of scattering factors from our best-fit disk model with the HG approximation.

Fig. 8 Disk around GSC 07396-00759, as seen with SPHERE/IRDIS in polarized light (top-left panel) and total intensity (top-right panel; Sissa et al. 2018) in the H-band. The bottom panels show the best-fit disk model found in this work (bottom left) and found by Sissa et al. (2018) in the same scaling as the according observations. The dash-dotted line is parallel to the disk semimajor axis and crosses the central star.The coronagraphic mask is indicated by the (shaded) circular region in each panel.

5 Discussion

5.1 Disk properties

Our results from the modeling of the dust distribution around GSC 07396-00759 indicate a radius of ~ 107 au for the radial distribution of the grain number density of the debris disk, and as shown in the residuals image in the right panel of Fig. 3, most of the signal coming from the disk has been successfully removed.

The modeling results are globally consistent with the results reported by Sissa et al. (2018). In particular, we found that the estimated full width at half maximum (FWHM) of the disk of au (δ_r ≈ 27 au) as well as the inclination (i ≈ 84°) and position angle (θ ≈ 149°) for the disk around GSC 07396-00759 are in good agreement with the results obtained by Sissa et al. (2018, i ≈ 83°, θ ≈ 149°, and FWHM ≈ 56 au), although their forward modeling of the angular differential imaging (ADI) observations of GSC 07396-00759 yielded a smaller disk radius (r₀ ≈ 70 au) than our model of the polarized light observations (r₀ ≈ 107 au).

Our disk radius estimate is also slightly larger than the radius obtained by Sissa et al. (2018) using a geometric approach, where the authors measured the position of the maximum brightness (disk spine) with respect to the apparent semimajor axis (i.e., the apparent disk radius obtained via forward modeling) as a function of the separation from the star after the removal of background stars in the five-mode principle component analysis IRDIS image (see Fig. 2 of Sissa et al. 2018). Under the assumption that the disk is not flared, Sissa et al. (2018) found that a ring of material with radius r = 1.34′′ (~ 98 au) and inclination i = 84.5° can properly describe the observed spine up to 1.2′′ (~ 89 au) in their ADI images. However, as also pointed out, the presence of at least three very bright stars may have altered the light distribution of the disk, and thus, the results have to be taken with caution as well.

It is surprising, nevertheless, that we found r₀ to be ≈107 au (≈ 1.5′′) in our model, while we only detected uniform emission up to about ≈ 70 − 80 au (1–1.1′′) in the image. This suggests that according to our model, we do not fully recover the major axis of the disk with our observations, because we also found that the maximum of the polarized phase function (see Fig. 6) is shifted to smaller scattering angle (θ ≈ 55°), while the phase function at the same time is relatively low at scattering angles of 90°.

We suggest that this is due to a degeneracy between the reference radius of the density distribution r₀ and the value of the anisotropic scattering factor g for the phase function. We emphasize, however, that this degeneracy is solely the case for highly inclined disks, because significant azimuthal information is lost in this particular, nearly edge-on configuration (see, e.g., Fig. 8 of Engler et al. 2017).

Our value of g = 0.60 therefore might be too large. Nevertheless, we can exclude very small g values, because, if g were close to zero, the major axis of the disk would become much brighter than the minor axis, both due to the phase function peaking closer to 90 degrees and the column density that is larger at the major axis of the disk (see, e.g., Olofsson et al. 2020, for discussion). This, however, is not what we observe in the images, because we see a more or less uniform distribution of polarized light along the disk, as opposed to two distinct lobes, seen for example in the modeling of the disk around HD 61005 (Olofsson et al. 2016). The observed degeneracy might also play a significant role for ADI images, since forward-scattering is even more pronounced in ADI observations, and therefore, the results of Sissa et al. (2018) may also be impacted bythis, as indicated by the side-by-side comparison of the two observations shown in Fig. 8.

To explore the possibility that our radius for the disk around GSC 07396-00759 is overestimated, we reran our grid of models with the reference radius fixed to r₀ = 69.9 au, obtained by Sissa et al. (2018) using forward modeling, to test whether we can find a solution that fits our observations. The resulting best-fit model in this scenario is shown in Fig. D.1. In particular, we found a standard deviation of the radial distribution of δ_r ≈ 15.3 ± 0.2 au, an inclination of i ≈ 78.5° ± 0.3 and a position angle ϕ ≈ 149.2° ± 0.1. For the analytical polarized phase function we found a coefficient g = 0.33. Although the model converged on a solution, the results clearly show that we under-estimate the reference radius, as in particular the outer regions are not well reproduced, when using the results from the forward modeling of the ADI observation by Sissa et al. (2018). On the other hand, the reference radius we inferred is in fact quite close to the radius of the disk of r₀ ≊ 98 au obtained by Sissa et al. (2018) from their geometric measurements. We thus suggest that the estimated reference radius (r₀ ≊ 107 au), while not definitive, is probably a better estimate than the one inferred by Sissa et al. (2018) using forward modeling. However, an independent confirmation of the disk radius is required to settle the discussion. One possible way to alleviate this degeneracy and to resolve the discrepancies between the different measurements is to obtain high angular resolution observations at longer wavelengths with the Atacama Large Millimeter Array (ALMA; Cronin-Coltsmann, in prep.). These observationsare more sensitive to larger dust grains, and would therefore allow us to place stronger constraints on the shape parameters of the debris disk, and consequently allow us to more reliably estimate the shape of the phase function.

5.2 Disk radii

Another important diagnostic of debris disks alongside the phase function and the size distribution of their dust are the dimensions of the debris disks (i.e., their extent and radii). Believed to be the by-product of the planet formation process, the dust observed in debris disks is thought to be produced by the grinding down of larger bodies, the planetesimals. In general, theory would suggest that the presence of a planetesimal belt, be it caused by failed growth to planets or enhanced planetesimal formation due to complex and dynamical interaction with already formed planets, should be related to distance to the central star (Matrà et al. 2018). Regardless of the underlying processes, it is clear that the radii of debris disks contain valuable information on the formation processes of planetesimals and planets.

Several studies of the radius–luminosity (R–L) relation have been conducted over the years on this topic, such as for Herschel PACS-resolved disk radii between 70 μm and 160 μm (Pawellek et al. 2014), for planetesimal belt central radii from thermal light imaging (Matrà et al. 2018), and recently on 25 GPIES-detected debris disks (L_⋆ = 2–14 L_⊙) in scattered light by Esposito et al. (2020). These studies found only a marginal correlation between the luminosity of the star and the radius of the disk, indicating that the morphology of debris disk is likely a result of temperature-independent processes, such as shaping of planetesimal populations by planets, stirring by various mechanisms, and long-term collisional depletion (Pawellek et al. 2014).

Despite the interesting results, the studies lack stars with L < 2L_⊙ at the lower end of the mass spectrum, which may have introduced a selection bias affecting their results. With the newly available scattered light data from resolved disks around low-mass stars with L_⋆ < 2L_⊙, we therefore compiled a list of stars with resolved debris disks, and observed in scattered light from the literature, and attempt to test for a possible correlation between the size of the disk and the properties of the host star, such as its spectral type or luminosity. We chose to primarily focus on scattered light observations, because the size and radius of the observed disk depends also on the wavelength in which they are observed, as well as the strength of radiation pressure and stellar winds depending on the evolutionary state of the star.

Our sample contains a total of 46 stars (see Table E.1), including GSC 07396-00759, with spectral types ranging from B9 to M3, ages between 10 Myr and 2 Gyr, and stellar luminosities reaching from 0.1 L_⊙ to 25 L_⊙. In case there are several publications for a given target, we focused on those publications with the highest S/N observations, and a similar modeling approach to ours. We then used the published parameters that describe the radial dust density distribution of the disk, such as r₀, Δr, α_in, and α_out, to estimate the FWHM as a function of the “peak” radius of the radial distribution r₀ along the mid-plane of the disk. Out of the total sample presented, for 11 stars only the detected inner and outer extent of the disk (r_min, r_max) were given. In these cases we equate the peak radius r₀ to the mean of r_min, and r_max, and the difference r_max − r_min as width of the disk, respectively. For the other 35 stars, the estimated FWHM is based on a radial distribution that was directly derived either from the standard deviation of the assumed Gaussian, or a radial distribution R(r) that uses two power laws and is given by the expression: (5)

where r is the radial coordinate in the equatorial plane, and r_c is a critical radius that divides the ring into inner and outer regions with separate density power-law indices of α_in and α_out, respectively (Augereau et al. 1999). As also pointed out by Augereau et al. (1999), the maximum of the dust density does not occur at r_c but at a peak radius calculated via (6)

where Γ_in = α_in + γ and Γ_out = α_out + γ. For this analysis we set γ = 1, thus assuming the scale height is a constant fraction of the radius throughout the disk.

To estimate the dependence of the disk radii on the stellar luminosity, we followed the description in Matrà et al. (2018) and used a power-law model in which the belt locations R_i (in au) are linked to their host star luminosities L_⋆,i (in L_⊙) through the form . The parameter ϵ_i represents the intrinsic scatter of the distribution, and is assumed to follow a Gaussian distribution with standard deviation σ_intr = f_ΔRR_i. The best solution to this power-law model is determined using an affine invariant ensemble sample MCMC (emcee package, Foreman-Mackey et al. 2013). We then draw random samples from an uninformative uniform prior on the free parameters , α and f_ΔR, using a likelihood function described by Eq. (24) in Kelly (2007) to determine the posterior probability distributions for the parameters. To facilitate the comparison with Esposito et al. (2020), we excluded HR 7012 from the fit as an outlier too.

In Figs. 9 and 10, we present the results of this analysis. Figure 9 shows the estimated FWHM (left panel) and the FWHM divided by r₀ (right panel), respectively, as a function of the radius of the disk r₀. To point out that the sample is not homogeneous, we mark the stars where only the inner and outer radius of the disk is given with a diamond, while the results obtained from density distributions using a Gaussian or power laws to describe the radial profile are marked by circles. The stars are color-coded according to their spectral type. We additionally indicated the limits for FWHM = r₀ and FWHM = 2r₀ to guide theeye. The largest disk contained in our sample is that of HD 15745, a F2V star whose disk was discovered by Kalas et al. (2007). The star is one of 10 stars that have a disk with a FWHM larger than the radius of the respective disk, with the largest width-to-radius ratios being found for AU Mic (Krist et al. 2005), HD 15745, and β Pic (Millar-Blanchaer et al. 2015). Most of the stars have FWHM smaller than the radius of the disk, and overall we found neither an obvious correlation nor a clustering of data between the radius of the disk, its radial extent, and the spectral type of the star.

From the analysis of the radius–luminosity relationship in our sample we found au, α = 0.02 ± 0.08, and f_ΔR = 0.07 ± 0.02, taking the percentiles of the posterior probability distributions for each parameter. The results are illustrated in Fig. 10. Shown are the scattered light disk radii as a function of the stellar luminosity of our sample stars, as well as our best-fit power-law function as solid (red) line and the 1-σ confidence interval estimated from a randomly drawn sample of the corresponding probability distributions shown in Fig. E.1. The data points that are also contained in the sample analyzed by Esposito et al. (2020) are marked as diamonds. Because of the inconsistent radius estimates between this work and the work of Sissa et al. (2018), we repeated the fit using the radius estimate by Sissa et al. (2018) and show the result as solid (black) line.However, it is apparent that the derived relation is not strongly dependent on which estimate of the radius of the disk around GSC 07396-00759 is included. For comparison, we furthermore plotted the radius-luminosity power laws for the planetesimal belt central radii from thermal light imaging by Matrà et al. (2018, α = 0.19 ± 0.04, dash-dotted line), and the GPIES-detected disks in scattered light by Esposito et al. (2020, α = 0.25 ± 0.09, dashed line), respectively.

As can be seen, our result differs from the results reported by Matrà et al. (2018) and Esposito et al. (2020), both reporting evidence of a statistical relation between the radius and the luminosity (although marginal), while our result suggest that there is no correlation between the stellar luminosity L_⋆ and the scattered-light radius r₀. This would be consistent with the Pawellek et al. (2014) finding of no significant correlation between L_⋆ and r₀ detected in the Herschel PACS survey, and further support for the idea that the dimensions of debris disks are likely set or influenced by temperature-independent processes, such as collisions or dynamical interaction with a planetary perturber. However, we want to emphasize that most of the Herschel detections are marginally resolved and therefore are not directly comparable to our results. The presented sample was also derived from different observations, with varying spatial resolution, contrast performance at short angular separations, sensitivities, and differing modeling strategies. One should furthermore note that, as it is the case for GSC 07396-00759, most of the disks resolved, in particular around low-mass stars, are highly inclined, making them favorable for detection in scattered light, but also susceptible to projection effects and inaccurately estimated phase functions. The influence of stellar winds on the shape of the disk, in particular on the radial distribution of the smaller dust grains, also remains rather unclear due to the low number of resolved disks around low-mass stars (L_⋆ <1 L_⊙). It is thus also possible that some radii are over- or under-estimated and may not necessarily measure the “peak” radius of the disk (i.e., the radius of the dust-producing planetesimal belt).

Because of these potential sources of uncertainty in the scattered-light radius measurements, the lack of correlation should be taken with caution, and better constraints on the radii of resolved disk around low-mass stars using independent observation methods are required to better understand the relation between the luminosity of the star, the peak of the radial dust distribution of debris disks observed in scattered light, and its implications for the formation and evolution of the dust-producing parent belt.

Fig. 9 FWHM as a function of the center of the radial distribution r0 (left panel), and FWHM normalized by the radius of the disk as a function of r0 (right panel). Objects for which the FWHM was estimated readily from a radial Gaussian or power-law distribution are shown as circles, while objects where only the inner and outer radius were given are marked by diamonds. The spectral types of the sample stars are indicated by the color bar. The radius and width of the disk around GSC 07396-00759, estimated by Sissa et al. (2018) and in this work, are marked by the cross and the asterisk, respectively. To guide the eye, we also show the limits for FWHM = r0 and FWHM = 2r0 as a gray line in the left panel.

Fig. 10 Scattered light disk radii as a function of the stellar luminosity. Shown are the best-fit power-law function, as a solid red line, and the 1σ confidence interval drawn from the corresponding probability distributions. Data points that are also part of the sample analyzed byEsposito et al. (2020) are marked as diamonds. Similar to Esposito et al. (2020), we excluded HR 7012 (blue open circle) from the fit as an outlier. The best-fit result using the radius estimate by Sissa et al. (2018) (black star) is plotted as solid black line. Also shown for comparison are the radius-luminosity power laws for planetesimal belt central radii from thermal light imaging by Matrà et al. (2018) (dash-dotted line) and the GPIES-detected disks in scattered light by Esposito et al. (2020) (dashed line).

5.3 Dust properties

Despite the aforementioned assumptions and caveats about our modeling strategy, the best fitting model does account for most of the observed signal, the probability density functions of their parameters appear well constrained, and the estimated disk properties are compatible with the previous result by Sissa et al. (2018). We therefore used the inferred parameterized phase function to further investigate the dust properties of the disk around GSC 07396-00759, as described in Sect. 3.2.3.

The χ² map obtained from our fit of the parameterized phase function is presented in Fig. 5. We found that the parameterized polarized phase function is relatively well reproduced by small dust grains of 0.34 μm (see Fig. 6). The plot, however, also revealed a degeneracy between the grain size and the porosity, which suggests that we are most likely only tracing the small dust particles (0.34–1 μm) that are dispersed throughout the disk. Composed of sub-micron-sized monomers these particles can also make up larger aggregates. As explained in Min et al. (2016), the polarization properties of such aggregates are intimately related to the size of the individual monomers and not to the overall size of the aggregate itself. This degeneracy suggests that we cannot fully reconcile all key aspects (e.g., phase function, porosity and grain size) with a single scattering theory, a well known problem in the study of debris disks (see, e.g., the review by Hughes et al. 2018, or the studies of HR 4796 A, or HD 191089 by Milli et al. 2017, 2019; Olofsson et al. 2019, and Ren et al. 2019, respectively). Since the dust grains can be aggregates composed of smaller particles, the inferred dust particle size of s ≈ 0.34μm is most likely a lower limit for the dust grain size in the disk of GSC 07396-00759.

Nevertheless, we can use this information to try to constrain the strength of the stellar winds. To put very rough constraints on the mass-loss rate of GSC 07396-00759, we used Eqs. (2) and (3) described in Sect. 3.1 to compute the total net pressure force acting on a grain (β = β_RP + β_SW, Sezestre et al. 2017), depending on the estimated minimum grain size s_min, the blow-out size of the dust grains, for which β ≤ 0.5 (assuming the parent bodies are on circular orbit). Consequently, all the grains that are smaller are no longer bound to the star and, thus, would be removed from the system rapidly. Q_RP, the radiation pressure efficiency averaged over the stellar spectrum in Eq. (3), was calculated using the asymmetry parameter g_sca, which can also be computed with the OpacityTool, and is equal to Q_ext(λ, s) − g_sca(s) × Q_sca(λ, s).

Assuming a stellar luminosity L_⋆ = 0.13L_⊙, a mass of M_⋆ = 0.62M_⊙, a minimum blow-out size of 0.34 μm for the dust grains, we estimate that the requirement of having a sufficiently strong stellar wind is obtained for average stellar mass-loss rates ranging from Ṁ_⋆ ≈ 10Ṁ_⊙ up to Ṁ_⋆ ≈ 500Ṁ_⊙, depending on the assumed porosity of the dust particles with this size. For completeness we show in Fig. 11 the total net pressure force β and the force exerted by the radiation pressure β_RP as a function of the grain size, exemplary for an assumed minimum blow-out size of s_min = 0.34μm, and a porosity of 70%, yielding a mass-loss rate of Ṁ_⋆ ≈ 250Ṁ_⊙. Although our results are consistent with findings for similar-type stars such as AU Mic (see, e.g., Schüppler et al. 2015; Sezestre et al. 2017), with mass-loss rate between Ṁ_⋆ ≈ 50Ṁ_⊙ and Ṁ_⋆ ≈ 300Ṁ_⊙ for grain sizes ranging from 0.04 μm up to 0.34 μm, respectively,we cannot derive a definitive conclusion due to the degeneracies in our modeling approach.

Fig. 11 β as a functionof grain size. βRP and β are shown as black and blue lines, respectively. The horizontal solid gray line (β = 0.5) is the upper limit for bound trajectories assuming zero eccentricity for the parent body (see Sezestre et al. 2017 for further details). The vertical dashed line indicates the corresponding minimum required grain size (s = 0.34μm) of dust particles to not be expelled from the disk, assuming a porosity of 70%.

6 Conclusion

In this paper, we have presented high angular resolution polarimetric observations of GSC 07396-00759 obtained with SPHERE/IRDIS in the near-infrared employing a broadband H filter. We have derived the stellar parameters from the SED, reconstructed using photometric data ranging from the UV to the mid-infrared wavelengths. From the position in the HRD combined with theoretical isochrones and mass tracks (Baraffe et al. 2015), we estimated an upper limit for the age of ≲ 20 Myr, a stellar mass of ~ 0.62M_⊙, and a stellar radius of R_⋆ ~ 0.71R_⊙.

With the newly determined stellar parameters, we then characterized the disk structure and modeled the dust distribution of the disk around GSC 07396-00759 using a radiative transfer model that takes into account the effects of radiation and stellar wind pressure, which is likely more efficient in low-mass stars such as GSC 07396-00759. We found that the polarized light observation is best described by our model by an extended disk with a dust distribution centered at a radius r₀ ≈ 107 ± 2 au, with a standard deviation of the radial distribution of δ_r ≈ 27 ± 1 au, which is highly inclined at an angle i ≈ 84.3° ± 0.3 and a position angle ϕ ≈ 148.7° ± 0.1. The derived scattering asymmetry parameter is g = 0.60 ± 0.03 and was estimated using a modified HG function that accounts for the angle dependence of the fractional polarization, produced by the particle scattering, by adopting the Rayleigh scattering function as an approximation.

We were not able to fully recover the major axis of the disk, because of its nearly edge-on configuration that we observe in the almost uniform emission along the disk and because of its measure in the polarized phase function. The anisotropic scattering factor g and consequently the radius of the disk may therefore be overestimated, although we argue that the observations are reasonably well reproduced by our model and that the geometric structure of the disk seen in polarized light is consistent with total intensity images taken with SPHERE-IRDIS by Sissa et al. (2018).

We further studied the dust properties, such as the shape and size of the dust particles, using a polarized phase function model that we calculated with the OpacityTool and assuming a DHS with a single grain size to fit the inferred parameterized polarized phase function. We found that the polarized phase function is reasonably well reproduced by small micron-sized dust grains with s ≈ 0.34μm. However, the modeling results also imply that the DPI observations may only trace the small monomers that are part of larger aggregates (Min et al. 2016), if they are present, and therefore the deduced dust grain size may only be a lower limit on the size of the particles dispersed throughout the disk. The results, therefore, suggest that we cannot reconcile all key aspects of the disk using a single scattering theory such as the DHS theory to explain, for example, the shape of the phase function or its dependence on the dust grain size and porosity. This is a longstanding problem in the analysis of debris disks, and further observations are required to solve this issue. Such observationscould be obtained, for example, using SPHERE/IRDIS with the star-hopping mode, which would allow both the polarized and scattered light phase functions – and therefore the degree of polarization, which contains crucial information on the shape of the phase function – to be retrieved, but they are typically challenging to measure.

Nonetheless, the observed extent of the disk around GSC 07396-00759 is comparable to the width of other disks around low-mass stars such as AU Mic or GJ 581 and appears less well constrained than the dust belts resolved in scattered light around higher-mass stars such as HR 4796 A. However, the comparison of the width as a function of the radius of the disk for 46 stars of spectral type A to M and with resolved disks in scattered light to a first approximation does not show any significant correlation between the extent of the disk, its radius, and the spectral type of the star. As a second diagnostic we also used the estimated radii and analyzed their dependence on the luminosity of the host star. From our sample of 44 disk-bearing stars with luminosities ranging from 0.1 L_⊙ to 25 L_⊙, we found that the R −L_⋆ relation is best explained by a power-law function of the form . This is consistent with no statistically significant correlation between belt radii, where most of the observed dust is released, and host star luminosities and lends further support to the idea that the dimensions of debris disks are likely set or influenced by temperature-independent processes, such as collisions or dynamical interaction with a planetary perturber. However, the peak distribution of the dust particle, and by extension the location of the planetesimal belt from which they originate, observed in scattered light may be overestimated due to projection effects, inaccuracies in the estimated phase function, or, in case of young low-mass stars such as GSC 07396-00759, the influence of the stellar activity, especially stellar winds. Thus, further studies are needed to better understand how and why the planetesimal belts arise, in particular around low-mass stars where the distribution of small dust particles may be significantly affected by stellar winds present in young, low-mass stars.

Overall, we conclude that if the disk around GSC 07396-00759 is indeed this extended (δ_r ≈ 30 au), and if the disk is dominated by small dust grains of 0.34 μm, then the stellar winds could be as strong as 500 times the solar mass-loss rate in GSC 07396-00759 and, thus, could play a dominant role in the transport of particles into the outer disk that would otherwise remain closer to their parent bodies. We note, however, that the mass-loss rate estimates should be taken with caution because of the observed degeneracies in our models. Furthermore, our mass-loss rate is primarily an averaged one over time because the considerations above do not take into account that β_SW for a given dust grain size may be time-variable because of fluctuations in the stellar flux caused by intense and frequent flares. This is common for young and low-mass stars such as GSC 07396-00759, although the light curves of GSC 07396-00759 obtained with TESS and during the ASAS-SN survey do not show signs of stellar flares. Nevertheless, the coronal activity of the star suggested by the observed variations in the UV and X-ray may not be negligible and, thus, should be taken into account in future studies in order to better understand the influence of the different pressure forces acting on the dust grains that are responsible for the observed morphology of the disk.

Acknowledgements

We are extremely grateful to the anonymous referee for carefully reading the manuscript, and providing very helpful comments to improve the quality of our manuscript. C. Adam, J.O., A.B., and M.M. acknowledge support by ANID, – Millennium Science Initiative Program – NCN19_171. J.O. acknowledges financial support from Fondecyt (grant 1180395). A.B. acknowledges financial support from Fondecyt (grant 1190748) and the Universidad de Valparaíso. A.Z. acknowledges support from the FONDECYT Iniciación en investigación project number 11 190 837. T.H. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of VOSA, developed under the Spanish Virtual Observatory project supported by the Spanish MINECO through grant AyA2017-84089. VOSA has been partially updated by using funding from the European Union’s Horizon 2020 Research and Innovation Programme, under Grant Agreement no 776403 (EXOPLANETS-A). This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by the NASA Explorer Program. This research made use of Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration 2018).

Appendix A Light curves

To measure the stellar rotation period of GSC 07396-00759 we analyzed available light curves of GSC 07396-00759 observed with TESS (Ricker et al. 2015) and during the ASAS-SN (Shappee et al. 2014; Kochanek et al. 2017), respectively. GSC 07396-00759 was observed with TESS (600–1000 nm), covering an observation period of about 27 days, between June 19, 2019, and July 17, 2019, and shown in Fig. A.1. The ASAS-SN light curve, on the other hand, was obtained in the V band over a time period of 2.5 years between March 10, 2016, and September 22, 2018, and is presented in Fig. A.2. The period estimate was obtained using Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration 2018). After each light curve was preprocessed (i.e., outlier removal and normalization was applied), we determined the stellar rotation period P_rot from the corresponding Lomb-Scargle periodogram (Lomb 1976; Scargle 1982). The uncertainty associated with each period was evaluated by resampling the light curve using the bootstrap method (with replacement), and corresponds to the 95% confidence interval of our sample estimates. It should be noted, however, that these are only statistical uncertainties and they do not represent the totality of the error budget of the light curves. From the TESS light curve of GSC 07396-00759, we estimate a stellar rotation period P_rot = 11.63 ± 0.02 d, whereas the ASAS-SN light curve yielded a period P_rot = 12.06 ± 0.02 d, respectively.

Fig. A.1 TESS light curve and periodogram of GSC 07396-00759. Left panel: Normalized TESS light curve of GSC 07396-00759, observed between June 19, 2019, and July 17, 2019. Right panel: Lomb-Scargle periodogram computed from these data, with an inset showing the light curve folded over the detected stellar rotation period Prot = 11.63 d.

Fig. A.2 ASAS-SN light curve and periodogram of GSC 07396-00759. Left panel: Normalized ASAS-SN light curve of GSC 07396-00759, observed between March 10, 2016, and September 22, 2018. Right panel: Lomb-Scargle periodogram computed from these data, with an inset showing the light curve folded over the detected stellar rotation period Prot = 12.06 d.

Appendix B Scattering angle

As mentioned in Sect. 3.2.3, we computed a grid of models to estimate whether we can reproduce the parameterized phase function using the grain size and porosity as the only free parameters. However, as a consequenceof the high inclination of the disk and the employed inner mask, not all scattering angles between 0 and 180 degrees could be sampled. To account for this, we used the best-fit parameters determined for the disk, namely r₀, δ_r, i, and ϕ, as well as the parameter of the employed masks to estimate range of valid scattering angle. In Fig. B.1 we show the results of this exercise, which yielded observable scattering angle from 9.14° to 170.86° (θ_min,max = 90° ± 80.86°).

Fig. B.1 Distribution of possible scattering angle as a function of the model geometry obtained from our best-fit model of the disk. The inner and outer model boundaries are indicated by the circle and the ellipse in the image. Values within the circular mask or outside the ellipse are not considered in the modeling. North is to the top and east to the left.

Appendix C Model convergence

The convergence and stability of our MCMC solution was assessed by estimating the maximum autocorrelation length among all parameters (i.e., the average autocorrelation time ) after each iteration. The fitting was considered converged when the number of iterations is larger than 100 times the average autocorrelation time and its changes, between subsequent iterations, are less than 1%. The evolution of the autocorrelation time as a function of the iterationstep over the course of the modeling is presented in Fig. C.1. At the end of the modeling, the average autocorrelation time was 79 steps and the mean acceptance fraction (Gelman & Rubin 1992) for our best fitting model was 0.47.

Fig. C.1 Autocorrelation time as a function of MCMC iteration steps. The estimated average autocorrelation time is shown in blue (solid line). Also shown are the (red), Niter∕50 (blue), and Niter∕100 (black) lines to indicate different possible levels of acceptance for the model parameter estimates. For this work, the model is consideredconverged if and the change in consecutive estimated autocorrelation times τ is less than 1%.

Appendix D Model comparison

To explore the possibility that our radius for the disk around GSC 07396-00759 is overestimated, we reran our grid of models with the reference radius fixed to r₀ = 69.9 au, obtained by Sissa et al. (2018) using forward modeling, to test whether we can find a solution that fits our observations. The resulting best-fit model in this scenario is shown in Fig. D.1. In particular, we found a standard deviation of the radial distribution of δ_r ≈ 15.3 ± 0.2 au, an inclination of i ≈ 78.5° ± 0.3 and a position angle ϕ ≈149.2° ± 0.1. For the analytical polarized phase function we found a coefficient g = 0.33.

Fig. D.1 Observation and model images for GSC 07396-00759. Shown from left to right are the measured Qϕ, best-fit model, and residual image obtained from our disk model using a linear scaling for each panel. The regions outside of the ellipse and within the circle (dashed green lines) are excluded from the χ2 calculation. The coronagraphic mask is indicated by the (shaded) circular region in each panel. North is to the top and east to left in each panel.

Appendix E Radius-luminosity relation

We estimated the dependence of the disk radii on the stellar luminosity following the description in Matrà et al. (2018). In particular we used a power-law model where the belt locations R_i (in au) are linked to their host star luminosities L_⋆,i (in L_⊙) through the form . The parameter ϵ_i represents the intrinsic scatter of the distribution, and is assumed to follow a Gaussian distribution with standard deviation σ_intr = f_ΔRR_i. The model was then applied to our sample, which contains a total of 46 stars, including GSC 07396-00759, with spectral types ranging from B9 to M3, ages between 10 Myr and 2 Gyr, and stellar luminosities reaching from 0.1 L_⊙ to 25 L_⊙, presented inTable E.1. The best solution to the power-law model, illustrated in Fig. 10, is determined using an affine invariant ensemble sample MCMC (emcee package, Foreman-Mackey et al. 2013). The 1σ confidence interval was estimated from a randomly drawn sample of the corresponding posterior probability distributions shown in Fig. E.1. The 1D histograms presented in this figure represent the probability distributions of each parameter marginalized over the other two, while the contour maps represent the central 68.3%, 95.5%, and 99.73% of the 2D probability distributions of different pairs of parameters, marginalized over the third.

Fig. E.1 Projected posterior distribution and density plots of the free parameter (slope α, intercept , and fractional intrinsic scatter fΔR) of the power law fitted to the sample data points. The plot additionally shows the 50 %, 16 %, and 84 % quartiles (vertical dashed lines), representing the distribution median and the 1 σ uncertainties (lower and upper bound), respectively. The 1D histograms represent the probability distributions of each parameter marginalized over the other two. The contour maps represent the central 68.3%, 95.5%, and 99.73% of the 2D probability distributions of different pairs of parameters, marginalized over the third.

References

All Tables

All Figures

Fig. 1 SPHERE/IRDIS polarimetric imaging observations of GSC 07396-00759. Shown from left to right are the total intensity Itot (log stretch) and the Qϕ and Uϕ images (both in linear stretch). The Qϕ and Uϕ images were convolved with a Gaussian point spread function (σ = 2 pixel) to increase the visibility in this plot. The coronagraphic mask is indicated by the green circular region in each panel. North is to the top and east to the left in each panel. It should be noted that the disk is not visible in the total intensity image (left panel). The structure appearing to extend from the NW to the SE is a point spread function artifact.

In the text

	Fig. 2 Signal-to-noise map estimated from the Qϕ and noise images (without convolution) of our IRDIS DPI observations. The plot is shown with a linear stretch between [− 3σ, 3σ]. The inner and outer model boundaries are indicated by the circle and the ellipse in the image. North is to the top and east to the left.
In the text

Fig. 3 Observation and model images for GSC 07396-00759. Shown from left to right are the measured Qϕ, best-fit model, and residual image obtained from our disk model using a linear scaling for each panel. The regions outside of the ellipse and within the circle (dashed green lines) are excluded from the χ2 calculation. The coronagraphic mask is indicated by the (shaded) circular region in each panel. North is to the top and east to left in each panel.

In the text

Fig. 4 Projected posterior distribution and density plots of the free parameters used in the disk modeling. The plot additionally shows the 50%, 16%, and 84% quartiles (vertical dashed lines), representing the distribution median and the 1 σ uncertainties (lower and upper bound), respectively. The 1D histograms represent the probability distributions of each parameter marginalized over the other free parameter. The contour maps represent the central 68.3%, 95.5%, and 99.73% of the 2D probability distributions of different combinations of parameters, marginalized over one another.

In the text

	Fig. 5 2D χ2 map constructed from the computed χ2 at each point (model) of the parameter grid. The free parameters of the model are the porosity and the grain size. The contours enclose the regions where the χ2 is smaller than the 5% (dashed line) and 15% (solid line) quartiles. To increase the contrast, the upper limit of the χ2 color map isset to the 50% quartiles.
In the text

Fig. 6 Polarized phase function (S12, HG) as a functionof scattering angle θ. Shown are the parameterized polarized phase functions from our best-fit disk model for g = 0.60 ± 0.09 (black line and gray area) and the range (red area, enclosing values between the 16% and 84% quartiles) of the best-fit model (s = 0.34μm, porosity = 0.01–0.99) of this phase function, considering a DHS. Additionally shown are the polarized phase functions resulting from the “low”-χ2 regions indicated in Fig. 5. The phase functions were sampled around grain sizes and porosities of about 1 μm and 0.80 (blue area) and grains larger than 100 μm and porosities below 0.6 (green area), respectively. The range of possible scattering angles (θ = 9.14° to θ = 170.86°) observable with the used mask is indicated by the blue shaded area.

In the text

	Fig. 7 2D map of the anisotropic scattering factor computed at each point of the grain-size, porosity parameter grid with the DHS model. The contours for (blue lines) and (red lines) indicate the range of scattering factors from our best-fit disk model with the HG approximation.
In the text

Fig. 8 Disk around GSC 07396-00759, as seen with SPHERE/IRDIS in polarized light (top-left panel) and total intensity (top-right panel; Sissa et al. 2018) in the H-band. The bottom panels show the best-fit disk model found in this work (bottom left) and found by Sissa et al. (2018) in the same scaling as the according observations. The dash-dotted line is parallel to the disk semimajor axis and crosses the central star.The coronagraphic mask is indicated by the (shaded) circular region in each panel.

In the text

Fig. 9 FWHM as a function of the center of the radial distribution r0 (left panel), and FWHM normalized by the radius of the disk as a function of r0 (right panel). Objects for which the FWHM was estimated readily from a radial Gaussian or power-law distribution are shown as circles, while objects where only the inner and outer radius were given are marked by diamonds. The spectral types of the sample stars are indicated by the color bar. The radius and width of the disk around GSC 07396-00759, estimated by Sissa et al. (2018) and in this work, are marked by the cross and the asterisk, respectively. To guide the eye, we also show the limits for FWHM = r0 and FWHM = 2r0 as a gray line in the left panel.

In the text

Fig. 10 Scattered light disk radii as a function of the stellar luminosity. Shown are the best-fit power-law function, as a solid red line, and the 1σ confidence interval drawn from the corresponding probability distributions. Data points that are also part of the sample analyzed byEsposito et al. (2020) are marked as diamonds. Similar to Esposito et al. (2020), we excluded HR 7012 (blue open circle) from the fit as an outlier. The best-fit result using the radius estimate by Sissa et al. (2018) (black star) is plotted as solid black line. Also shown for comparison are the radius-luminosity power laws for planetesimal belt central radii from thermal light imaging by Matrà et al. (2018) (dash-dotted line) and the GPIES-detected disks in scattered light by Esposito et al. (2020) (dashed line).

In the text

Fig. 11 β as a functionof grain size. βRP and β are shown as black and blue lines, respectively. The horizontal solid gray line (β = 0.5) is the upper limit for bound trajectories assuming zero eccentricity for the parent body (see Sezestre et al. 2017 for further details). The vertical dashed line indicates the corresponding minimum required grain size (s = 0.34μm) of dust particles to not be expelled from the disk, assuming a porosity of 70%.

In the text

	Fig. A.1 TESS light curve and periodogram of GSC 07396-00759. Left panel: Normalized TESS light curve of GSC 07396-00759, observed between June 19, 2019, and July 17, 2019. Right panel: Lomb-Scargle periodogram computed from these data, with an inset showing the light curve folded over the detected stellar rotation period Prot = 11.63 d.
In the text

	Fig. A.2 ASAS-SN light curve and periodogram of GSC 07396-00759. Left panel: Normalized ASAS-SN light curve of GSC 07396-00759, observed between March 10, 2016, and September 22, 2018. Right panel: Lomb-Scargle periodogram computed from these data, with an inset showing the light curve folded over the detected stellar rotation period Prot = 12.06 d.
In the text

	Fig. B.1 Distribution of possible scattering angle as a function of the model geometry obtained from our best-fit model of the disk. The inner and outer model boundaries are indicated by the circle and the ellipse in the image. Values within the circular mask or outside the ellipse are not considered in the modeling. North is to the top and east to the left.
In the text

Fig. C.1 Autocorrelation time as a function of MCMC iteration steps. The estimated average autocorrelation time is shown in blue (solid line). Also shown are the (red), Niter∕50 (blue), and Niter∕100 (black) lines to indicate different possible levels of acceptance for the model parameter estimates. For this work, the model is consideredconverged if and the change in consecutive estimated autocorrelation times τ is less than 1%.

In the text

Fig. D.1 Observation and model images for GSC 07396-00759. Shown from left to right are the measured Qϕ, best-fit model, and residual image obtained from our disk model using a linear scaling for each panel. The regions outside of the ellipse and within the circle (dashed green lines) are excluded from the χ2 calculation. The coronagraphic mask is indicated by the (shaded) circular region in each panel. North is to the top and east to left in each panel.

In the text

Fig. E.1 Projected posterior distribution and density plots of the free parameter (slope α, intercept , and fractional intrinsic scatter fΔR) of the power law fitted to the sample data points. The plot additionally shows the 50 %, 16 %, and 84 % quartiles (vertical dashed lines), representing the distribution median and the 1 σ uncertainties (lower and upper bound), respectively. The 1D histograms represent the probability distributions of each parameter marginalized over the other two. The contour maps represent the central 68.3%, 95.5%, and 99.73% of the 2D probability distributions of different pairs of parameters, marginalized over the third.

In the text