Multiband polarimetric imaging of HD 34700 with SCExAO/CHARIS

We present Subaru/SCExAO + CHARIS broadband (JHK) integral field spectroscopy of HD 34700 A in polarized light. CHARIS has the unique ability to obtain polarized integral field images at 22 wavelength channels in broadband, as the incoming light is first split into different polarization states before passing though the lenslet array. We recover the transition disk around HD 34700 A in multiband polarized light in our data. We combine our polarized intensity data with previous total intensity data to examine the scattering profiles, scattering phase functions and polarized fraction of the disk at multiple wavelengths. We also carry out 3D Monte Carlo radiative transfer simulations of the disk using MCFOST, and make qualitative comparisons between our models and data to constrain dust grain properties. We find that in addition to micron-sized dust grains, a population of sub-micron grains is needed to match the surface brightness in polarized light and polarized fraction. This could indicate the existence of a population of small grains in the disk, or it could be caused by Mie theory simulations using additional small grains to compensate for sub-micron structures of real dust aggregates. We find models that match the polarized fraction of the data but the models do not apply strong constraints on the dust grain type or compositions. We find no models that can match all observed properties of the disk. More detailed modeling using realistic dust aggregates with irregular surfaces and complex structures is required to further constrain the dust properties.


INTRODUCTION
Protoplanetary disks are the birthplaces of planetary systems.The initial conditions, stellar radiation, and dynamical interactions between planet formation and disk evolution shape the disk structures and influence the coagulation and growth of dust particles (Tazaki et al. 2019;Benisty et al. 2022).Observations of disk spectral energy distributions (SEDs) and spatially resolved imaging of disks at various stages of the disk evolution can inform us about the physical processes governing the disk evolution, and help develop self-consistent theories in disk evolution and planet formation (Dullemond et al. 2007;Andrews 2020).
Modern high contrast imaging cameras and techniques have become powerful tools to observe circumstellar disks in scattered light.For example, scattered-light observations of PDS 70 with the Very Large Telescope (VLT) in total and polarized intensity resolved the system's disk ring and inner disk and helped to distinguish these ★ E-mail: minghan@ucsb.edusignals from those of protoplanets (Keppler et al. 2018;Müller et al. 2018;Mesa et al. 2019;Haffert et al. 2019); similarly, analysis of Subaru Telescope and VLT observations of AB Aurigae resolved a morphologically complex disk with evidence for an embedded protoplanet(s) (Currie et al. 2022;Boccaletti et al. 2020, Dykes et al. 2024, submitted).The Gemini Planet Imager (GPI) recently observed a large sample of 44 bright Herbig Ae/Be stars and T-Tauri stars, yielding polarized scattered light detections in 80% of them and revealing a wide range of morphologies (Rich et al. 2022).High contrast imaging in polarized, scattered light has enabled new tests to validate radiative transfer models (Brauer et al. 2019;Arriaga et al. 2020;Bruzzone et al. 2020) and dust aggregate models of disks (Tazaki & Dominik 2022;Tazaki et al. 2023;Ginski et al. 2023).
HD 34700 A is a young spectroscopic binary system (Aa and Ab) with an estimated period of 23.5 days, originally thought to consist of two solar-mass stars and a debris disk, with an age ≳ 10 Myr (Torres 2004).However, the Gaia EDR3 distance places the system at a much further 350.5 pc (Gaia Collaboration et al. 2021) compared to the original ∼125 pc estimate, and the system was re-evaluated to  Torres (2004); 3. Mora et al. (2001); 4. Skrutskie et al. (2006); 5. Monnier et al. (2019); have an age of ∼5 Myr and a transition/protoplanetary disk (Monnier et al. 2019).Past observations of the disk include near-infrared (J and H) broadband total intensity and polarized intensity by the Gemini Planet Imager (GPI) Monnier et al. (2019), millimeter images by the Submillimeter Array (Benac et al. 2020), and integral field spectroscopy images in total intensity in JHK by the Coronagraphic High Angular Resolution Imaging Spectrograph (CHARIS) coupled with the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) (Uyama et al. 2020).These observations have revealed that the disk around HD 34700 A is rich with sub-structures: numerous spiral arms, a cavity at the north, and many darkening features.Most recently, the Spectro-Polarimetric High-contrast Exoplanet REsearch facility (SPHERE; Beuzit et al. 2019) at the VLT was able to fully resolve the inner disk around HD 34700 at a separation of ∼65 au to ∼117 au (Columba et al. 2024).
In this paper, we present new SCExAO + CHARIS integral field spectroscopic observations of the disk around HD 34700 A in polarized light.In Section 2, we describe our unique spectropolarimetric dataset and the data reduction processes.In Section 3, we present the reduced images of the disk, describe the disk morphology and features, and measure the surface brightness profiles and phase functions of the disk.In Section 4 we carry out 3D radiative transfer modeling of the disk at multiple wavelengths in JHK, and discuss our results and their implications.Finally, we summarize our work in Section 5.

CHARIS Observation Data
We present the polarimetric imaging results of HD 34700 A at J, H and K bands from CHARIS (Peters et al. 2012;Groff et al. 2014) coupled with SCExAO (Minowa et al. 2010;Jovanovic et al. 2015b).We also use an earlier total intensity dataset presented and analyzed by Uyama et al. (2020) for comparison and to calculate polarized fraction.The stellar properties of the HD 34700 A system are listed in Table 1.

Total Intensity Data
We use the total intensity data of HD 34700 A taken by Subaru/SCExAO + CHARIS on 2019 January 12 UT with a fixed pupil and a Lyot coronagraph mask to suppress starlight (Principal Investigator [PI]: Thayne Currie).The data were taken in broadband integral field spectroscopy (IFS) mode.In this mode, CHARIS covers a wavelength range of 1.16 − 2.37 m at a spectral resolution of  ∼ 19 with a pixel scale of 16.15 mas pixel −1 (Currie et al. 2023;Chen et al. 2023).The total intensity data was previously analyzed and published by Uyama et al. (2020).HR 2466 was also observed for a PSF reference for reference-star differential imaging (RDI; Lafrenière et al. 2009).The data were taken under very good seeing conditions (  ∼ 400 mas); the typical Full Width at Half Maximum (FWHM) was ∼ 30, 45, 55 mas in JHK bands, respectively.The total exposure time was 2168.6 seconds (30.98-second exposures × 70 cubes) for HD 34700 A and 2952.95seconds (1.475second single exposure × 14 coadds × 143 cubes) for HR 2466.We re-reduce the total intensity data for this work using an improved RDI PSF subtraction algorithm described in Section 2.2.3.

