Quadrant polarization parameters for the scattered light of circumstellar disks. Analysis of debris disk models and observations of HR 4796A

This paper introduces the quadrant polarization parameters $Q_{000}$, $Q_{090}$, $Q_{180}$, $Q_{270}$ for Stokes $Q$ and $U_{045}$, $U_{135}$, $U_{225}$, $U_{315}$ for Stokes $U$ for the characterization of the azimuthal dependence of the scattering polarization of spatially resolved circumstellar disks illuminated by the central star. These parameters are based on the natural Stokes $Q$ and $U$ quadrant pattern produced by circumstellar scattering. They provide a simple test of the deviations of the disk geometry from axisymmetry and can be used to constrain the scattering phase function for optically thin disks without detailed model fitting of disk images. The parameters are easy to derive from observations and model calculations and are therefore well suited to systematic studies of the dust scattering in circumstellar disks. It is shown for models of optically thin and rotationally symmetric debris disks that the quadrant parameters normalized to the integrated azimuthal polarization or quadrant ratios like $Q_{000}/Q_{180}$ depend only on the disk inclination $i$ and the polarized scattering phase function of the dust, and they do not depend on the radial distribution of the scattering emissivity. Because $i$ is usually well known for resolved disk, we can derive the shape of the phase function for the sampled scattering angle range. This finding also applies to models with vertical extensions as observed for debris disks. Diagnostic diagrams are calculated for normalized quadrant parameters and quadrant ratios for the determination of the asymmetry parameter $g$ of the polarized Henyey-Greenstein scattering phase function. We apply these diagrams to measurements of HR 4796A, and find that a phase function with only one parameter does not reproduce the data well, but find a better solution with a three-parameter phase function.


Introduction
Circumstellar disks reflect the light from the central star, and the produced scattered intensity and polarization contain a lot of information about the disk geometry and the scattering dust particles. The scattered light of disks is usually only a contribution of a few percent or less to the direct light from the central star and therefore requires observations with sufficiently high resolution and contrast to resolve the disk from the star. This was achieved in recent years for many circumstellar disks with adaptive optics (AO) systems at large telescopes using polarimetry, a powerful high-contrast technique, to disentangle the scattered and therefore polarized light of the disk from the direct and typically unpolarized light of the central star (e.g., Apai et al. 2004;Oppenheimer et al. 2008;Quanz et al. 2011;Hashimoto et al. 2011;Muto et al. 2012).
With AO systems, the observational point spread function (PSF) depends to a large extent on the atmospheric turbulence and is highly variable (e.g., Cantalloube et al. 2019). For this reason, the disks are often only detected in polarized light and it is not possible to disentangle the disk intensity signal from the strong, variable PSF of the central star (see e.g., Esposito et al. 2020). Therefore, analyses of the scattered light from the disk are often based on the differential polarization alone and only in favorable cases can one combine this with measurements of the disk intensity. For the data analysis, the observed polarization must first be corrected for instrumental polarization effects, and this is relatively difficult for complex AO systems (e.g., Schmid 2021, and references therein). For this reason, the first generation of AO systems with polarimetric mode provided useful qualitative polarimetric results but hardly any quantitative results. This situation has changed with the new extreme AO systems GPI (Macintosh et al. 2014) and SPHERE (Beuzit et al. 2019), which, in addition to better image quality, also provide polarimetrically calibrated data for the circumstellar disk (Perrin et al. 2015;Schmid et al. 2018;de Boer et al. 2020;). Thus, quantitative polarization measurements are now possible for many circumstellar disks. However, the technique is not yet well established, and detailed studies have only been made for a few bright, extended disks; for example for HR 4796A (Perrin et al. 2015;Milli et al. 2019;Arriaga et al. 2020), HIP 79977 (Engler et al. 2017), HD 34700A (Monnier et al. 2019), HD 169142 (Tschudi & Schmid 2021), and HD 142527 (Hunziker et al. 2021). There also exist a few polarimetric studies based on HST polarimetry, such as those for the disks of AU Mic (Graham et al. 2007) or AB Aur (Perrin et al. 2009). The large majority of polarimetric disk observations in the literature focus their analysis on the high-resolution disk geometry (e.g., Benisty et al. 2015;Garufi et al. 2016;Avenhaus et al. 2018), and therefore the measurements are not calibrated and polarimetric cancelation effects introduced by the limited spatial resolution are not taken into account. Even for the detailed polarimetric studies mentioned above, the derived photo-polarimetric parameters are relatively heterogeneous. Convolution effects are only sometimes taken into account, and measurements are given as observed maps, azimuthal or radial profiles, or as dust parameters of a well-fitting disk model. This makes a comparison of results between different studies difficult and inaccurate, in particular because measuring uncertainties and model ambiguities are rarely discussed in detail.
Therefore, this paper introduces Stokes Q and U quadrant polarization parameters, which are simple to derive but still well defined, and facilitate systematic, quantitative studies of larger samples of circumstellar disks in order to make polarimetric measurements more valuable for our understanding of the physical properties of the scattering dust. Quantitative polarimetry of disks might be very useful to clarify whether dust properties are different or similar for systems with different morphology, age, level of illumination, or dust composition. The quadrant polarization parameters should also be useful for model simulations describing the appearance of a given disk for different inclinations so that intrinsic properties can be disentangled from the effects of a particular disk inclination. The same parameters can also be used to quantify the impact of the convolution of the intrinsic disk signal or of simulated images with instrumental PSF profiles to correct the observable polarization of small and large disks for the effects of smearing and polarimetric degradation (Schmid et al. 2006).
Using Stokes Q and U quadrant parameters for the description of the polarization of disks is a new approach and to the best of our knowledge this is the first publication using this type of analysis. These quadrant polarization parameters are introduced in Sect. 2 using the published data of Milli et al. (2019) for the bright debris disk HR 4796A as an example. The interpretation of the measured values is illustrated with corresponding model calculations for optically thin debris disks, which are described in Sect. 3. The models are similar to the classical simulation for the scattered intensity of debris disks of Artymowicz et al. (1989) and Kalas & Jewitt (1996), but they also consider the scattering polarization as in the models of Graham et al. (2007) and Engler et al. (2017). It is shown in Sect. 4 that the relative quadrant polarization values for optically thin, axisymmetric (and flat) debris disks are independent of the radial dust distribution and they depend only on the disk inclination i and the shape of the polarized scattering phase function f φ (θ). Therefore, in Sect. 5 we construct simple diagnostic diagrams for quadrant polarization parameters which constrain the f φ (θ) function for a given i, or directly yield the scattering asymmetry parameter g if we adopt a Henyey-Greenstein scattering function for the dust. The diagnostic diagrams are applied to the quadrant polarization measurements of HR 4796A and the obtained results are compared with the detailed, model-free phase-curve extraction of Milli et al. (2019). Finally, in Sect. 6 we discuss the potential and limitations of the quadrant polarization parameters for the analysis of disk observations and model simulations in a broader context.

Polarization parameters in sky coordinates
Polarimetric imaging of stellar systems with circumstellar disks provides typically sky images for the intensity I obs (α, δ) and the Stokes linear polarization parameters Q obs (α, δ) and U obs (α, δ). For dust scattering, the circular polarization is expected to be much smaller than the linear polarization and is usually not measured; it is therefore neglected in this work. Q and U are differential quantities for the linear polarization components Q = I 0 − I 90 and U = I 45 − I 135 , which can also be expressed as polarization flux P = p × I = (Q 2 + U 2 ) 1/2 and polarization position angle θ P = 0.5 atan2(U, Q) 1 . The polarized flux p × I is by definition a positive quantity, which, for noisy Q or U imaging data, suffers from a significant bias effect (Simmons & Stewart 1985). The azimuthal Stokes parameter Q φ can be used as an alternative for P for circumstellar disks. Q φ measures polarization in the azimuthal direction with respect to the central star (α 0 , δ 0 ) and U φ in a direction rotated by 45 • with respect to azimuthal. Circumstellar dust, which scatters light from the central star, mostly produces polarization with azimuthal orientations θ P ≈ φ αδ + 90 • , while U φ is almost zero. Therefore, Q φ can be considered as a good approximation for the polarized flux of the scattered radiation from circumstellar disks p × I = (Q 2 φ + U 2 φ ) 1/2 ≈ Q φ , and this approximation also avoids the noise bias problem (Schmid et al. 2006). For the model calculations of optically thin disks presented in this work, there is strictly p × I = Q φ , while small signals U φ < 0.05 Q φ can be produced by multiple scattering in optically thick disks (Canovas et al. 2015) or disks with aligned aspherical scattering particles. The U φ signal is often much larger, that is, U φ ∼ > 0.1 Q φ , in observations because of the PSF convolution problem for poorly resolved disks and polarimetric calibration errors. Both effects should be taken into account and corrected for the polarimetric measurements of disks. The azimuthal Stokes parameters are defined by with according to the description of Schmid et al. (2006) for the radial Stokes parameters Q r , U r , and using Q φ = −Q r and U φ = −U r .

Disk integrated polarization parameters
Quantitative measurements of the scattered radiation from circumstellar disks were obtained in the past with aperture polarimetry, which provided the disk-integrated Stokes parameters Q, U or the polarized flux P = (Q 2 + U 2 ) 1/2 usually expressed as fractional polarization relative to the system-integrated intensity Q/I tot , U/I tot , or P/I tot and the averaged polarization position angle θ p = 0.5 atan2(U, Q) (e.g., Bastien 1982;Yudin & Evans 1998). This polarization signal can be attributed to the scattered light from the disk if the star produces no polarization and if the interstellar polarization can be neglected, or if these contributions can be corrected. Usually, the disk intensity I = I tot − I star cannot be distinguished from the stellar intensity with aperture polarimetry, apart from a few exceptional cases, such as that of β Pic. Aperture polarimetry provides only the net scattering polarization, but misses potentially strong positive and negative polarization components +Q, −Q and +U, −U which cancel each other in unresolved observations. Therefore, Q and U, or P and θ p provide only one value and one direction for circumstellar disks, which agglomerate all possible types of deviations from axisymmetry of the scattering polarization of a circumstellar disk.
Disk-resolved polarimetric imaging avoids or strongly reduces the destructive cancellation effect and provides therefore much more information about the scattering polarization of disks. The most basic polarization parameter for the characterization of a resolved disk is the disk integrated azimuthal polarization Q φ which can be considered as equivalent to the polarized flux for resolved observations P 2 . The integrated polarized flux Q φ depends on the spatial resolution of the data; this aspect is not considered in this work because for well-resolved disks, like that of HR 4796A, the effect of limited spatial resolution is small and can be corrected with modeling of the instrumental smearing.
If the disk intensity I is measurable, one can also determine the disk-averaged fractional polarization p φ = Q φ /I. Unfortunately, it is still often very difficult to measure the disk intensity I(α, δ) with AO-observations because the signal cannot be separated from the intensity of the variable point spread function I star (α, δ) of the much brighter central star. In these cases, the fractional polarization can only be expressed relative to the total system intensity p φ,tot = Q φ /I tot . Integrated or averaged quantities are well defined but not well suited to characterizing the polarimetric features and the azimuthal dependence of the polarization for circumstellar disks. Therefore, we introduce new polarization parameters to quantify the individual positive and negative polarimetric components +Q, −Q and +U, −U of spatially resolved disk observations and models.

Quadrant polarization parameters in disk coordinates
For polarimetric imaging of circumstellar disks, the geometric orientation of the disk and the Stokes Q and U parameters are often described in sky coordinates. This is inconvenient for the characterization of the intrinsic scattering geometry of the disk and therefore we define new polarization parameters Q d and U d in the disk coordinate system (x, y), where the central star is at x 0 = 0, y 0 = 0 and the x and y are aligned with the major and  Milli et al. (2019). Left: Azimuthal polarization Q φ and Stokes Q and U in relative α, δ-sky coordinates; Right: Q φ , Q d , and U d in x, y-disk coordinates. minor axis of the projected disk, respectively. The positive xaxis is pointing left to ease the comparison with observations in relative sky coordinates α − α 0 and δ − δ 0 and to get the same convention for the Q d and U d orientations in x, y and sky images. Thus, the relative α − α 0 , δ − δ 0 sky coordinates of the observed images I, Q φ , Q, and U must be rotated according to as shown in Fig. 1 for the imaging polarimetry of HR 4796A using ω = 242 • (−118 • ). We use the convention that ω aligns the more distant semi-minor axis of the projected disk with the positive (upward) y-axis. Also, the Stokes parameters for the linear polarization must be rotated from the Q, U sky system to the Q d , U d disk system using the geometrically rotated (x, y)-frames We define for the Q d (x, y) and U d (x, y) polarization images the quadrant parameters Q 000 , Q 090 , Q 180 , and Q 270 for Stokes Q d and U 045 , U 135 , U 225 , and U 315 for Stokes U d , which are obtained by integrating the Stokes Q d or U d disk polarization signal in the corresponding quadrants as shown in Fig. 1 (e) and (f). This selection of polarization parameters is of course motivated by the natural Q and U quadrant patterns for circumstellar scattering where the signal in a given quadrant typically has the same sign everywhere and is almost zero at the borders of the defined integration region. This is strictly the case for all the model calculations for optically thin debris disks presented in this work and the same type of quadrant pattern is also predominant for the scattering polarization of proto-planetary disks. Multiple scattering and grain alignment effects can introduce deviations from a "clean" quadrant polarization pattern which may be measurable in high-quality observations (Canovas et al. 2015). However, the smearing and polarization cancellation effects introduced by the limited spatial resolution are typically much more important in affecting the quadrant pattern. For poor spatial resolution or very small disks, the quadrant pattern disappears (Schmid et al. 2006) or is strongly disturbed for asymmetric systems (Heikamp & Keller 2019).
Eight quadrant polarization values seems to be a useful number for the characterization of the azimuthal distribution of the polarization signal of disks. The parameters provide some redundancy to check and verify systematic effects, or to allow alternative disk characterizations if one parameter is not easily measurable or is affected by a special disk feature. Let us consider the redundancy in the context of the geometrical symmetry of disks and the dust scattering asymmetry.
For an inclined, but intrinsically axisymmetric disk, the Stokes Q d quadrants have the symmetry and the Stokes U d quadrants have the anti-symmetries Special cases generate additional equalities; for example the models with isotropic scattering (see Fig. 5) have a front-back symmetry and therefore there is also Q 000 = Q 180 and U 045 = −U 135 = U 225 = −U 315 . For an axisymmetric disk seen pole-on, all quadrant parameters have the same absolute value. If a disk deviates from an intrinsically symmetric geometry, for example a brighter +x side, then this would result in |Q 090 | > |Q 270 | for Stokes Q d and |U 045 | > |U 315 | or |U 135 | > |U 225 | for Stokes U d . Dust, which is predominantly forward scattering, produces more signal for front side polarization quadrants compared to the backside quadrants and this is equivalent to |Q 180 | > |Q 000 | or |U 135 | > |U 045 | and |U 225 | > |U 315 |. Properties that can be deduced from parameter ratios derived from the same data set are important because this reduces the impact of at least some systematic uncertainties in the measurements.
The quadrant parameters are also linked to the integrated polarization parameters Q φ , Q d , and U d . For the sum of all four Stokes Q d and Stokes U d quadrants, there is For pole-on systems, there is Q d = U d = 0, because of the symmetric cancellation of positive and negative quadrants and intrinsically axis-symmetric systems have U d = 0 for all disk inclinations because of the left-right antisymmetry. Axisymmetric but inclined systems have Q d 0 in general. For disks with larger inclination i, a smaller fraction of the disk is "located" in the quadrants Q 000 and Q 180 and a larger fraction is located in Q 090 and Q 270 near the major axis because of the disk projection. Edge-on disks are almost only located in the quadrants Q 090 and Q 270 and quadrant sums approach Q 000 + Q 180 → 0 and For the sums of the four absolute quadrant parameters for Stokes Q d , there is of course and the equivalent exists for sums of the Stokes U d quadrants.
For a pole-on view, the system is axisymmetric with respect to the line of sight and there is where each quadrant has the same absolute value of Q φ /2π = 0.159 Q φ . These sums will converge for edge-on disks i = 90 • to Σ |Q xxx | → Q φ for Stokes Q d and to Σ |U xxx | → 0 for Stokes U d .