Polarized Intensity Data
We use two sets of polarized intensity data taken on two different dates.We will call them the primary polarized intensity data and the reference polarized intensity data.
Our primary data were taken on UT December 16 2019 (PI: Tim Brandt) with CHARIS broadband in polarimetry mode.In this mode, the incoming light is separated into two orthogonal polarization states by a Wollaston prism, and split into two beams by a beamsplitter, before entering the CHARIS lenslet array.The light then passes through a prism to be dispersed into a grid of microspectra on the detector.As a result, two 1 ′′ × 2 ′′ sections corresponding to the two orthogonal polarization states of the same field of view are dispersed onto the CHARIS detector, which originally covers a field of view of 2 ′′ ×2 ′′ .As the light is first separated into two polarization states and then passed through the integral field spectrograph, we can extract Stokes parameter images at the full spectral resolution of CHARIS broadband: 22 wavelengths spanning JHK for every spatial element in the 1 ′′ × 2 ′′ field of view.This is a unique advantage to CHARIS data, as most other polarimetry instruments can only observe at broad filters (e.g.GPI; Perrin et al. 2015;Monnier et al. 2019) SPHERE;de Boer et al. 2020) or at a specific wavelength (e.g.ZIMPOL; Columba et al. 2024).
The CHARIS polarimetry mode allows us to obtain full-frame polarized spectral data for HD 34700.The integration time for each exposure was 50.16 seconds for a total of 146 exposures (∼122 minutes).A half-waveplate was inserted and cycled through angles 0°, 45°, 22. • 5 and 67.• 5, with two exposures at each angle in every cycle.A Lyot coronagraph was used to block out the central binary.
In addition to our primary data, we use a reference CHARIS PDI dataset taken in the same polarimetry mode on UT 2019 Feb 26 (PI: Thayne Currie), but with the standard 4-satellite-spot pattern induced by SCExAO for astrometric and spectrophotometric calibration.The reference dataset has 284 exposures and an integration time of 7 seconds for each exposure.We use this reference dataset for flux calibration purposes to correct for systematics potentially caused by the "alternating grid" calibration pattern of our primary data (described in Section 2.2.2).We describe the systematics and the flux calibration in Section 2.2.4.The standard reduction and calibration of both PDI datasets follow the procedure described in Section 2.2.3.Then, our primary polarized intensity data is further corrected for systematics based on the reference polarized intensity data.The reference data has much less total integration time than the primary data, and covers a parallactic rotation of only ∼ 3°.Without enough field rotation, the reference data only imaged the disk partially due to the small field of view.Therefore, we use the reference data only for calibration purposes.The calibration procedure is described in detail in Section 2.2.4.

Polarized Flat Field Data
After cube extraction, the integral field spectropolarimetry data were adjusted with additional flat-field images, termed "PDI flat-fields," to correct effects introduced by the extra polarimetry optics (which are not included in the optical path for the standard CHARIS flat fields).PDI flat-field images for CHARIS are collected by observing SCExAO's halogen lamp calibration source (Jovanovic et al. 2016) with the Wollaston prism in place.For this purpose, the PDI flatfield images are extracted into cubes exactly as for the science data.The flat-field cubes are then median combined to create a single flat-field image cube.This cube is then normalized by dividing each wavelength slice by the median value across the slice.Subsequently, each science cube is flat-fielded by dividing it by this master PDI flat-field cube.Compared to the distributed flat fields in the CHARIS Data Reduction Pipeline (CHARIS DRP) (Brandt et al. 2017), these polarized flat fields show different structures for the two 1 ′′ ×2 ′′ field of views, and have variations in transmission of up to 10% which would bias the final Stokes parameters if not corrected.
As no PDI flat-fields were taken contemporaneously with the data, the nearest available calibration data from UT 2021 March 20 were adopted for this purpose.Comparison with more recent PDI flat-field images shows only small changes over time, such that the impact of time difference is expected to be negligible.

Cube Extraction and Selection Criteria
The raw data were first extracted into data cubes using the CHARIS DRP (Brandt et al. 2017).Each data cube consists of 2D images at 22 wavelength channels that span the range of J, H and K bands.The data cubes are spectrophotometrically calibrated to the bright primary stars using Kurucz model atmospheres (Castelli & Kurucz 2003), adopting spectral type G0V for HD 34700 A and A2V for HR 2466 Cutri et al. (2003).For the polarized intensity data, we experienced technical issues with optics alignment, a period of time with bad seeing, and varying adaptive optics (AO) system performance during the observing sequence.We visually inspect each exposure and manually remove ones with poor AO correction or seeing, and ones with instrumental artifacts.As a result, 70 exposures were discarded and 76 exposures were kept for the analysis.

Image Registration Using Alternating Satellite Spots
To locate the centroid of the central star, align the images in coronagraphic data, and calibrate the flux of the disk's surface brightness, a diffractive grid generated by sine waves on the deformable mirror (DM) of SCExAO in the pupil plane projects the primary star onto the focal plane with a fixed contrast at fixed separations and angles (Jovanovic et al. 2015a,b;Sahoo et al. 2020).Usually, there are four satellite spots at a fixed / separation that are roughly at the vertices of a square centered at the central star.However, the satellite spots contaminate a significant area of the disk in this case.To mitigate this contamination, we adopted an alternating pattern of satellite spots (Sahoo et al. 2020), where only two spots along one diagonal are projected in each exposure, and it alternates between the two diagonals over the sequence.This allows all parts of the disk to be free of satellite spots and therefore uncontaminated for half of the exposures, potentially providing better recovery of the disk signal, as well as better isolation of the satellite spots from the disk signal for a better flux calibration.The satellite spots are fitted and used to triangulate the centroid of the primary star in order to register the images.This is done using a centroiding algorithm described in Chen et al. (2023).The re-registered images are then used for PDI reduction.

RDI and PDI
Once registered, the PDI data cubes are further reduced following a similar procedure as that described in Lawson et al. (2021), but with the updates and caveats summarized hereafter.Initial spectrophotometric calibration is conducted following Lawson et al. (2021), but using only the cubes with visible satellite spots in every channel.Because the original field of view is split into two halves, the satellite spots in one of the two alternating diagonals move beyond the field of view after the seventh wavelength channel.Therefore, we only measure the fluxes of the satellite spots along the diagonal that stay within the field of view, and adopt this calibration for the exposures with the other configuration.The remaining cubes adopt the flux calibration scaling factors for the nearest calibrated cube in time.Results for this spectrophotometric calibration strategy were comparable to those achieved when first subtracting the alternating cubes to better isolate the flux of the satellite spots.Matching of half-waveplate cycles and the PDI procedure itself are carried out following the procedure outlined in Lawson et al. (2021).In place of the first order correction for instrumental polarization described therein, instrumental polarization is corrected using the full Mueller Matrix model summarized in Joost 't Hart et al. (2021).
To carry out PSF subtraction of the total intensity data from Uyama et al. (2020), we used the polarimetry-constrained reference star differential imaging (PCRDI) technique (Lawson et al. 2022).In PCRDI, available polarized intensity images are used to estimate the total intensity of the circumstellar signal (CSS) in order to suppress its impact on the creation of a stellar PSF model -and thus to eliminate the over-subtraction that typically plagues PSF-subtracted total intensity products.Our procedure for optimizing the CSS estimate in PCRDI diverges slightly from the one introduced in Lawson et al. (2022).Rather than optimizing a function of five parameters -four governing the scattering surface and an overall scaling factor -only the four parameters of the scattering surface are allowed to vary.For each trial solution, the resulting polarized intensity-based CSS estimate is treated like a disk model in RDI forward modeling and then compared with the residuals for a standard (unconstrained) RDI reduction for the real data.This allows the optimal scaling factor for each wavelength channel to be determined analytically by minimizing the sum of the squared residuals between the (scaled) forward modeled image and the data.Once the optimal parameters are identified, the corresponding CSS estimate is used to carry out PCRDI for the total intensity data.A comparison of PCRDI with the standard RDI is shown in Figure 1.The flux is scaled by the de-projected radius of the disk surface to make the difference more distinguishable visually.In PCRDI, the disk is visibly brighter and suffers less over-subtraction.The surrounding area is also much closer to zero, compared to over-subtracted negative areas in standard RDI.