Quadrant polarization parameters for HR 4796A
We use the high-quality differential polarimetric imaging (DPI) of the bright debris disk HR 4796A from Milli et al. (2019) shown in Fig. 1 as an example for the measurement of the quadrant polarization parameters. The data were taken with the SPHERE/ZIMPOL instrument (Beuzit et al. 2019;Schmid et al. 2018) in the very broad band (VBB) filter with the central wavelength λ c = 735 nm and full width of ∆λ = 290 nm. These data provide "only" the differential Stokes Q(x, y) and U(x, y) signals or the corresponding azimuthal quantities Q φ (x, y) and U φ (x, y), but no intensity signal because it is difficult with AO-observations to separate the disk intensity from the variable intensity PSF of the much brighter central star. The scattered light of the disk around HR 4796A was previously detected in polarization and intensity with AO systems in the near-IR (e.g., Milli et al. 2017;Chen et al. 2020;Arriaga et al. 2020), and in intensity in the visual with HST (e.g., Schneider et al. 2009).
Quadrant polarization parameters Q xxx and U xxx for HR 4796A derived from the data shown in Fig. 1(e) and (f) are given in Table 1 as relative values using the integrated polarized flux Q φ as reference. The quadrant values were obtained by integrating the counts in the annular apertures sections as illustrated in Fig. 2, which avoid the high noise regions from the PSF peak in the center. The uncertainties indicated in Table 1 account for the image noise, but do not account for systematic effects related to the selected aperture geometry or polarimetric calibration uncertainties. The noise errors are particularly large for the quadrants Q 000 and Q 180 because of the small separation of these disk sections from the bright star and additional negative noise spikes which are particularly strong for the Q 000 quadrant. Because of this noise, the formal integration gives Q 000 /Q φ = −0.04. From the decreasing trend of the signal in for example Q φ towards the back side of the disk and the measured disk signals of about 0.03 for U 045 /Q φ or U 315 /Q φ , an absolute signal of less than |Q 000 |/Q φ < 0.01 is expected for a smooth dust distribution in the ring and any reasonable assumptions for the dust scattering. Therefore, we do not consider the noisy Q 000 -measurement in the quadrant sums ΣQ xxx . The noise in the Q 180 -quadrant is also significantly enhanced when compared to other quadrants despite the relatively strong signal and the small integration area. A detailed analysis of the noise pattern might improve the measuring accuracy, but this is beyond the scope of this paper.
Because of the noise, the apertures for Q φ and the quadrants U 045 and U 135 were restricted to avoid the noisy and essentially signal-free region around φ xy ≈ 0 • . One should note that uncertainties for the quadrant measurements, for example for Q 180 , also affect the Q φ value and a dominant noise feature therefore has an enhanced impact on relative parameters such as Q xxx /Q φ , which include the disk integrated azimuthal polarization Q φ . These are typical problems for high-contrast observations of inclined disks and it can be very useful to select only high-signal-to-noise quadrant values for the characterization of the azimuthal dependence of the disk polarization.
Asymmetries between the left and right or the positive and negative sides of the x-axis can be deduced from the absolute quadrant parameters |Q 090 | and |Q 270 |, |U 045 | and |U 315 |, and |U 135 | and |U 225 |. Table 1 gives relative left-right brightness differences ∆ aaa bbb calculated according to for Stokes Q d parameters and equivalent for ∆ 045 315 and ∆ 135 225 for Stokes U d .
The asymmetry values ∆ 090 270 and ∆ 135 225 both yield more flux on the left side of the (x, y)-plane, which is the SW-side for the HR 4796A disk. This does not agree with previous determinations, including even the analysis of the same data by Milli et al. (2019) who measured more flux on the NE side. A more detailed investigation reveals that the peak surface brightness is indeed higher for the disk on the NE side and that the left-right asymmetry depends on the width of the annular apertures used for the flux extraction. If we integrate only a narrow annular region with a full width of ∆x = 0.1 at the location of the major axis then we also get more flux for the NE side or negative ∆ 090 270 -values of −3 ± 2 % but still less than the "negative" (SW-NE) asymmetry measured by Milli et al. (2019) for Q φ or Schneider et al. (2018) for the intensity. We therefore conclude that there are subtle asymmetries at the level of ∆ ≈ 10 % present for the disk HR 4796A, which are positive for a wide flux extraction and negative for a narrow flux extraction. Measurement uncertainties in the left-right differences are at the few percent level for bright quadrant pairs. The back-side to front-side brightness contrast can be expressed with quadrant ratios such as and their equivalents for other quadrant ratios. These ratios are small ∼ < 0.5 for HR 4796A which is indicative of dust with a strong forward scattering phase function. For isotropic scattering in an axisymmetric disk, the ratios would be Λ 000 180 = Λ 045 135 = 1 for all inclinations. Because the polarization flux in the backside quadrants is small, one can also assess the disk forward scattering with a comparison of the brighter Q 090 and Q 270 quadrants with the Stokes Q d front quadrant Q 180 or the Stokes U d front quadrants U 135 and U 225 as given in Table 1.
As demonstrated in Table 1, the polarization quadrants provide a useful set of parameters for the quantitative characterization of the geometric distribution of the polarization signal in circumstellar disks. Comparisons with model calculations are required to assess the diagnostic power of the derived values for the determination of the scattering properties of the dust or for the interpretation of the strength of disk asymmetries.

Disk model calculations
The simple disk models presented in this work follow the basic calculations for the scattered intensity from debris disks (Artymowicz et al. 1989;Kalas & Jewitt 1996) and the scattering polarization (e.g., Bastien & Menard 1988;Whitney & Hartmann 1992;Graham et al. 2007;Engler et al. 2017) but focus on the azimuthal dependence of the scattering polarization and the determination of the quadrant polarization values. The model disks are described by a dust-density distribution in cylindrical coordinates ρ(r, ϕ d , h) where r is the radius vector r = (x 2 d + y 2 d ) −1/2 and ϕ d the azimuthal angle in the disk plane 3 , where x d coincides with x for the major axis of the projected, inclined disk. We consider in this work very simple, optically thin τ 1, axisymmetric models ρ(r, h), where the dust scattering cross-section per unit mass σ sca is independent of the location. The total dust scattering emissivity (integrated over all directions) (r, h) of a volume element is given by the stellar flux F λ (r, h) = L λ /(4π(r 2 + h 2 )), the disk density ρ(r, h), and σ sca : The incident flux decreases as F λ ∝ 1/R 2 with R 2 = r 2 + h 2 , because in our optically thin scattering model we neglect the extinction of stellar light by the dust and the addition of diffuse light produced by the scatterings in the disk. The resulting images for the scattered light intensity I(x, y) and the azimuthal polarization Q φ (x, y) are obtained with a line of sight or z-axis integration of the scattering emissivity and the scattering phase functions f I (θ xyz , g) for the intensity and f φ (θ xyz , g, p max ) for the polarized intensity The transformations from the disk coordinate system (r, ϕ d , h) with x d = r sin ϕ d and y d = r cos ϕ d to the sky coordinate system (x, y, z) is given by This also defines the scattering angle θ xyz and the radial separation to the central star R xyz for each point (x,y,z).

Scattering phase functions
The scattering phase function f I (θ, g) for the intensity is described by the Henyey-Greenstein function or HG-function (Henyey & Greenstein 1941), where θ is the angle of deflection The asymmetry parameter g is defined between −1 and +1 and backward scattering dominates for negative g, forward scattering dominates for positive g, while the scattering is isotropic for g = 0.
The scattering phase function f φ (θ, g, p max ) for the polarized flux adopts the same angle dependence for the fractional polarization as Rayleigh scattering, but with the scale factor p max ≤ 1 for the maximum polarization at θ = 90 • . This description is often used (e.g., Graham et al. 2007;Buenzli & Schmid 2009;Engler et al. 2017) as a simple approximation for a Rayleigh scattering-like angle dependence but reduced polarization induced by dust particles (e.g., Kolokolova & Kimura 2010;Min et al. 2016;Tazaki et al. 2019).
The angle dependence of the fractional polarization of the scattered light is 3 the azimuthal angle ϕ d defined in the disk plane is different from the azimuthal angle φ xy , which is defined in the sky plane. The scattered intensity I sca can be split into the perpendicular I ⊥ and parallel I polarization components with respect to the scattering plane, so that I sca = I ⊥ + I , Q sca = I ⊥ − I and I ⊥ = (I sca + Q sca )/2, I = (I sca − Q sca )/2. Together with p sca , this yields or, expressed as scattering phase functions, f ⊥ = f I k ⊥ (θ, p max ) and f = f I k (θ, p max ), where k ⊥ and k are the expressions in the square brackets and p sca = k ⊥ − k . The scattering plane in a projected image of a circumstellar disk has always a radial orientation with respect to the central star. Therefore the induced scattering polarization Q sca , which is perpendicular to the scattering plane, translates into an azimuthal polarization Q φ for the projected disk map. The scattering phase functions are related by where f φ can be separated into a normalized part f n φ for p max = 1 and the scale factor p max . Figure 3 illustrates the scattering phase functions for g = 0.2 for the total intensity f I and the corresponding polarization components f ⊥ and f for p max = 0.5 and 0.2. The differential phase function f φ = f ⊥ − f has the same θ-dependence f n φ (θ, g) for both cases, and only the amplitude scales with p max . In this formalism, the total intensity phase function f I (g, θ) does not depend on the adopted polarization parameter p max . Expected values for approximating dust scattering with HG scattering functions are p max ≈ 0.05 − 0.8, but we often set this scale factor in this work to p max = 1 because this allows us to plot the intensity and polarization on the same scale. A value p max = 1 applies for Rayleigh scattering but the HG-function with g = 0 (isotropic scattering) differs from the Rayleigh scattering function for the intensity 4 . Hereafter, f φ (θ, g, p max ) and f n φ (θ, g) are also referred to as the HG pol -function and normalized HG pol -function, respectively.
The polarimetric scattering phase function for the azimuthal Stokes parameters f φ can be converted into phase functions f Q and f U for the Stokes Q d and U d parameters where φ xy = atan2(y, x) is the azimuthal angle in the sky coordinates aligned with the disk as illustrated in Fig. 1. These simple relations are valid in our optical thin (single scattering) models because U φ (x, y) = 0 and the Stokes Q d (x, y) and U d (x, y) model images can then be calculated as in Q φ in Eq. 18 but using the phase functions for f Q and f U .

Flat disk models and azimuthal phase functions
The calculations for the scattered flux I(x, y) and polarization Q φ (x, y) with Eqs. 17 and 18 are strongly simplified for a flat disk because the integrations for a given x-y coordinate along the z-coordinate can be replaced by single values for the separation r xy from the star, for the scattering emissivity (r xy ), and for the scattering angle θ xy . The volume density ρ must be replaced by a vertical surface density Σ(r) or a line-of-sight surface density Σ(r)/ cos i. Of course, the scattering in the disk plane must still be treated as in an optically thin disk, where κ Σ = a Σ + σ Σ is the disk extinction coefficient composed of the contributions from absorption a Σ and scattering σ Σ .

Projected flat disk image
The scattered intensity and polarization for a flat disk are given by and similar for the Stokes Q d (x, y) and U d (x, y) using the phase function from Eqs. 24 and 25. The scattering emissivity is proportional to the line-of-sight surface density Σ(r xy )/cos i, which considers the disk inclination, and where r 2 xy = x 2 + (y/ cos i) 2 and θ xy = acos(−z xy /r xy ) = acos(−y tan i/r xy ), because y = y d cos i and z xy = y d sin i.
In these equations for I(x, y), Q φ (x, y), Q d (x, y), and U d (x, y), the radial dependence of the scattering emissivity (r xy ) is separated from the azimuthal dependence described by the scattering phase functions f I (θ xy ), f φ (θ xy ), f Q (θ xy ), and f U (θ xy ). This is very favorable for the introduced quadrant parameters, which describe the polarization of the scattered light of disks with an azimuthal splitting of the signal Q d (x, y) → Q 000 , Q 090 , Q 180 , Q 270 and U d (x, y) → U 045 , U 135 , U 225 , U 315 by integrating the polarization in the corresponding quadrants outlined by the black lines in the Q d and U d panels of Figs. 4 and 5.
Disk images for I(x, y), Q d (x, y), U d (x, y), and Q φ (x, y) for flat disk models are shown in Fig. 4 for the asymmetry parameter g = 0.3, the scale factor p max = 1, and three different inclinations  i = 0 • , 45 • , and 75 • . For the radial surface scattering emissivity, a radial dependence (r) ∝ 1/r is adopted extending from an inner radius r 1 to the outer radius r 2 = 2 r 1 . Along the x-axis, the scattering angle is always θ xy = 90 • and the surface brightness increases for higher inclinations as in the inclined surface emissivity ∝ 1/ cos i.
The scattering asymmetry parameter g is relatively small and therefore the front-back brightness differences are not strong in Fig. 4. The forward scattering effect is much clearer in Fig. 5, where two disks are plotted with the same i and (r), but for isotropic scattering g = 0 and strong forward scattering g = 0.6.
For the disks shown in Figs. 4 and 5, the absolute signal drops in the radial direction for all azimuthal angles φ xy = atan2(x, y) by exactly a factor of two from the inner edge (r 1 ) xy to the outer edge (r 2 ) xy according to the adopted radial dependence of the scattering emissivity (r) ∝ 1/r. This is equivalent to the statement that the (relative) azimuthal dependence along the ellipse r xy describing a ring in the inclined disk is the same for all separations in a given disk image I(r xy , ϕ d ), Q φ (r xy , ϕ d ), Q d (r xy , ϕ d ), or U d (r xy , ϕ d ). This is also valid if the azimuthal an- gle ϕ d for the disk plane is replaced by the on-sky azimuthal angle φ. Therefore, it is possible to determine the azimuthal dependence of the scattered light in a very simple way using azimuthal phase functions defined in the disk plane, without considering the radial distribution of the scattering emissivity.

Scattering phase functions for the disk azimuth angle
The azimuthal dependence of the scattered light can be calculated easily for flat, rotationally symmetric, and optically thin disks as a function of the azimuthal angle ϕ d in the disk plane. For this, we have to express the scattering angle θ as a function of ϕ d and the disk inclination i according to where ϕ d = 0 for the far-side semi-minor axis of the projected disk.
The dependence of the scattered intensity and polarization flux with disk azimuthal angle ϕ d follows directly from the scattering phase functions f I and f φ and the change of the scattering angle θ as a function of ϕ d and i Figure 6 shows the azimuthal scattering functions 4π f I (ϕ d , i, g) for the intensity and 4π f n φ (ϕ d , i, g) for the polarized flux for g = 0.3 and for different inclinations. The factor 4π normalizes the isotropic scattering case 4π f I (ϕ d , i,g=0) = 1 and scales all other cases g 0 accordingly.
The enhanced forward scattering for f I (ϕ d ) around ϕ d = 180 • is clearly visible for inclined disks. The polarization function f n φ (ϕ d ) has for backward scattering around ϕ d = 0 • and forward scattering around ϕ d = 180 • strongly reduced values in highly inclined disks when compared to the intensity as can be seen for the green and red curves in Fig. 6. At ϕ d = 90 • and 270 • , the functions f I and f n φ have (for given g) the same value for all inclinations because the scattering angles are always θ ϕ,i = 90 • (cos ϕ d = 0 in Eq. 30).
The azimuthal dependence of the fractional polarization is given by and this function does not depend on the asymmetry parameter g. This dependence is equivalent to Rayleigh scattering and is shown in Fig. 6 for completeness.
For the phase function f Q and f U one needs to consider that the Stokes Q d and U d parameters are defined in the disk coordinate system projected onto the sky while f φ is given for the azimuthal angle ϕ d for the x d , y d -coordinates of the disk midplane. For the splitting of f φ (ϕ d , i, g, p max ) into f Q and f U , the azimuthal angle φ xy defined in the x-y sky plane must be used. The relation between ϕ d and φ xy is The phase functions for Q d and U d are then equivalent to the conversion given in Eqs. 24 and 25: These azimuthal function of the "on-sky" Stokes parameters f Q and f U is shown in Fig. 7. The functions are characterized by their double wave, which for i = 0 • are exact double-wave cosine f Q (ϕ d ) ∝ − cos 2ϕ d and double-wave sine f U (ϕ d ) ∝ − sin 2ϕ d functions. Deviations from the sine and cosine function become larger with increasing i, particularly for large asymmetry parameters g. The positive and negative sections of the f Q (ϕ d ) and f U (ϕ d ) functions correspond to the positive and negative polarimetric quadrants. The f Q and f U functions can also be expressed as normalized functions f n Q and f n U for p max = 1 and as fractional polarization p n Q and p n U equivalent to p n φ given above. For the adopted HG pol dust scattering phase function, these phase functions f ϕ provide a universal description of the azimuthal flux and polarized flux dependence as a function of inclination i for all flat, optically thin, rotationally symmetric disks.

Disk-averaged scattering functions
The total intensity I for our disk model can be conveniently calculated in r, ϕ d -coordinates because the integration can be separated between the r-dependent scattering emissivity and the ϕ d -dependent scattering phase function f I according to The first term represents the total scattering emissivity of the disk, and the second term is the disk averaged scattering phase function for the intensity f I (i, g) .  For the integrated polarization parameters, the same type of relation can be used The scale factor p max can be separated from the normalized versions of the disk-averaged scattering functions as and similar for f Q and f U . The disk-averaged fractional polarization follows from and similar for p Q or p U , where the latter is always zero for rotationally symmetric disks. Unlike for the azimuthal dependence of the fractional polarization p n (ϕ d , g), the disk-averaged parameters p n φ (i, g) depend on the g-parameter because in this average g shifts the flux weight between disk regions producing higher or lower levels of scattering polarization.