Polarized Intensity Calibration
For coronagraphic data, the brightness of the astrophysical source is usually calibrated using the satellite spots and the known contrast of the spots to the central host star.Then, using the known magnitude of the host star we can convert the data units to physical units for the astrophysical source, which is the disk in this case.Since we are observing an extended disk signal rather than a point source, a fixed four-spot pattern would permanently overlay parts of the disk, potentially introducing biases the current data reduction pipeline is not equipped to correct for.Therefore, we opted to use an alternating twospot pattern described in Section 2.2.2, which leaves all parts of the disk unobscured for at least half the exposures.The contrast between the satellite spots and the host star for SCExAO is calibrated using an internal source over a narrow bandpass.This contrast measured at a specific DM amplitude and wavelength, and scales as  2  −2 , where  is the modulation amplitude of the DMs, and  is wavelength.The spot contrasts were measured to a precision of a few percent for both the regular four-spot pattern as well as the alternating two-spot pattern, and the grid amplitude and wavelength dependence are found to agree with the expected scaling (Currie et al. 2020;Sahoo et al. 2020).However, by adopting the calibration results in Sahoo et al. (2020), we find the calibrated surface brightness of the disk of our primary polarized intensity data to be systematically fainter by about a factor of 3 fainter than previously published polarized intensities in J and H bands Monnier et al. (2019), and also fainter by a similar factor than the reference polarized intensity data, described in Section 2.1.2.This systematic difference in the flux calibration is likely due to the unconventional alternating two-spot pattern.For the fixed fourspot pattern, the spot contrasts remain relatively stable over time between different tests.The stability of this mode has been quantified in Currie et al. (2020), where an ∼ 8% variation in the contrast was found between two calibration tests at different times.However, for the alternating two-spot pattern, no such investigation on contrast stability over different epochs has been carried out to our knowledge.We therefore assume the combined effect of hardware and software handling of the rarely used two-spot pattern is likely responsible for this large systematic shift in our flux calibration.To address this shift, we use the reference polarized intensity data described in Section 2.1.2to fit for a scaling factor.We use the mean-collapsed images of J and H bands of both datasets for calibration, as K band images have ∼ 5 − 10 times lower signal-to-noise ratios (SNR).For each collapsed image, we find the best-fit global scaling factor that scales our primary data to the reference data in polarized intensity.When fitting for this scaling factor, we consider a 10-pixel wide annulus centered at the best-fit ellipse to the disk's peak brightness, described in Section 3.1).To estimate the mean and error of the scaling factor, we bootstrap resample the wavelengths in J and H and create data cubes consisting of images at the bootstrapped wavelengths in J and H, then collapse the these data cubes and fit for the scaling factors for these bootstrapped samples.We find an overall calibration factor of 3.18 ± 0.06.This factor scales our polarized intensity data to the reference CHARIS data described in Section 2.1.2,observed in the same polarimetry mode, but with a diffractive grid configuration that has a robust flux calibration.The uncertainty represents only the scatter of the scaling factor among the bootstrapped samples, and is added to the other uncertainties of the flux calibration in quadrature.A small uncertainty here reassures that this systematic is not wavelength dependent.

IMAGING ANALYSIS
The double-differencing procedure described in Lawson et al. (2021) gives Q and U Stokes parameter images at each wavelength.These are collapsed in to J, H, and K broadband images shown in Figure 2. The first two columns show the typical quadrant-patterned Q and U Stokes images in J, H and K.The last column shows the radial Stokes images,   .Radial Stokes parameters,   and   , are defined as: (1) Radial Stokes parameters are first introduced in Schmid et al. (2006) as   and   , which are projections of Q and U onto the azimuthal and radial directions (±  ), and 45°angles from those directions (±  ).The radial Stokes parameters are useful for two main reasons: 1.In the case of single scattering events, the polarization would be entirely in the azimuthal direction and appear as +  , while   provides an estimate of the noise (Monnier et al. 2019); and 2. For noisy data, the polarized intensity will have a positive bias, while   does not.The noise in Q and U, originally centered around zero for a well-calibrated dataset, would now become strictly positive after taking the sum of the squares (Schmid et al. 2006).These properties make   a good alternative for the total polarized flux in Equation 3, and   a good proxy for the noise when analyzing disk signals.However, large optical depths, inclination, asymmetric scattering combined with PSF convolution can all influence the polarization and produce signals in   in many situations (Monnier et al. 2019;Canovas et al. 2015;Dong et al. 2016).For HD 34700, Monnier et al. (2019) found significant signals in   , which we also found in our CHARIS polarized dataset.The average SNR of the disk signal per pixel area (∼ 261 mas 2 ) in   is around 3-5 in the collapsed J, H, K images, with the brightest region in J band reaching a maximum SNR of around 10. Therefore, we do not use   as a proxy for the noise estimate, and instead use a more typical standard deviation approach described in Section 3.3.Furthermore,   is no longer equivalent to total polarized light,   , for this disk; however, it remains a good alternative to   , provided that we also project our models into   .Whether our signal in   results from real dust scattering physics or from geometric/PSF convolution effects requires careful analysis and modeling of the   component of the disk, which is beyond the scope of this paper.We restrict our analysis and modeling to the   component only.For simplicity, for the remainder of this paper, we will use "  " and "polarized intensity" interchangeably, and refer to "  /" as the polarized fraction, where  is the total intensity.We note that we adopt the radial Stokes conventions in Appendix A in Monnier et al. (2019), where polarization along the azimuthal direction corresponds to a positive signal in   .
Our polarized intensity images obtained simultaneously at 22 wavelengths from 1.16 m to 2.37 m, combined with the total intensity data from Uyama et al. (2020) taken in the same CHARIS broadband mode, allow us to examine the wavelength dependence of the polarization properties in more detail than was previously possible.We show the total and polarized intensity images collapsed into J, H, and K band in Figure 3, as well as the polarized fraction images obtained by taking the ratio between corresponding polarized and total intensity images.We also show false color composite renderings using these images, with J, H, and K images acting as the blue, green and red color channels in an RGB image, respectively.
To reduce the cost of computation in disk modeling and at the same time boost the SNR without losing too much resolution in wavelength, we binned our wavelength channels into slightly larger wavelength bins.First, we discard the two wavelength bins at 1.37 m and 1.87 m in which atmospheric transmission is low, and also the two reddest wavelengths at 2.29 m and 2.37 m which have SNR< 3 for most areas of the disk.We bin the remaining 18 wavelengths by groups of three to obtain 6 wavelength bins centered at 1.2, 1.3, 1.5, 1.7, 1.9, and 2.1 m.The data are mean-collapsed within each bin.Subsequent analysis and modeling are all carried out at these wavelengths.