Normalized quadrant polarization parameters
The Stokes Q xxx and U xxx quadrant polarization parameters correspond to the individual positive and negative sections of the f Q (ϕ d ) and f U (ϕ d ) disk phase function shown in Fig. 7. The relation between quadrant parameters and phase function follow the same scheme as for the disk-integrated quantities Q d and U d described by Eq. 37 but the integration is limited to the azimuthal Table 2. Integration ranges for Stokes Q d and U d polarization quadrants for on-sky azimuthal angles φ x,y = atan2(y, x) and disk azimuthal angles The first term is the disk scattering emissivity integrated for the quadrant and the second term is the averaged f n Q scattering phase function for this quadrant. Because is independent of the azimuthal angle, the first term can be expressed as a fraction of the disk integrated emissivity · (ϕ 2 − ϕ 1 )/2π and the equation takes the form where we introduce the normalized quadrant polarization Q n xxx . The same scaling factors and p max are involved as for the equations for the integrated polarization parameters Q φ , Q d , and U d , and therefore Q n xxx and f Q (i, g) are related with the same factors and p max to observed quantities Q xxx and Q d . The same formalism applies to the Stokes U d quadrants.
According to Eq. 44, the normalized polarization quadrants Q n 000 , Q n 090 , Q n 180 , Q n 270 and U n 045 , U n 135 , U n 225 , U n 315 depend only on i and g as in the disk-averaged functions f n Q or f n U . The azimuthal integration range ϕ 1 to ϕ 2 is of course different for each quadrant, as summarized in Table 2. For the Q d -quadrants, the geometric projection effect introduces an inclination dependence for the integration boundaries ϕ 1 and ϕ 2 . For increased i, the fraction of the sampled disk surfaces s xxx = (ϕ 1 (i) − ϕ 2 (i))/2π in the left and right quadrants Q n 090 and Q n 270 is enhanced s 090 = s 270 = 0.5 − atan(cos i)/π and the disk surface fraction located in the front and back quadrants Q n 180 and Q n 000 is reduced s 000 = s 180 = atan(cos i)/π, respectively (we note that atan(cos 0 • ) = π/4 and atan(cos 90 • ) = 0). There is no i-dependence for the splitting of the Stokes U d quadrants, because the "left-side" back and front quadrants and the "right-side" back and front quadrants sample the same disk surface fraction s 045 = s 135 = s 225 = s 315 = 0.25 for all inclinations.

Calculation of disk polarization parameters
The previous section shows that the disk-integrated radiation parameters such as I, Q φ and the quadrant polarization parameters Q xxx and U xxx can be expressed as disk-averaged phase functions f (i, g) and normalized polarization parameters Q n xxx (i, g) Fig. 8. Disk-averaged scattering phase functions for the intensity 4π f I , the normalized (p max =1) azimuthal polarization 4π f n φ , and the normalized Stokes Q d flux 4π f n Q versus the disk inclination for HG-asymmetry parameters g = 0 (black), 0.3 (blue), 0.6 (red), and 0.9 (green). Lines give the results for flat disks and crosses for vertically extended disk rings (Sect. 4.3, Fig. 11). and U n xxx (i, g) describing the azimuthal dependence of the scattering, and scale factors for the total disk-scattering emissivity and the maximum scattering polarization p max . This simplicity allows a concise but comprehensive graphical presentation covering the full parameter space of the disk model parameters for the intensity and polarization, including the newly introduced quadrant polarization parameters. The Appendix gives a few IDL code lines for calculation of the numerical values. In addition, we explore the deviation of the results from vertically extended disk models and from the flat disk models.

Calculations of the disk-averaged intensity and polarization scattering functions
The disk-averaged scattering phase functions f (i, g) are equivalent to the basic integrated quantities I, Q φ , and Q d if normalized scale factors = 1 and p max = 1 are used (Eqs. 37 to 39). The functions f I (i, g) , f n φ (i, g) , and f n Q (i, g) are plotted in Fig. 8 as a function of inclination for the four asymmetry parameters g = 0, 0.3, 0.6, and 0.9. The plot includes the results from calculation of vertically extended, three-dimensional disk rings (crosses) for comparison, which are discussed in Sect. 4.3. The results for f (i, g) are multiplied by the factor 4π because this sets the special reference case for isotropically scattering dust g = 0 for all inclinations to 4π f I (i, 0) = 4π f n φ (i, 0) = 1 and simplifies the discussion. For pole-on disks i = 0 • , there is f n φ (0 • , g) = f I (0 • , g) for all g parameters because the scattering angle is θ = 90 • everywhere. For enhanced scattering asymmetry parameters g > 0 (but also for g < 0), the disk intensity is below the isotropic case 4π f I < 1 for lesser and moderately inclined disks i ∼ < 60 • , because the enhanced forward scattering (or backward scattering) produces enhanced scattered flux in directions near to the disk plane and reduced flux for polar viewing angles.
The forward scattering enhances the scattered intensity to 4π f I > 1 for i > ∼ 60 • and this effect becomes particularly strong for g → 1 and i → 90 • . This behavior is well known and produces a strong detection bias for high-inclination debris disks (e.g., Artymowicz et al. 1989;Kalas & Jewitt 1996;Esposito et al. 2020).
For the polarized flux f n φ (i, g) , an enhanced inclination does not produce an enhancement of the f φ signal because the strong increase in scattered flux from the forward scattering direction (or backward direction for g < 1) is predominantly unpolarized. For moderate asymmetry parameter |g| ∼ < 0.6, this causes an overall decrease of f φ with inclination (Fig. 8) while for extreme values |g| ∼ > 0.9 the huge flux increase compensates for the lower fractional polarization for forward and backward scattering.
The phase function for the Stokes Q d parameter f n Q (0 • , g) is zero for the pole-on view because of the symmetric cancellation of +Q d and −Q d signals. The function increases steadily with i (Fig. 8) and for i = 90 • or edge-on disks there is f n Q (90 • , g) = f n φ (90 • , g) , because all dust is aligned with the major axis and produces polarization in the +Q d direction.
Fractional polarization. The fractional polarizations p n φ (i, g) and p n Q (i, g) in Fig. 9 can be deduced from the ratio of the phase functions shown in Fig. 8. Observationally, a fractional polariza-tion determination requires a measurement of the integrated disk polarization and the integrated disk intensity.
As shown in Fig. 9, the fractional azimuthal polarization p n φ (i, g) only depends to a very small extent on g for small inclinations i ≤ 35 • with deviations < ±0.02 and the i-dependence is well described by A measurement of the fractional azimuthal polarization for lowinclination disks is therefore equivalent to a determination of the p max -parameter. The lower panel of Fig. 9 includes also the purely polarimetric ratio p Q / p φ = Q d /Q φ , which includes no scaling factor p max and systematic uncertainties from the polarimetric measurements might be particularly small. Therefore, the ratio Q d /Q φ can be used to determine the scattering asymmetry parameter g, if polarimetric cancellation effects for Q φ are taken into account for poorly resolved disks. For the extended disk HR 4796A, cancellation can be neglected, and we can use Q d /Q φ = 0.652 ± 0.030 (ΣQ xxx /Q φ from Tab. 1) and the inclination i = 75 • , which yields a value of about g = 0.7 from Fig. 9. This method is useful for systems with i ≈ 30 • − 80 • because the separations between the g-parameter curves are quite large. The curves in Fig. 9 for flat disk models are not applicable for edge-on disks i > ∼ 80 • with a vertical extension (see Sect. 4.3).

Calculations of the normalized quadrant polarization parameters
The normalized quadrant polarization parameters Q n 000 , Q n 090 , Q n 180 , and U n 045 U n 135 are plotted in Fig. 10 as a function of i for different g parameters. The "right-side" quadrants Q n 270 , U n 315 , U n 225 have the same absolute values as the corresponding "leftside" quadrants because of the disk symmetry. All values are multiplied by the factor 2π so that the reference case g = 0 and i = 0 • is set to 2π |Q n xxx (0 • , 0)| = 2π |U n xxx (0 • , 0)| = 1, similar to the normalization for the disk-averaged scattering function 4π f n φ (i, 0) . For pole-on disks, all the normalized quadrant polarization values have the same absolute value |Q n xxx (0 • , g)| = |U n xxx (0 • , g)| for a given g . This value is lower for larger g-parameter because less light is scattered perpendicularly to the disk plane in the polar direction, as in the f n φ function. For i > 0 • the quadrant values show different types of dependencies on i and g.
As described in Sect. 3.3, for the Stokes Q d quadrants the disk inclination introduces a geometric projection effect which increases the sampled disk area for the "left" and "right" quadrants for larger i , and therefore the Q n 090 and Q n 270 -values reach a maximum of (Q n φ /2) for i = 90 • . The areas for the front and back quadrants go to zero for i → 90 • and therefore so do the values Q n 180 and Q n 000 , but with a difference which depends strongly on g. The effect of the scattering asymmetry is already clearly visible for relatively small values, g ≈ 0.3, and inclinations, i ≈ 10 • , and becomes even stronger for larger g and i as can be seen from the enhanced |Q 180 | values relative to |Q 000 |. A similar front-back quadrant effect also occurs for the U d -components with enhanced absolute values for the front-side quadrant |U 135 | and reduced values for the backside quadrant |U 045 |.

Comparison with three-dimensional disk ring models
The polarization parameters derived in the previous sections are calculated for geometrically flat disks because this simplifies the calculations enormously. Of course, real disks have a vertical extension but observations of highly inclined debris disks typically show a small ratio h/r < ∼ 0.1 (see e.g., Thébault 2009). Therefore, the flat disk models could serve as an approximation for 3D disks, and we explore the differences. For this, we calculated models for a rotationally symmetric, optically thin disk ring with a central radius r ring and a Gaussian density distribution for the ring cross-section, with full width at half maximum (FWHM) of ∆ FWHM = 2.355 δ. Figure 11 compares images of such 3D disk rings with ∆ FWHM = 0.2 r ring , g = 0.6, and different i with flat disk models. The vertical extension of the 3D model is most obvious for the edge-on (i = 90 • ) case for which the flat disk model gives only a profile along the x-axis.
For i < 90 • , the differences between flat disks and vertically extended disks are already small for i = 75 • and hardly recognizable for lower inclination. In particular, the values for the diskaveraged scattering functions f and the normalized quadrant polarization parameters Q n xxx and U n xxx are equal or very similar as can be seen in Figs. 8 and 10 where the results from the 3D disks are plotted as small crosses together with those from the flat disk models. The agreement is typically better than ±0.01. An example of a systematic difference between 3D disks and flat disks is a slightly lower value (about 0.01) for the azimuthal polarization f φ in Fig. 8 for pole-on (i = 0 • ) 3D disks. For disks with a vertical extension, not all scatterings are occurring exactly in the disk midplane, but also slightly above and below where the scattering angle is smaller or larger than 90 • and therefore (1 − cos 2 θ)/(1 + cos 2 θ) in Eq. 20 is smaller than one. Another example is the reduced f I for edge-on (i = 90 • ) 3D disks, because less material lies exactly in front of the star where forward scattering would produce a strong maximum for large g-parameter. This is most visible for g = 0.6 (the effect is even stronger for g = 0.9 but this point is outside the plotted range). Similar but typically also very small effects are visible for the normalized quadrant polarization parameters in Fig. 10. The presence or absence of the vertical extension produces a strong difference for the Stokes U d which is also clearly apparent in Fig. 11 for i = 90 • .
These comparisons show that the quadrant polarization parameters derived from flat, optically thin, rotationally symmetric disk models are for most cases essentially indistinguishable from those derived from 3D disk models with a vertical extension typical for debris disks. Only for edge-on or close to edge-on disks, i > ∼ 80 • , can the vertical extension introduce significant differences, which needs to be taken into account.

Diagnostic diagrams for the scattering asymmetry g
The normalized quadrant polarization parameters and quadrant ratios for disk models using the HG pol -function depend only on the scattering asymmetry parameter g and the disk inclination i as shown in Fig. 10. From imaging polarimetry of debris disks one can often accurately measure the disk inclination, and therefore the quadrant polarization parameters are ideal for the determination of g. The method described in this work for the single Article number, page 11 of 21 A&A proofs: manuscript no. DiskQuadrants5 Fig. 10. Normalized (p max = 1, = 1) quadrant polarization parameters for 2π Q n 000 , 2π Q n 090 , 2π Q n 180 (left) and 2π U n 045 , 2π U n 135 (right) as function of disk inclination and for the HG-asymmetry parameter g = 0 (black), 0.3 (blue), 0.6 (red), and 0.9 (green). Lines give the results for flat disks and crosses for vertically extended disk rings (Sect. 4.3, Fig. 11) Fig. 11. Comparison of the 3D-ring model with the flat disk model. Left: I, Q d , U d and Q φ for a 3D disk ring with a full width at half maximum density distribution of ∆ FWHM /r ring = 0.2. Scattering parameters are g = 0.6 and p max = 1 for inclinations i = 0 • , 45 • , 75 • and 90 • . The same gray scale from +a (white) to −a (black) is used for all panels. The black lines in the Q d and U d panels indicate the polarization quadrants. Right: Same but for the flat disk model. parameter HG pol -function can be generalized to other, more sophisticated parameterizations for the polarized scattering phase function of the dust.
This diagnostic method is based on the strong assumption that the intrinsic disk geometry is rotationally symmetric, which is often a relatively good assumption for debris disks, but there are also several cases known with significant deviations from axisymmetry (e.g., Debes et al. 2009;Maness et al. 2009;Hughes et al. 2018). This can affect the determination of the g-parameter and one should always assess possible asymmetries in the disk symmetry using for example the left-right symmetry parameters ∆ 090 270 or ∆ 135 225 (Sect. 2.4).

Relative quadrant parameters
The azimuthal dependence of the scattering polarization can be described by the relative quadrant parameters Q xxx /Q φ and U xxx /Q φ . These parameters are directly linked to the normalized scattering phase functions f n Q , f n U , and f n φ for the polarization according to and equivalent to U xxx /Q φ . The relative quadrant polarization parameters |Q xxx |/Q φ and |U xxx |/Q φ are plotted in Figure 12 as a function of g for disk inclinations i = (0 • ), 15 • , 30 • , 45 • , 60 • , and 75 • . For pole-on disks, all relative quadrant values are equal to (2 π) −1 independent of the g-parameter, because of the normalization with Q φ .
With increasing g and i, the quadrant parameters show the expected steady increase for the front-side quadrants |Q 180 |/Q φ and |U 135 |/Q φ (red lines) and the steady decrease for the backside quadrants |Q 000 |/Q φ and |U 045 |/Q φ (blue lines). The separation between red and blue lines produces particularly large ratios for high inclination because the range of scattering angles extends from strong forward scattering to strong backward scattering.
For high inclination disks, the substantial contribution of forward and backward scattering strongly reduces the fractional scattering polarization in the front-side and back-side quadrants and the relative quadrant values are therefore < (2 π) −1 for small g. For the Q 000 and Q 180 quadrants, the reduction is further accentuated by the disk projection which reduces the sampled disk area for inclined disks. Strong forward scattering g → 1 compensates these two effects to a certain degree for the front-side quadrants (red curves), while it further diminishes the flux in the (blue) back side quadrants.
The two positive Q d quadrant parameters Q 090 /Q φ and Q 270 /Q φ depend mainly on the disk inclination. For higher inclination, a greater area of the disk is included in these two quadrants and therefore their relative contribution to the total polarized flux Q 090 /Q φ increases from (2 π) −1 for i = 0 • to 0.5 for i = 90 • . For edge-on systems, there is 2 · Q 090 = Q φ = Q d or all polarized flux of a disk is located only in the left and right quadrants Q 090 and Q 270 . Figure 12 includes the relative quadrant parameters measured for HR 4796A from Table 1. This disk has an inclination of about 76 • (e.g., Chen et al. 2020;Milli et al. 2019) and the i = 75 • -lines match this value well. For Stokes Q d the front side value |Q 180 |/Q φ is very sensitive for the determination of the gparameter because the corresponding curve in Fig. 12 is steep. This results in a value of g = 0.72 ± 0.03 with a relative uncertainty of only about ±4 %, despite the relatively large measuring error of about ±13 % for |Q 180 |/Q φ . The polarization signal is strong near the major axis and the corresponding quadrant values |Q 090 |/Q φ and |Q 270 |/Q φ can be measured with high precision of roughly ±3 %. However, because the |Q 090 |/Q φ curve is rather flat in Fig. 12, the resulting uncertainty on the derived g-values is also about ±3 %. The obtained g-values for |Q 090 |/Q φ = 0.62 and |Q 270 |/Q φ = 0.72 differ significantly because of the described deviation of the HR 4796A disk geometry from axisymmetry.
For the backside quadrant, only an upper limit of about |Q 000 |/Q φ < 0.04 could be measured. This limit is not useful for constraining the g-value because expected values for an inclination of 75 • are very low, namely |Q 000 |/Q φ ≈ 0.004 for g = 0.0 and ≈ 0.001 for g = 0.6, and the corresponding curve is below the plot range covered in Figure 12.
The relative quadrant values of HR 4796A for Stokes U d yield low asymmetry parameters g ≈ 0.4 for the back-side quadrants |U 045 |/Q φ and |U 315 |/Q φ , and larger values of g ≈ 0.5 and g ≈ 0.7 for the front side, with again a significant left-right asymmetry.