Morphology of The Dominant Ring
The disk around HD 34700 A is dominated in scattered light by a bright ring ∼0.′′ 5 in radius and rich with structure.et al. 2020).This difference (∼ 2) is modest considering that the datasets are taken by different instruments and the ellipses are fit using different methods.In this section we revisit the basic morphology and orientation of the ring using our CHARIS data.
We adopt three approaches to fit for the position and morphology of the dust ring.For our first fitting method (abbreviated "lsq distance" in Table 2), we take a very similar approach to that used in Monnier et al. (2019) and Uyama et al. (2020).We divide the ring into small azimuthal sections and compile a list of the pixel coordinates of the local radial maxima for each section.We use the Python least-squares ellipse fitting software given by Hammel & Sullivan-Molina (2020), the same software used by Uyama et al. (2020).It performs a least-squares fit by minimizing the sum of the squares of the algebraic distances of the provided data points to the ellipse, and is shown to be robust against noise.Our second and third fitting methods maximize the sum of the fluxes of the pixels the ellipse passes through instead of minimizing the least-squares distance of a collection of pixels.In our second method ("max flux"), we take the average intensity of all the pixels a trial ellipse passes through, and maximize this average.The third method ("max flux interp") takes a set of evenly spaced (in terms of arclength) points on the ellipse, and interpolates the values at the exact coordinates of those points.The second and third methods maximize similar metrics, and thus produced similar results as expected (< 1).They agree with the first method to ∼ 2.
There are two major features of the disk that produce notable de-viations from an elliptical ring.One is the discontinuity on the north side of the disk; the other is an arc-like feature on the south side of the disk that is fainter in polarized intensity than its surroundings (an under-dense spiral arc).These features are visible in Figure 3, labeled F1 and F5 respectively.When selecting the local radial maxima around the ring for the first ellipse fitting method, the maxima at these features clearly deviate from the ellipse.Additionally, the flux changes at these features could also potentially bias the other two fitting methods that use flux maximization.Therefore, for all three methods, we mask out a 60°section at the north side and a 30°region at the south side.The north side mask covers the large cavity (F1) and a darkening feature (F3).The south side mask covers another darkening feature F5 features.Compared to the fits without the mask, the best-fit inclination and PA change between 0 − 2.The errors of the best-fit parameters for the first method are calculated via bootstrap resampling of the points selected for fitting.For the second and third methods, while we mask out two major disk features when fitting, the many arcs of the disk could have small effects on the radial profile that are not obvious by eye.Therefore, to estimate the uncertainties, we divide the remaining unmasked disk into 90 azimuthal sections with 3 degrees per section.We then carry out jackknife resampling by masking out one section at a time and fitting for the ellipse for each Jackknife sample, and estimate the uncertainties using these fits.
To account for systematic differences between these fitting methods, we take the error-weighted average of the three methods to derive our final reported disk parameters, summarized in Table 2.We also include the uncertainty of the CHARIS instrument's north pointing (Chen et al. 2023) in our errors on the inferred PA, and add all uncertainties in quadrature for our final error.The best-fit inclination and semi-major axis are in good agreement (< 1) with the best-fit values in Monnier et al. (2019) and Uyama et al. (2020).The best-fit PA, 74.• 3 ± 1. • 2, agrees with that in Uyama et al. (2020) within 1, but differs by ∼ 2 from the best-fit value in Monnier et al. (2019), 69.• 0 ± 2. • 3.

Disk Features
In Figure 3 we annotate the various features of the disk, including the cavity at the north side of the disk (F1), the darkening features (F2-F5), an arc that is part of the inner disk (F6), and an instrumental/reduction artifact (F7).The features labeled F1-F6 have been discussed in detail in past works Monnier et al. (2019); Uyama et al. (2020); Benac et al. (2020) and we do not focus on studying them in more detail in this paper.F6 shows a faint inner arc the J band images, which is more clearly seen in Monnier et al. (2019) and is confirmed to be part of a mis-aligned inner disk (Columba et al. 2024).The inner disk could cast shadows onto the outer disk, changing the intensity and scattering profiles Brauer et al. (2019).This poses major difficulties in the radiative transfer modeling and subsequent interpretation of the disk, which we discuss in more detail in Section 4. We observe that the wavelength dependence of the brightness of most of these features remain consistent with the rest of the disk and appear blue in the RGB image, suggesting broad consistency in the dust properties.The exceptions are the darkening features F4 and F5, at PA∼ 120°a nd ∼ 180°, respectively.In these regions, the wavelength dependence is noticeably shallower than the rest of the disk.F6 is an arc that is part of the inner disk, which is both too faint and too close to the inner working angle to be fully resolved by CHARIS.Finally, F7 is an arc-like feature that can be seen in the J band.It is more visible in the J band polarized fraction panel.In an un-collapsed cube with five wavelength channels in J band, this feature can be seen radially expanding through wavelengths as would be expected for a diffrac-  tion feature originating from the star.We consider F7 to most likely be an instrumental/data reduction artifact.