Quadrant ratios
The scattering asymmetry g can also be derived from quadrant ratios describing the brightness contrast between the front and back sides of the disk as shown in Fig. 13, which includes the measurements from HR 4796A.
High-quality determinations of g are achieved if the frontand the back-side quadrant polarizations can be accurately measured. For low-inclination systems, this should be possible for ratios like |Q 000 |/|Q 180 | or |U 045 |/|U 135 | where both the front side and back side are bright. For high-inclination systems, as in HR 4796A, the back side can be faint and the ratios |Q 000 |/|Q 180 | or |U 045 |/|U 135 | are small ( ∼ < 0.1) and therefore difficult to measure accurately. As an alternative, one can use ratios based on the bright quadrants like |Q 090 |/|Q 180 | or |Q 090 |/|U 135 | or equivalent ratios using the right-side quadrants Q 270 and U 225 . Many aspects of the diagnostic diagrams plotted in Fig. 13 are similar to the description of the relative quadrant parameters in the previous section.

Comparison of different g determinations
The measured quadrant polarization parameters for HR 4796A can be used to strongly constrain the asymmetry parameter g of the adopted HG pol scattering phase function f φ (θ), but only for the θ-range sampled by the used quadrant parameters (see also Hughes et al. 2018). The flux-weighted distribution of scattering angles θ sampled by a quadrant strongly depends on inclination i, but also on g as illustrated in Fig. 14. Indeed, all quadrants of a nearly pole-on disks probe f φ (θ) only near the scattering angle of 90 • , while some quadrants probe a large θ-range for strongly inclined disks.
We define the median angle θ med for the angle that represents the 50th percentile of a cumulative polarized intensity distribution covered by one quadrant. The back and front quadrants Q 000 and Q 180 sample a narrow range, the median angle θ med is close to the most extreme backward and forward scattering angle θ med ≈ 90 • ± i for a disk with inclination i, and θ med are essentially identical for different g. The ranges of scattering angles θ covered by the quadrants Q 090 , U 045 , and U 135 are very broad for larger i and the θ med depend significantly on g as shown in Fig. 14. For example, θ med (Q 000 ) is 90 • for isotropic scattering, and becomes smaller for i → 90 • and g → 1 as indicated by the colored θ med points (diamonds) and the dotted lines.
The quadrants U 045 and U 135 sample the back-and frontside parts of the disk and their θ med -angles lie between those of the Q d -quadrants. The front-side quadrant also shows a strong Fig. 12. Relative quadrant polarization values |Q 000 |/Q φ (blue), |Q 090 |/Q φ (black) |Q 180 |/Q φ (red) on the left and |U 045 |/Q φ (blue), |U 135 |/Q φ (red) on the right for flat disk models with different i and as a function of the scattering asymmetry parameter g. The measured values for HR 4796A are given on the right side in each panel and the corresponding g-parameter is given at the top or bottom using the i = 75 • curves for these derivations. tendency towards smaller θ med (U 135 ) values for larger g and i, as in the Q 090 -quadrant, while the g-dependence of θ med (U 045 ) is much smaller. Figure 15 shows the g-parameters obtained for HR 4796A as a function of the θ-angle probed by the used quadrant parameters. The angle θ corresponds for relative quadrant parameters to θ med for i = 75 • and g = 0.6 as given in Fig. 14 and the horizontal uncertainty bar spans two-thirds of the plotted θ-distribution (from the 16.6 to the 83.3 percentiles). For the quadrant ratios, the adopted θ-values are the mean of the θ med of the two quadrants and the horizontal bars illustrate their separation. In principle, one should consider for the θ med -values the systematic trend of g(θ med ) from higher g-values (≈ 0.7) for forward-scattering quadrants to lower values (≈ 0.4) for the backward scattering quadrants. We neglect this effect which would introduce θ med shifts of about −5 • for U 135 and U 225 , shifts of about +5 • for Q 090 and Q 270 , and smaller shifts for the other quadrants. Figure 15 shows for HR 4796A a systematic dependence of the derived g-parameters with scattering angle θ. The results from the relative quadrant parameters and the quadrant ratios are roughly consistent. The colors indicate measurements for the left (blue) and right (red) sides of the disk and g-values differ significantly between the two sides for θ med = 35 • and 70 • because of the left-right disk brightness asymmetry. On the fainter side, this effect reduces the derived g-value for U 225 /Q φ , while g is enhanced for Q 270 /Q φ because the intrinsic faintness of Q 270 mimics a disk with relatively little 90 • -scattering because of the normalization with Q φ . The intrinsic left-right brightness asymmetry of HR 4796A has less impact on the g determination based on quadrants ratios from the same side. This redundancy helps to disentangle the effects of the scattering asymmetry g from geometric or left-right disk brightness asymmetries.

A "mean" asymmetry parameter g for HR 4796A
The clear trend of the derived g-parameter with scattering angle θ for HR 4796A in Fig. 15 reveals that the used HG pol -function is an oversimplified description of the polarized scattering phase function for this object. For fainter or less well resolved disks, and for those with low inclination, it may not be possible to recognize such systematic deviations from a HG pol function, and for all these cases the derived g-value from the HG pol -function could serve as a good starting point for the analysis of quadrant polarization parameters.
Therefore, for HR 4796A we also derive a "mean" value for the HG asymmetry parameter g despite the discussed trend. To this end, for the seven measured relative quadrant values Table 1, we determine the bestfitting g-asymmetry parameter for the adopted disk inclination i = 75 • . This yields HG pol (i = 75 • , g = 0.65) with a weighted sum of squared deviations of χ 2 = 15.4, and the corresponding calculated and measured values are plotted in Fig. 16a. The differences between the models are more visible in Fig. 16b, where the deviations of data points and calculations from the best-fit model are shown. The large χ 2 -value indicates that the adopted HG pol -fit does not describe the data well because of the significant left-right asymmetry between |Q 090 | and |Q 270 |, or |U 135 | and |U 225 | at the level of about 4 σ (σ: standard deviations), which cannot be described with an axisymmetric disk model. Additionally, the best fit underestimates the quadrant values at small (θ med = 16 • ) and large (124 • ) scattering angles. Panel (b) also includes the best-fit results for slightly different disk inclinations HG pol (i = 73 • , g = 0.60) and HG pol (i = 77 • , g = 0.71), which differ very little and produce deviations between fit and data that are similar to the i = 75 • solution. The corresponding Henyey-Greenstein scattering phase functions HG pol (or f φ (θ)) are given in panel (c), while panel (d) shows the HG intensity function ( f I (θ)) for the best quadrant solution for i = 75 • and the solution from Milli et al. (2019).
Fortunately, we can compare the result from the quadrant parameter fitting with the analysis of the same HR 4796A data by Milli et al. (2019). They extracted from the polarimetric imaging data a detailed phase curve shown in Fig. 16c covering the scattering angle range θ = 13 • to 145 • . They also fitted the extracted phase curve with a single parameter HG pol -function and obtained an asymmetry parameter of g = 0.43 which is much smaller than our value of g = 0.65 (Fig. 16c). An important reason for this discrepancy is the sampling of the scattering angles of the data used for the fitting. In this work, the fitting is based on seven quadrant polarization values -for a disk with i = 75 • and g > 0.4-which are strongly biased towards small θ-values because of the forward "distorted" distribution of the polarized flux. In addition, the back-side quadrants U 045 and U 315 are weak and the corresponding measurements have a low signal-to-noise ratio of S /N < 5 and therefore a small weight. Thus, the fit to the quadrant values predominantly samples the range θ ≈ 16 • to 67 • of the scattering phase function. The analysis of Milli et al. (2019) samples a much broader range and particularly also more backward scattering angles. Therefore, the result of these latter authors of g = 0.43 closely matches the gvalues derived in this work by the quadrant ratios |U 045 |/|U 135 | and |U 315 |/|U 225 | (Fig. 13). On the other hand, the HG pol -fit of Milli et al. (2019) underestimates their extracted phase function in the forward-scattering range θ ≈ 16 • to 35 • . This comparison illustrates the bias effect that can be introduced by different kinds of phase curve sampling, if the adopted model curves f φ (θ) do not match well the real scattering phase function of the dust.
Compared to the case of the best single parameter HG pol function (Fig. 16), this fit does not underestimate the relative quadrant values at θ med = 16 • and 124 • and passes in the middle of the discrepant quadrant values at 35 • and 67 • for the left and right disk sides (see Fig. 17a and b). The corresponding polarized scattering phase function f φ (θ) in Fig. 17c has a much . For quadrants Q 000 and Q 180 only the case g = 0 is shown, because differences for other g-parameters are very small. The distributions are normalized for each quadrant individually. For the quadrants Q 000 , U 045 , and U 135 , the median angles θ med (diamonds) for large g-values increase with inclination as illustrated with the dotted lines. Fig. 15. Scattering asymmetry parameters for HR 4796A derived from relative quadrant parameters (open symbols) and quadrant ratios (filled symbols) as a function of the θ med or the mean of θ med , respectively. Colors indicate measurements from the left or SW (blue) and the right or NE (red) disk sides. wider peak extending from θ ≈ 20 • to ≈ 90 • , closely matching the directly extracted phase curve from Milli et al. (2019).
It is interesting to compare our results with the double HG pol function obtained by Milli et al. (2019) from the fit to the detailed phase-curve extraction which is included in Fig. 17  The good agreement between the double HG pol -fits of Milli et al. (2019) and the solution found for the quadrant polarization values shows that the selection of a more appropriate scattering phase function strongly reduces the large difference in the deduced g-determination described in Sect. 5.3.2 using only the single HG pol function. It should also be noted that the polarized phase function fit f φ (θ) of Milli et al. (2019) does not consider the left-right asymmetry of the disk in HR 4796A and therefore the phase curve uncertainty attained by these latter authors is larger than their measurement uncertainties. It seems likely that the azimuthal polarization signal extracted by Milli et al. (2019) would probably allow the determination of a better constrained empirical f φ (θ)-function for HR 4796A if the significant azimuthal dependence on the dust density is included in the fitting. Using a more detailed disk model for the fitting of the derived quadrant polarization parameters seems to be less useful because of the small number of measured values, a matter that is discussed further in Sect. 6.3. This example shows that selecting a good model fit function is important for the analysis of the quadrant polarization parameters and this should be investigated in more detail. The double HG pol is probably not an ideal choice, because in the intensity scattering function significant weight is given to the forward and backward scattering angles, which produce less polarization and therefore contribute less to the signal in the corresponding polarization quadrants Q 000 or Q 180 . This could explain the substantial differences for θ < 16 • or for θ > 124 • between the two derived best-fitting functions in Fig. 17d where the quadrant parameters provide no or only weak constraints on the shape of the scattering phase function.
In a future study, alternative polarized scattering phase functions f φ (θ) should be investigated for the fitting of polarimetric data, which give more weight to intermediate scattering angles θ ≈ 90 • ± i. Such a curve should also consider deviations of the fractional scattering polarization from a symmetric curve (Rayleigh-like) with respect to θ = 90 • as already derived from observations of HR 4796A by Perrin et al. (2015) and Arriaga et al. (2020). Considering this could be particularly important when constraining f φ (θ) for the dust in debris disks with smaller inclinations and a more limited observable range of θ-angles.

New polarization parameters for circumstellar disks
In recent years, the scattering light of many proto-planetary and debris disks has been spatially resolved with high-resolution polarimetric imaging using modern AO systems at large telescopes . Unfortunately, the presented results for the measurements of the polarized light from circumstellar disks are highly heterogeneous and are rarely flux calibrated, and it is therefore very difficult to compare the results from different studies for a systematic investigation of disks.
The main motivation of the present paper is the promotion of a photo-polarimetric parameter system which should help to ho- mogenize the polarimetric measurements for circumstellar disks and allow more straightforward comparisons between measurements of different disks and model results. The introduced quadrant polarization parameters Q 000 , Q 090 , Q 180 , Q 270 and U 045 , U 135 , U 225 , U 315 are defined for the Stokes Q d and U d parameters aligned with the apparent major and minor axis of the projected disk; they are based on the "natural" quadrant pattern produced by circumstellar scattering and measure within these quadrants the integrated Stokes Q d and Stokes U d flux, respectively.
These eight quadrants are very well suited for the description of the azimuthal dependence of the polarization signal of disks, except for edge-on or nearly edge-on systems. Furthermore, they can be used to quantify geometric deviations of the disk from axisymmetry from differences between left and right quadrants or characterize the disk inclination effects and the dust scattering asymmetry from ratios between back-side and front-side quadrants.
This disk characterization only requires differential polarization measurements, like relative quadrant parameters Q xxx /Q φ and U xxx /Q φ , or quadrant ratios like Q 000 /Q 180 . No absolute flux calibration with respect to the intensity of the star I star or the disk I are required and therefore one can also use polarimetric imaging of disks obtained in coronagraphic mode or with the central star saturated. In addition, the eight quadrant parameters are partially redundant and offer multiple options for the characterization of a disk, meaning that problems with a particular quadrant, for example because of the peculiarities of a disk or observational effects, can be mitigated.
The quadrant polarization measurements should be particularly well adapted for well-resolved, extended, low-surfacebrightness debris disks, which are relatively common (Espos-caused by the polarized scattering phase function f φ (θ) of the dust.

Limitations
The quadrant polarization parameters are designed for a simple description and analysis of the azimuthal dependence of the scattering polarization of circumstellar disks. The method is well defined but it has limitations, which must be taken into account in the interpretation of the results.
Importantly, one should be aware that real disks are often quite complex and a description using only eight parameters or less only yields rough information about the left-right disk asymmetry and the differences in the front-back brightness distribution. There are various effects that can cause significant departures from axisymmetry in the disk geometry: an intrinsic ellipticity introduced by noncircular orbits of dust particles, different types of hydrodynamic instabilities introducing spiral structures, azimuthal density features, lobsided disks, shadows cast by unresolved dust structures near the central star, dynamical interactions with proto-planets or other gravitating bodies in the system, and probably other effects.
Such asymmetries can be identified easily as left-right differences but they can also produce brightness effects between the disk front-and backside which are then blended with frontback brightness effects caused by the dust-scattering asymmetry in optically thin disks or the angle-dependent surface reflectivity in optically thick disks.
For a disk with complex morphology based on a small number of measured polarization quadrant parameters, it can be difficult to recognize whether the asymmetries are caused by the disk geometry, the scattering phase function, optical depth effects, or an observational problem. Therefore, it is certainly always useful to examine the disk polarization images for the presence of strong azimuthal structures which can be taken into account for disentangling the effects of the disk geometry from those of the dust scattering phase function for an interpretation of the data.
Using only eight quadrant polarization parameters for the characterization of the detailed structure of a well-observed disk can of course only provide limited information as demonstrated for HR 4796A. In such a case, a detailed extraction of the polarized flux or a two-dimensional model fitting to the data as in Milli et al. (2019) or Arriaga et al. (2020) will provide more accurate results and a less ambiguous interpretation. For example, a detailed model analysis for HR 4796A could consider two or three parameters for the intensity scattering phase function f I (θ), one or two parameters for the shape of the fractional scattering polarization p sca (θ) (and not only the fixed Rayleigh scattering like curve given in Eq. 20), a description of the ring geometry, and three or more parameters for the azimuthal dust density distribution.
Therefore, the quadrant polarization parameters are less suitable for a detailed investigation of well-observed disks, and are more suitable for the exploration and the approximate description of the global properties of the scattering dust in many disks. However, it is still useful to derive these parameters for wellobserved disk prototypes for a comparison with disks for which a detailed analysis is hardly possible, or for multi-wavelength studies of a given disk where a few well-defined parameters are sufficient to recognize and quantify wavelength dependencies for the polarized dust scattering phase function.

Conclusions
The quadrant polarization parameters introduced in this work seem to be very useful for a simple description of the azimuthal dependence of the polarization signal of circumstellar disks. These parameters can be determined from observations of many different types of circumstellar disks, for example debris disks around young or old stars, with or without strong illumination or dust blow-out signatures, or for proto-planetary disks with small or large central cavities and different kinds of hydrodynamical features.
The measured quadrant parameters can be compared with disk models that take the PSF smearing and cancelation effects into account and can explore the expected polarization signatures introduced by different descriptions for the scattering dust. Accumulating such data for a larger sample will allow a search for systematic trends in dust scattering properties for different disk types and for different wavelengths and inform us about the homogeneity or heterogeneity of dust-scattering properties in circumstellar disks. This can be achieved with relatively small uncertainties when compared to circumstellar shells or clouds, because the scattering angles θ, which have an important impact on the produced polarization signal, are typically very well known for resolved circumstellar disks. Investigations of the dust in circumstellar disks are also very attractive because many studies indicate that the dust evolves strongly in these systems and this could produce systematic trends for different disk types, which could be measurable with the new generation of AO polarimeters.