Surface Brightness Profiles
To measure surface brightness profiles of the disk, we first convolve the collapsed images in J, H, and K with a normalized Gaussian kernel with  = 1.0 pixels.The Gaussian kernel is smaller than the size of the instrumental PSF, and thus can reduce the noise and preserve the disk structures.This effectively measures the local surface brightness.We measure the surface brightness along the azimuthal direction.We take all the pixels intersected by the best-fit ellipse and measure the surface brightness and PAs of these pixels, where the PA of a pixel refers to the PA of the position of the pixel relative to the central star.The polarized fraction is defined as   /, where  is the total intensity.
To measure the errors of each data point, we first estimate a corresponding noise image for every collapsed image.We subtract from each image the median filtered image with a filter size of 3 pixels, which provides an estimate for the noise assuming there is negligible structure on the scale of the size of the filter.We divide this noise image into 11 concentric annuli around the central binary, from a radius of 10 pixels to a radius of 65, roughly corresponding to the radii of the inner working angle and the extent of the field of view.Thus, each annulus has a width of 5 pixels.We measure the standard deviation of the residuals for each annulus and use that as the noise estimate at the radius of the bin center.This gives us radially dependent noise estimates at 11 radii.Finally, we interpolate the noise at each pixel's radius to obtain the noise image.We then add the uncertainties of the flux calibration described in Section 2.2.4 to this noise image in quadrature.To properly propagate the noise for the convolved intensity scattering profiles, the squares of the noise images (i.e.variance images) are convolved with the squares of the Gaussian kernels.The noise image is about a factor of 2-3 times higher than the flux calibration uncertainty for the bright part of the disk, but they are on the same order of magnitude, ∼ 2%−5% relative to the bright part of the disk.This gives a typical SNR of ∼20-50 per pixel area in the bright parts of the disk.In the first two cases, the blue points are the pixels selected to calculate distance or flux.For the last case, the blue points are the equally spaced points on the ellipse for interpolation.The green lines are the best-fit ellipses to the blue points in each case.The best-fit values for semi-major axis, inclination and PA agree within 1 between the second and third methods, and they agree with the those of the first method to ∼ 2; these values are reported in Table 2.
bins described in Section 3. The shaded regions correspond to the disk features described in Section 3.2.The total intensity is measured from the dataset described in Section 2.1.1;this is the same dataset presented in Uyama et al. ( 2020), but processed with a different RDI algorithm described in Section 2.2.3 that tries to minimize over-subtraction of the disk signal.Our total intensity scattering profiles share the same overall shapes as those in Uyama et al. ( 2020), showing higher brightness in the forward scattering direction (north side) and a sharp dip at the north side cavity of the disk.However, our profiles show a bright portion at PA∼ −20°to −80°, as opposed to the more extended "dip" region in Uyama et al. (2020).This is likely due to the improved RDI reducing over-subtraction, as the brightness now goes back up to the model level beyond the cavity without a slightly extended dip shown in Uyama et al. (2020).The surface brightness monotonically decreases as a function of wavelength for both total intensity and polarized intensity.From J band to K band, the drop in brightness is about a factor of two, while the degrees of polarization remain similar.This trend is consistent with the results in Uyama et al. (2020), but differs from GPI data where the H band is brighter in both total and polarized intensity (Monnier et al. 2019).This is not too surprising considering the heavy dependence of peak brightness on the Strehl ratio.As a result, the polarized fraction is a better quantity for comparison.The polarized fraction as a function of PA (bottom panel of Figure 5) shows a similar shape as the polarized intensity, and is slightly higher at shorter wavelengths for the entire disk except in the regions with darkening features or the artifact.We do see agreement in wavelength dependence of the polarized fraction with the results in (Monnier et al. 2019).The peak polarized fraction usually appears near a scattering angle of 90° (Benisty et al. 2022) which should align with the semi-major axis under the assumption of a simple inclined circular disk.Our data suggests a different PA of the disk than the one inferred from ellipse fitting.Monnier et al. (2019) also found that the PA of the axis of peak polarized fraction is more indicative of the true PA of the disk based on radiative transfer modeling.We describe the PA inferred from the polarized fraction in more detail in Section 4.

Scattering Phase Functions
We measure the total intensity, polarized intensity and polarized fraction as a function of scattering angle using Diskmap (Stolker et al. 2016).These correspond to the apparent scattering phase functions of the disk.The scattering angle is defined as the angle between the incident light (from the star) and scattered light (towards the observer) in the scattering plane, implicitly assuming an optically thin disk.A scattering angle < 90°would represent forward scattering, while > 90°represents backward scattering.For a specified geometry of the scattering surface, Diskmap calculates the deprojected radius and scattering angle at every pixel coordinate of the disk.The observed brightness is then corrected by a factor of  2  , where   is the deprojected radius of a point on the scattering surface.Then, the scattering angles are binned into 60 bins over the available range of scattering angles.The east and the west half of the disk would roughly cover the same range of scattering angles.We do not treat them separately.This means that each scattering angle bin often consists of pixels from both the east and west side of the disk.The average intensity of the pixels in each bin and their standard deviation make up a data point and its error.For HD 34700, we measure these intensities over a narrow strip over the bright dominant ring.The strip is centered at the ellipse which is fit to the brightest part of the ring and is 6 pixels in width.We also mask the north side cavity from PA = 60°to 120°f or our measurements.We use the surface scale height of a flared disk: where  0 = 172.4au is the reference radius defined by the best-fit semi-major axis, ℎ 0 = 17 au is the disk scale height at the reference radius, and  = 1.125 is the disk flaring index.The ℎ 0 and  values are the best-fit disk model in Monnier et al. (2019).For the inclination, we use the inferred value from the best-fit ellipse from this work of 41°.For the PA, we use 92°which corresponding to our axis of maximum polarized intensity, instead of the 74.• 3 from our ellipse fitting.We discuss further in Section 4.2 how the models compare to the data and why the models favor the PA corresponding to the axis of maximum polarization instead of the semi-major axis of the  best-fit ellipse.Given these parameters, the scattering angle that we can probe for the disk ranges from ∼ 40°to ∼ 125°.
For single scattering events on a collection of particles, the intrinsic scattering phase functions represent the probability densities of scattering occurring at every angle, and would correspond to the scattering matrix elements (Bohren & Huffman 1998).However, for an optically thick disk, multiple scattering and limb brightening effects make the apparent scattering phase functions deviate from the intrinsic ones (Tazaki et al. 2023).Because of this, we are measuring the apparent scattering phase functions using Diskmap.We cannot measure the intrinsic ones without knowing the ratio of single scattering vs multiple scattering that we see at every point in the disk together with the angles that multiply scattered photons arrived from.But we can account for the effect of multiple scattering and limb brightening when simulating scattered disk images using Monte Carlo radiative transfer modeling.By measuring the apparent scattering phase functions of both the data and model images, we can compare them to gain insights on dust grain properties.We describe this modeling process further in Section 4. For simplicity, we will use the term "scattering phase functions" to refer to the apparent ones henceforth.
Figure 6 shows the final results of our phase function estimation.The top two panels show the measurements for irradiation-corrected total intensity  and azimuthal Stokes   , both normalized such that the peak scattering probability is 1.This is purely for aesthetic purposes.As we do not have access to all scattering angles, there is no way to properly normalize the integrated probability to 1.The bottom two panels both show the polarization fraction (  /) as a function of scattering angle, but assume different PAs of the disk relative to celestial north when defining the PA of each pixel coordinate.The third panel assumes a disk PA of 92°, corresponding to the PA of the axis of maximum polarized fraction.The last panel assumes a PA of 74.• 3, corresponding to the PA of the semi-major axis of the best-fit ellipse to the disk ring's peak brightness.The peak polarized fraction occurs at larger scattering angles when PA = 74.• 3. We discuss how our models compare to the data at these two PAs in Section 4.2.
The inner mis-aligned disk can cast shadows that produce features and alter the brightness of the outer disk (Brauer et al. 2019).Without knowledge of the inner disk properties and an estimate of its influence on the radiative transfer, we can only make limited observations of the total and polarized intensity phase functions.On the other hand, the polarized fraction primarily depends on the scattering properties of the local dust grains responsible for the scattering and is very weakly dependent on other common disk model parameters (e.g.dust mass and scale height; Dykes et al. 2024 submitted).Therefore, compared to the intensity phase functions, the polarized fraction is more informative of the properties of the dust grain population for our case.
The total intensity phase function is expected to favor forward scattering (small scattering angles), drop steeply in the first ∼ 30°, and have either a plateau or a slightly rising trend towards larger angles Benisty et al. (2022); Min et al. (2016).Due to the limited range of scattering angles we are probing, our total intensity measurement only captures the plateau part of the phase function.We see a high degree of polarization, with maximum polarized fractions of around 0.6-0.8across J, H, and K bands.We also see that the polarized fractions are higher at longer wavelengths, peaking at ∼0.8 in K band, ∼ 0.7 in H band, and ∼0.6 in J band.We discuss in more detail how our models match these observations of the data in Section 4.

DISK MODELING WITH MCFOST
In this section, we carry out modeling of the outer disk, resolved in CHARIS images in scattered light, at the wavelengths corresponding to the data wavelength bin centers.We use the 3D Monte Carlo radiative transfer code with ray-tracing, MCFOST (Pinte et al. 2006(Pinte et al. , 2009(Pinte et al. , 2022)), to produce modeled images in scattered light.
MCFOST computes the optical properties of dust using Mie theory, typically treating the grains as homogeneous spheres (Pinte et al. 2006;Bohren & Huffman 1998).It also has the capability to add in porosity, volume mixing of different components, or to change the grain type to a distribution of hollow spheres (DHS).An effective optical index of the grains is computed using the Bruggeman rule before any Mie theory computation.
Compared to more exact methods such as the superposition Tmatrix (Mackowski & Mishchenko 1996) or the discrete dipole approximation (DDA; Purcell & Pennypacker 1973;Draine & Flatau 1994), Mie theory has limitations in how accurately it can produce scattered light images using its approximation of optical properties of dust grains Min et al. (2016).But in exchange, the reduced com- The scattering angle is defined as the angle between the incident light from the star and the scattered light to the observer, with forward scattering having angles < 90°.We show two polarized fraction phase functions measured assuming different PAs in the third and fourth panels.PA=92°c orresponds to the PA of the axis of maximum polarization, and PA=74.• 3 corresponds to the PA of the semi-major axis of the best-fit ellipse described in Section 3.1.The shaded regions indicate the approximate scattering angles covered by the labeled features in Figure 3 except F6, which is not part of the outer disk.
putational cost allows us to explore a parameter grid with more disk parameters.
The inner disk is only partially detected in CHARIS images and we cannot infer its properties from our data.While the inner disk contributes a negligible amount (∼ 1%) in thermal emission compared to the starlight at our observing wavelengths (1 − 2.1), we cannot account for the effects of shadowing and scattered star light by the inner disk.Because of this, we focus on examining the polarized fraction of our data and models.The degree of polarization is insensitive to shadowing effects; in even more morphologically complex disks than HD 34700 A's (e.g.AB Aurigae), the degree of polarization is also only very weakly dependent on other common disk model parameters (e.g.dust mass and scale height; Dykes et al. 2024 submitted).Therefore, the degree of polarization primarily depends on the dust grains responsible for the scattered light, and Mie theory simulations still offer us insights on dust properties (Min et al. 2005;Mulders et al. 2013;Arriaga et al. 2020).In addition, an examination of the degree of polarization at finer wavelength resolutions was not possible for J, H, and K broadband datasets in the past, and is now made possible with the addition of our CHARIS PDI dataset combined with a previous CHARIS total intensity dataset, which we binned into 6 wavelengths from 1.2 to 2.1.In this section, we carry out modeling under Mie theory of the outer disk resolved in our data and discuss the results and their implications.

Model Parameters and Modeling Strategy
We adopt the stellar properties of the central binary from Monnier et al. (2019), including stellar radii, effective temperatures, and masses.We update the distance of the system according to the Gaia EDR3 measurement: 350.5 +2.5 −2.4 pc (Gaia Collaboration et al. 2021).We model the disk as a circular, single-component, axisymmetric, flared disk around the central binary, observed at a non-zero inclination.The disk's surface density is defined as: Where  in and  out are the radii of the inner and outer edges of the disk.The scattering scale height is defined by Equation 4 in Section 3.1.And finally, the number density of the dust population is defined as: where  is the radius of the dust grain, and  is the grain size powerlaw index (Dohnanyi 1969;Mathis et al. 1977).Due to the aforementioned limitations of Mie theory and simplifications in the disk model, a conventional  2 goodness-of-fit metric of the model image to the data image would not be a reliable way to evaluate the model and does not provide further information beyond that gained from qualitative comparisons.Therefore, we simply select by eye the models that best match the data qualitatively, excluding the cavity region on the north side and the darkening features.This essentially requires the scattering profiles of the models to have shapes resembling those of the data, while leaving a generous room for the unaccounted effects.Then, we consider models with polarized fractions within ∼ 20% across all scattering angles to be an acceptable fit.Due to the qualitative nature of the comparison, instead of finding a best-fit model, we present in section 4.2 one representative model that is the best qualitative match to the data.Alongside the qualitative match, we also show and discuss models that are qualitatively different, what parameters led to these qualitative changes and what the implications are on dust grain properties.
We carry out the modeling by constructing a coarse parameter grid for the disk parameters, as a high-dimensional fine grid is prohibitively expensive in computation time.We fix the parameters that are constrained by other independent observations: stellar properties and dust mass are fixed at the values in Monnier et al. (2019).Distance is fixed at the Gaia EDR3 mean value, 350.5 pc (Gaia Collaboration et al. 2021).We also fix parameters that can be well constrained by imaging analysis alone, such as the inner radius,  in , and inclination, .Finally, to reduce the number of dimensions and computational cost, some parameters are fixed at their canonical or approximate values after testing that they negligibly influence the simulated images within a physically reasonable range.These fixed parameters include the surface density exponent, , and outer radius,  out .Our final parameter grid is shown in Table 3.We use two dust populations for our models following previous efforts by Monnier et al. (2019); Uyama et al. (2020), one for small grains and one for large grains.We label the minimum and maximum grain sizes  min and  max , and the transition size between the two populations   .Both populations follow the same number density distribution given by Equation 6 with the same  (i.e. the two grain populations share the same number density slope, with a discontinuous transition at   ).For each model, the two populations have a fixed mass fraction.This means that for each model, a fixed amount of dust mass will be distributed between  min and   , while the remaining dust mass will be distributed between   and  max .We simulate both Mie Bohren & Huffman (1998) and DHS Min et al. (2003, 2005) grains, and include two different compositions in our models: astronomical silicates Draine & Lee (1984) and an amorphous carbon mixture following the prescription given by Min et al. (2011) with a carbon partition parameter  = 1.

Model Results
We find that we cannot identify a model that fits well (within a few percent in regions without disk features) to all three quantities: 1. total intensity; 2. polarized intensity; 3. polarized fraction.However, this is not surprising given the limitations we have discussed in Section 4. We summarize the qualitative influence of the model parameters on the scattering profiles in this section.
For every set of model parameters in the parameter grid, MCFOST outputs model images in total and polarized intensity.We convolve each model with the instrumental PSF, and then measure the surface brightness profiles and scattering phase functions following the same methods in Section 3.3 and Section 3.4, respectively.Because the flaring index of the disk surface and the PA for each model can vary, the geometry of the scattering surface can change for each model.Since the scattering phase functions presented in Section 3.4 assume a specific scattering surface scale height and PA, they cannot be directly compared to those measured from the model images, which have different scattering surfaces.Therefore, to properly compare the scattering phase functions between each model and our data, we re-measure the scattering phase functions for the data using the scattering geometry of each model.We find that the models can closely reproduce the polarized fraction phase function of the data when assuming PA = 92°, for which the polarized fraction peaks between scattering angles of 95°− 110°.If we use PA = 74.• 3, the polarized fraction of the data skews towards large scattering angles and peak around 100°−115°as shown in Figure 6, which is a slightly worse fit to the models.Because the polarized fractions from radiative transfer models are grounded in scattering theory and are determined by the scattering phase functions of the dust population, it is more robust than the PA inferred from a simple ellipse fit to the disk ring's peak brightness, which cannot capture the complexities of the various disk features.This agrees with the conclusion in Monnier et al. (2019) that the axis of peak polarization is more indicative of the disk's true PA, as opposed to the PA of the semi-major axis of the best-fit ellipse.The solid lines are the measurements for the PSF-convolved model images, and the data points with error bars are measurements for our datasets using the scattering surface geometry of the corresponding models.Each row shows a different set of parameter models from the model grid, with key parameters annotated in the right column plots.We find no model that can simultaneously fit the total intensity and polarized intensity, which is not surprising considering the significantly simplified model compared to the feature-rich disk, in addition to the effect of the inner disk that we ignored.However, as the polarized fraction reflects the local properties of the dust particles responsible for the scattering, we are able to find models that fit the polarized fraction and they suggest significant masses in sub-micron as well as micron sized grains.
In Figure 7, we show a set of models with distinct qualitative features in their scattering profiles.We find that the number density exponent had a small effect on the models within the range [3.2, 3.8], and its effect is degenerate with changing the mass fraction and   to produce a dust population with similar effective optical indices.Therefore, in Figure 7, we choose to present models with a number density parameter of 3.5 (Dohnanyi 1969;Mathis et al. 1977) among the ones with similar qualitative features.The top-row panel shows a model that is the best qualitative match to the data, which still cannot fit the total and polarized intensity very well.The subsequent rows show models that have distinct qualitative features different from the data, which inform us about one or multiple physical properties of the dust population of the disk.

Qualitative Fit
The top row of Figure 7 shows a model that qualitatively best matches the data by eye, which still cannot match the total and polarized intensity of the data well.This is not surprising given the simplified model compared to the feature-rich disk, in addition to the effect of the inner disk that we ignored.However, as the polarized fraction reflects the local properties of the dust particles responsible for the scattering, it still provides meaningful insights on the grain properties.Indeed, the polarized fraction of the model and data match in terms of the overall degree of polarization, the wavelength dependence, as well as the scattering angle dependence.We find that significant mass fractions in both sub-micron sized grains and micron-sized grains are required to reproduce the polarized fraction as a function of scattering angle.As a result, a single population with a constant number density power-law cannot reproduce the shape and the level of surface brightness seen in the data.The small grains are responsible for more isotropic Rayleigh scattering and are needed to produce enough brightness at larger scattering angles.

Models with Insufficient Sub-micron Dust Grains
The second row of Figure 7 shows a model with much less mass in the small grain population.While we can still match the polarized fraction with a higher porosity, we see that the lack of small grains produces much less polarized light at larger PAs corresponding to 90°+ scattering angles.The large grains are responsible for skewing the surface brightness towards extreme forward scattering and insufficient brightness for the south side.

Models with Insufficient Micron Sized Grains
The third row of Figure 7 shows a model where the dividing line between small and large grains is at 2.Similar to the model shown in the second row, the excess in forward scattering remains a problem.In addition, with a number density power-law of −3.5, this leads to a much smaller fraction of grains with sizes comparable to the observing wavelengths ∼ 0.5 − 2.0.Grains much larger than our observing wavelengths are less efficient at scattering incident light, resulting in much lower surface brightness, as well as a lower polarized fraction.The first three rows show that it is necessary to have significant mass fractions of grain sizes that are sub-microns and microns to simultaneously produce the high degree of polarization, bright peaks at east and west sides of the disk, and sufficient overall brightness in total intensity and polarized intensity.It is difficult to explain the physical processes that can generate this double population with a discontinuous transition at the grain size  min .It has been found in laboratory experiments that simulations using Mie/DHS particles can suggest a sub-micron grain population that is not found in the measured dust distribution in the lab (Min et al. 2005).It is possible this reflects the rough surfaces with sub-micron-scale structures of real, irregular dust grains (Min et al. 2005).

Models with Different Grain Types and Compositions
The fourth and fifth rows of Figure 7 show models using DHS grains and Mie grains with a carbon-rich composition, respectively.Within the model parameter space we explored, we find that these two grain types/compositions can produce matching polarized fractions, but struggle with forward scattering excess more than Mie grain models and produce less scattered light in general.The silicate DHS grain models otherwise share similar results as the silicate Mie grain models.The Mie grain model with amorphous carbon generates a higher degree of polarization, such that fewer small grains and less porosity are needed to match the polarized fraction.Given the limitations of our models and the small parameter space we are able to explore, the relatively worse fits to the total and polarized intensities do not offer strong arguments against DHS grains or grains with amorphous carbon.Future work with more quantitative modeling and access to potentially fully resolved images of both the inner and outer disk may be able to provide more insights regarding grain types and composition.

SUMMARY AND DISCUSSION
We have presented in this paper spectro-polarimetry of the transition disk around HD 34700 A obtained using Subaru/SCExAO +CHARIS broadband (J, H, and K).We fit a simple ellipse to the disk ring using two different methods; assuming a circular disk, we obtain an inclination of 41°± 1. • 1 and a PA of 74.• 3 ± 1. • 2. Combining our polarized light observation with previous observations in total intensity (Uyama et al. 2020), and leveraging the unique capability of CHARIS to simultaneously take polarimetry data at all wavelength channels (22 wavelengths across JHK), we examine the surface brightness and polarized fraction as a function of wavelength and scattering angle.We see that shorter wavelengths are brighter in both total intensity and polarized intensity, where the brightness drops by a factor of ∼ 2 from 1.2 to 2.1.The degree of polarization is higher at longer wavelengths, where the polarized fraction peaks around 0.6 in J band, 0.7 in H band, and 0.8 in K band.The maximum degrees of polarization in J and H are both higher compared to Monnier et al. (2019), where they found peak fractions of around 0.5 and 0.6 in J and H bands, respectively.Compared to other planet-forming disks, HD 34700A already has one of the highest degrees of polarization based on the measurements in Monnier et al. (2019), as compiled in Tazaki & Dominik (2022).Our measurements would shift it slightly higher, bringing HD 34700A to a very similar level to the disk, UX Tau A (Tanii et al. 2012).We also detect an inner arc in our data, which is part of the mis-aligned inner disk of the system Columba et al. (2024).
To understand what the observations imply about the dust population in the outer disk of HD 34700 A, we carry out radiative transfer modeling using the 3D Monte Carlo radiative transfer code, MCFOST (Pinte et al. 2006).We make only qualitative comparisons between our models and the data, due to several limitations: 1.We lack knowledge of the inner disk properties, which could cast shadows and re-scatter light from the stars onto the outer disk, and alter its surface brightness.2. Optical properties of dust grains are calculated using Mie theory instead of more exact methods that are more expensive for producing imaging simulations at many wavelengths.3. The disk is simplified as a circular, single-component, axisymmetric ring with a flared surface.Our models struggle to produce a close match to the total and polarized intensity, which is not surprising given the aforementioned limitations.Even for optically thin disks, Mie theory calculations following Bruggeman mixing rules have been shown to struggle to simultaneously match total intensity and polarized fraction Arriaga et al. (2020).Instead, we focus our comparison on the polarized fraction as a function of scattering angle, which only depends on the scattering properties of the dust grain populations and the scattering geometry.We find that a two-population distribution with a broken number density power-law is necessary to reproduce the surface brightness of the data, where significant mass fractions in both sub-micron grains and micron-sized large grains are required.Such a distribution appears somewhat unphysical.This could be because Mie models are using an extra population of small grains to simulate the effects of sub-micron structures and monomers of large dust aggregates.This is consistent with the findings in Tazaki & Dominik (2022), where they investigated the degree of polarization of more complex dust aggregates.They found that micron-sized dust aggregates consisting of 0.1 and 0.2 monomers can have polarized fractions beyond 0.6 at near-infrared wavelengths, consistent with our models that contain a large population of sub-micron dust particles.
We also run models with DHS grains and Mie grains containing amorphous carbon.While these two grain types struggle a bit more to produce the surface brightness profiles, we also find models that match the polarized fraction of the data within the parameter space we explored.We find that DHS particles share similar results as Mie particles, and interestingly, grains containing amorphous carbon require much less porosity to produce a high level of degree of polarization.
For future work, it would be intriguing to employ more detailed grain scattering models of individual monomers and larger aggregates (e.g., Purcell & Pennypacker 1973;Draine & Flatau 1994;Mackowski & Mishchenko 1996;Tazaki & Dominik 2022;Tazaki et al. 2023).This would allow for a deeper investigation into whether realistic grains with rough surfaces and intricate aggregate structures can effectively explain the observed broken power-law distribution and provide a better constraint on the grain sizes.Another welcome improvement to the disk model would be to include the inner disk.This would be possible, for example, at wavelengths corresponding to H  filters (∼ 0.6) where the inner disk has been fully resolved (Columba et al. 2024).Finally, an excess in the spectral energy distribution (SED) was observed at ∼ 3 − 10, which suggested a potential warm disk component very close to the central binary (Monnier et al. 2019).Therefore, an additional source of incident light that accounts for such a component to match the SED would be another welcome improvement.

Figure 1 .
Figure 1.A comparison of the standard RDI with polarimetry-constrained RDI (PCRDI).The flux is in arbitrary units, scaled by the deprojected radius of the disk surface for better visual clarity.PCRDI is shown to visibly improve the over-subtraction of the disk signal and of the background.

Figure 2 .
Figure 2. The Stokes Q (first column), U (second column) and radial Stokes   (third column),   (fourth column) images of the disk, collapsed into J (top row), H (middle row), and K bands (bottom row).Astrophysical signal is present in   with an average SNR of around 4 to 5 depending on wavelength, and up to an SNR of 10 at the brightest regions in J band.

Figure 3 .
Figure 3. First rows: collapsed J, H, K band images.Left column: total intensity from Uyama et al. (2020), re-reduced with PCRDI that minimizes over-subtraction.Middle column:   of our polarized intensity data.Right column: polarized fraction,   /, where  is the total intensity, in the bright disk ring region.The bottom row contains false color composite renderings of the images as RGB images (R=K, G=H, B=J), with the maximum flux in J band normalized to the maximum value of blue.The various disk features are annotated in the very first panel, including the cavity (F1), darkening features (F2-F5), an inner arc F6, and an instrumental artifact (F7)

Figure 5
Figure5shows the azimuthal surface brightness profiles in total intensity (top panel),   (middle panel) and polarized fraction (bottom panel).The different colors correspond to the wavelength

Figure 4 .
Figure 4. Ellipse fitting to the disk.The three panels correspond to three different fitting methods: least-squares distance (left), maximum flux (middle) and interpolated maximum flux (right).In the first two cases, the blue points are the pixels selected to calculate distance or flux.For the last case, the blue points are the equally spaced points on the ellipse for interpolation.The green lines are the best-fit ellipses to the blue points in each case.The best-fit values for semi-major axis, inclination and PA agree within 1 between the second and third methods, and they agree with the those of the first method to ∼ 2; these values are reported in Table2.

Figure 5 .
Figure5.The azimuthal surface brightness profiles in total intensity (top panel) and   (middle panel), and   /total intensity (bottom panel), measured along the best-fit ellipse.The PAs refer to position angles relative to the celestial north.The shaded regions correspond to the approximate PAs covered by the labeled features in Figure3except F6, as F6 is not part of the outer disk.

Figure 6 .
Figure6.The normalized apparent scattering phase functions for total intensity  and azimuthal Stokes   are shown in the top and second panels, respectively.The scattering angle is defined as the angle between the incident light from the star and the scattered light to the observer, with forward scattering having angles < 90°.We show two polarized fraction phase functions measured assuming different PAs in the third and fourth panels.PA=92°c orresponds to the PA of the axis of maximum polarization, and PA=74.• 3 corresponds to the PA of the semi-major axis of the best-fit ellipse described in Section 3.1.The shaded regions indicate the approximate scattering angles covered by the labeled features in Figure3except F6, which is not part of the outer disk.

Figure 7 .
Figure 7.Comparison of models to the data.The azimuthal surface brightness profiles in total intensity and   are shown in the left and middle columns, respectively.The polarized fraction as a function of scattering angle is shown in the right column.The solid lines are the measurements for the PSF-convolved model images, and the data points with error bars are measurements for our datasets using the scattering surface geometry of the corresponding models.Each row shows a different set of parameter models from the model grid, with key parameters annotated in the right column plots.We find no model that can simultaneously fit the total intensity and polarized intensity, which is not surprising considering the significantly simplified model compared to the feature-rich disk, in addition to the effect of the inner disk that we ignored.However, as the polarized fraction reflects the local properties of the dust particles responsible for the scattering, we are able to find models that fit the polarized fraction and they suggest significant masses in sub-micron as well as micron sized grains.

Table 2 .
Best-fit ellipse parameters for different fitting methods

Table 3 .
Monnier et al. (2019);he transition grain size is the maximum size of the small grain population and also the minimum size of the large grain population.The two populations share the same grain size power-law distribution,    References: 1.Monnier et al. (2019); 2. Gaia Collaboration et al. (2021).