© The Authors 2024

1. Introduction

The study of dust properties in protoplanetary disks is crucial for understanding planet formation (Birnstiel 2024; Ginski et al. 2023; Andrews et al. 2018). Multi-band polarimetric imaging offers a powerful tool to probe these properties, as it reveals the scattering behavior of dust grains in the disk’s surface layers. Scattered light properties, such as disk reflectivity and colors (Krist et al. 2005; Debes et al. 2013; Mulders et al. 2013; Tazaki et al. 2019; Hunziker et al. 2021; Ma et al. 2023, 2024), as well as scattering phase functions (Ginski et al. 2023; Tazaki et al. 2023) have been used to constrain the dust properties in the disk surface region. However, our ability to understand the dust properties via measurements from multi-wavelength polarised and/or total intensity observations often relies on the premise of stable, uniform illumination of the outer disk over different observation epochs. Variability in scattered light, particularly when it is due to shadowing effects, can significantly impact our observations and interpretation of dust properties. This requires accurately calibrated polarized reflectivities from multi-epoch imaging polarimetry. These measurements have been scarcely available for protoplanetary disks until now.

PDS 70, a photometrically variable young K7 T-Tauri star of 5.4 Myr located 113.4 pc away in the Upper Centaurus-Lupus subgroup of the Sco-Cen star-forming region (Müller et al. 2018; Pecaut & Mamajek 2016; Gaia Collaboration 2018), serves as a key object in this context. The system’s well-defined inner and outer disks, separated by a clear gap, make it an ideal target for studying the influence of variable illumination on dust properties. The presence of two planets within the gap (Keppler et al. 2018; Haffert et al. 2019; Benisty et al. 2021) further adds to its significant role as a laboratory for planet formation. The spectral energy distribution (SED) fitting initially indicated the existence of these disks (Metchev et al. 2004; Riaud et al. 2006), with direct imaging in scattered and polarized light later resolving the outer disk and the large gap. The surface brightness of the ring-shaped outer disk peaks at 54 au, with an inclination of i = 49.7° and a position angle of 68.6° (Hashimoto et al. 2012; Keppler et al. 2018). The inner disk, although it remains unresolved, is estimated to have a radius of less than 18 au (Keppler et al. 2018; Benisty et al. 2021). Comparison between the dust continuum emission maps and reflected images reveals a similar morphology with dust size segregation (Long et al. 2018; Portilla-Revelo et al. 2022).

Variability is seen in the PDS 70 system across various wavelengths and observing techniques. The optical photometry shows both short-term variations in days and long-term on a yearly scale. The day-scale variation evolves between three-day periodic and stochastic “dipper-like” behavior, attributed to aperiodic occultation by micron-sized dust superposed on periodic stellar rotation (Thanathibodee et al. 2020; Gaidos et al. 2024). Infrared(IR) photometry dominated by inner disk emission also shows short and long-term variations, together with a one-year disappearance of emission at 3.4 μm and 4.6 μm in WISE data (Gaidos et al. 2024). In addition, strong mid-IR variation of 25% is revealed when comparing JWST and Spitzer SED with a 15-year interval (Perotti et al. 2023; Jang et al. 2024), and a comparable level of near-IR variability is seen in comparing JWST and SpeX (Christiaens et al. 2024), pointing to a dynamic inner disk geometry due to variable accretion. The scenario is supported by variable line emissions across multi-epoch (Thanathibodee et al. 2020; Campbell-White et al. 2023; Gaidos et al. 2024). Moreover, multi-epoch ALMA dust continuum emission also suggests variable inner disk structure (Casassus & Cárcamo 2022). The changing inner disk casts variable shadows on the outer disk, seen as azimuthally variable surface brightness in the scattered light and polarized light (Juillard et al. 2022; Christiaens et al. 2024).

In this letter, we combine multi-epoch polarimetric imaging observations of PDS 70, revealing time variability in polarized reflectivity and azimuthal brightness profiles across 0.6–2.2 μm and exploring the potential chromatic dependence of the near-to-far side brightness contrast. Our work underscores the need for the careful interpretation of the measurements of scattered light.

2. Observations and data reduction

For the first time in the narrow band N_R filter (λ = 645.9 nm, Δλ = 56.7 nm), PDS 70 was observed on March 1, 2023. This observation, using the SPHERE instrument on UT3 (Beuzit et al. 2019), is part of the optical scattered light survey (ID: 10.240P.001, PI: C. Ginski), aimed at characterizing the dust properties in young planet-forming disks. The observing sequence is designed to switch between long-exposures and short-exposures on target, combined with star-hopping to the reference star for long-exposures (Wahhaj et al. 2021). In this study, we selected all 11 polarimetric cycles of long-exposures, totalling up to 47 min of integration on target. The observation was performed in instrument-stabilized (P1) and the slow read-out mode (Schmid et al. 2018).

We retrieved the archival data for PDS 70 observed with the ZIMPOL and IRDIS subsystem on the SPHERE instruments in polarimetric differential imaging mode. The observing condition and the data selections are summarized in Table A.1 and detailed in Appendix A.1. To summarize, we collected seven epochs of observations spanning seven years, with five color filters (N_R, VBB, J, H, and Ks) from 0.6 to 2.2 μm. For simplicity, we labeled them by their band filter name, adding the epoch when two observations exist for the same band. Each epoch of observation contains full polarimetric cycles and non-saturated total intensity frames. For non-saturated polarimetric observations NR, VBB, and J17, the non-saturated total intensity is obtained simultaneously. For coronagraphic polarimetric observations J16 and Ks, two separate flux frames were taken before and after the polarimetric cycles. For coronagraphic H19, only one flux frame was taken after all polarimetric cycles. For the saturated polarimetric observation H21, eight pairs of non-saturated Q⁺ and Q⁻ frames are taken after the full polarimetric cycles.

The ZIMPOL and IRDIS data are reduced with their dedicated pipelines sz-software (Schmid et al. 2018) and IRDIS Data reduction for Accurate Polarimetry (IRDAP, van Holstein et al. 2020), respectively. On top of the processed data for each cycle, we conducted an extra step to minimize the integrated U_φ-signal within the disk region. This step should correct the telescope polarization and intrinsic stellar polarization for ZIMPOL and correct the stellar polarization for IRDIS. The disk region is defined as an elliptical aperture with r = [0.3″, 1.0″], with the inclination i = 49.7° and position angle PA = 68.6° adopted from Hashimoto et al. (2012), shown as the gray dashed lines in J16 panel of Fig. 1. The estimated and subtracted polarization are summarized in Table A.2. As the last step, we calculate the azimuthal polarization parameters Q_φ and U_φ from the reduced Stocks polarizations Q and U:

$$\begin{array}{cc}\hfill Q_{\phi }& =-Qcos(2\phi )-Usin(2\phi ),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill U_{\phi }& =Qsin(2\phi )-Ucos(2\phi ),\hfill \end{array}$$(2)

Fig. 1. Mean polarized intensity images Qφ(x, y) for N_R, VBB, J, H, and Ks band. All images are normalized to the stellar intensity I⋆ and divided by their pixel areas (3.6 × 3.6 mas2 for ZIMPOL and 12.25 × 12.25 mas2 for IRDIS) to provide the surface brightness contrast in mas−2, as indicated by the color bar. The first and the second rows present coronagraphic and non-coronagraphic observations respectively. The shared color bar is indicated on the right. The epochs are indicated on the top left in each panel. The grey dashed line in the J-band panel indicates the integrated region for disk signals. The black dotted line in the H panel indicates the deprojected region used to calculate the azimuthal profile. Up is to the north and left is to the east.

where φ is the position angle to the central star from north over east. The mean Q_φ(x, y) for each epoch are plotted in Fig. 1 and the corresponding U_φ(x, y) are plotted in Fig. A.1.

3. Results

Figure 1 shows temporal and chromatic variation of disk scattered light around PDS 70. Starting with the Ks image, which has the appearance of an archetypal transition disk, the disk near side is brighter than the far side, while it is nearly mirror symmetric to the minor axis. The near to far side brightness contrast is usually attributed to preferential forward scattering by grains large compared to the observing wavelength (e.g., Duchêne et al. 2004; Mulders et al. 2013. However, this typical appearance is absent in other datasets. In the H bands (H19 and H21), the disk is not mirror symmetric with its southern ansa being brighter. Additionally, comparing H19 and H21, the disk flux relative to the stellar flux varies significantly, hinting at temporal variations in the system. More strikingly, at optical wavelengths (NR and VBB), the far side flux dominates over the near side. Consequently, the common picture of the bright side being the near side of the disk does not hold in this case, perhaps because the disk surface on the near side is at least partly hidden, since the illuminated inner disk wall is too steep to be visible (Ma & Schmid 2022). In the following sections, we will further quantify and discuss these temporal and chromatic variations.

3.1. Integrated polarized reflectivity

As the first step, we integrated the total intensity I_⋆ = ∫I(x, y) within a 3″-diameter round aperture, to monitor the short-term stellar flux variability and calibrate the polarized disk flux. No significant temporal change in flux is observed for a given observing run lasting typically one hour. The integrated polarized flux Q_φ is measured from the polarized intensity images Q_φ(x, y) cycle by cycle within the elliptical aperture shown in Fig. 1 J16 panel. Table A.3 summarizes the integrated system flux I_⋆, disk polarized flux Q_φ, and their ratio and the detailed derivation is available in Appendix A.3.

For accurate photo-polarimetric measurements, PSF smearing effects must be considered. The PSF smearing is crucial because it causes a blurring of the intrinsic disk signal and introduces polarimetric cancellation between regions with positive and negative Q and U signals (Schmid et al. 2006; Tschudi & Schmid 2021). We followed the correction method described in Ma et al. (2024) to calculate the intrinsic disk-integrated polarized flux. Figure 2 shows the corrected disk-integrated polarized reflectivity ${\widehat{Q}}_{\phi }/I_{\star }$ as a function of wavelength and the values are summarized in Table A.3. Our results are consistent with the ${\widehat{Q}}_{\phi }/I_{\star }$ measurements of VBB, J17, H19, Ks data in Ma et al. (2024) considering we used a smaller integration region.

Fig. 2. Corrected integrated polarized reflectivity ${\widehat{Q}}_{\phi }/I_{\star }$ as a function of wavelengths. The solid and hollow markers represent the values derived from non-coronagraphic and coronagraphic observations, respectively. The horizontal errorbars indicate the width of the band filters. The observing epochs are labeled beside each scatter.

At given wavelength, the ${\widehat{Q}}_{\phi }/I_{\star }$ measurements of J16 and J17 are consistent. In contrast, ${\widehat{Q}}_{\phi }/I_{\star }$ is about 30% higher in H21 than in H19, and 140% higher in VBB than in NR, not explained by the smearing. This inconsistency is also not caused by the coronagraph, as we would expect similar behavior in the J band. The likely explanation for this discrepancy is a structural change in the inner disk, changing the overall illumination on the outer disk, as discussed in Sect. 4.2.

3.2. Azimuthal profiles of polarized light

For all seven observations, we derived the azimuthal profiles by averaging Q_φ(x, y) over an annulus with a deprojected radius of 35–70 au in 10° bins. The deprojected radius is calculated using diskmap (Stolker et al. 2016b), assuming a power-law function of the disk surface h(r) = 0.041 au⋅(r/1 au)^1.25 (Ginski et al. 2023). The annulus is shown in Fig. 1 H19 panel. The azimuthal profiles are plotted in Fig. 3, starting from the near side PA = 246°. The profiles are normalized to the peak value at given wavelengths and shifted for clarity. The NR profile is normalized to the peak value at VBB without offset. The J-(J16, J17), H-(H19, H21), and Ks-band profiles are offset by 1, 2, and 3, respectively. The J profiles use the J17 peak for normalization, and H profiles use the H21 peak. The vertical dashed lines indicate the minor and major axis. Our Q_φ(φ) profiles for VBB and J17 agree with the analysis of Keppler et al. (2018) across all position angles within the error margin for VBB and most angles for J17 except the near side. The comparisons are detailed in Appendix A.2.

Fig. 3. Azimuthal profile of the Qφ(φ) images averaged over an annulus of 35–70 au. The colored dots are calculated from the observations, with the uncertainties being the standard deviation over all selected cycles. The overplotted curves are calculated from the convolved fitted model. The vertical dashed lines indicate the position angles of the minor axis, smin(W) and smin(E), in the west and east directions, and the major axis, smaj(N) and smaj(S), in the north and south directions.

Comparing the two azimuthal profiles at J band, J16 and J17 marginally agree at most azimuthal angles within the errorbars. Except the J17 lacks the flux at PA ∼ 200°, likely a local shadow, which leads to the integrated flux difference. H21 agrees with H19 only near PA ∼ 130° and 250°, while for all other azimuthal directions, H21 is brighter than H19, with the surface brightness maximas being approximately 80% and 50% brighter on the N and S sides, respectively. Such an overall brightness change also exists between NR from 2023 and VBB from 2015. The surface brightness maximas in VBB are about 120% and 80% brighter on the N and S sides than in NR. Furthermore, when comparing the azimuthal profiles across all wavelengths, we observe a tentative trend that the position of the surface brightness maxima shifts closer to the major axis for longer wavelengths, though this trend is partly mixed with time variability.

The PSF smearing effects also affect the azimuthal variance Q_φ(φ) and need to be corrected. We describe the correction in Appendix B with an analytical disk model. The azimuthal profiles calculated from the convolved best-fitting model are overplotted in Fig. 3. To conclude, the overall azimuthal profile shape is preserved after the convolution and the PSF smearing effect is not substantial, especially when comparing the contrast between position angles symmetric to the minor axis. We introduce ${\widehat{\alpha }}_{X-Y}$ to quantify the brightness contrast between two azimuthal directions X and Y,

$$\begin{array}{c}\hfill {\widehat{\alpha }}_{X-Y}=\frac{{\widehat{Q}}_{\phi }(X)-{\widehat{Q}}_{\phi }(Y)}{{\widehat{Q}}_{\phi }(X)+{\widehat{Q}}_{\phi }(Y)},\end{array}$$(3)

with ${\widehat{Q}}_{\phi }(X)$ being the averaged corrected polarized flux at X direction in a 20° wedge with deprojected radius 35-70 au.

First, we measured the contrast between the west (W, near side, PA = 248.6°) and east (E, far side, PA = 68.6°) directions, denoted as ${\widehat{\alpha }}_{W-E}$. The results are plotted in Fig. 4 upper panel and the values are given in Table B.2. This contrast tentatively increases with wavelength, exhibiting a faint near side in the visible, which becomes equally bright near the J band, and brighter than the far side at the H and Ks band. We also calculated α_W − E without smearing effect correction for comparison with Keppler et al. (2018). Our VBB α_W − E = −0.45 ± 0.17 aligns closely with Keppler et al. (2018) measurement of −0.51. However, we report for J17 α_W − E = 0.11 ± 0.13, suggesting a slightly brighter W-side, whereas Keppler et al. (2018) reported −0.17. This small discrepancy is probably due to different normalization methods and uncertainty definitions. Despite this, the tentative trend of increasing near-to-far contrast with wavelength still holds. This well-defined effect can be caused by the wavelength dependence of the dust scattering properties or disk scattering geometry.

Fig. 4. Brightness contrast of PDS 70. Upper: Brightness contrast between the W (near side) and the E (far side) along the minor axis. Lower: Brightness contrast between the brightness maximas on the S and the N close to the major axis. The errorbar on wavelength is omitted.

The second contrast was derived between the two brightest regions near the major axis N and S (PA = φ₁ and φ₂ are epoch-dependent, values are given in Table B.1), denoted as ${\widehat{\alpha }}_{S-N}$. The results are plotted in Fig. 4 lower panel. This contrast shows no systematic trend with wavelength. Ideally, a geometrically symmetric disk expects ${\widehat{\alpha }}_{S-N}=0$ across all wavelengths. While equal brightness is seen for VBB and Ks, the W side is brighter for NR, J16, J17, and the E side is brighter for H19 and H21. Moreover, J16 and J17 have similar contrasts of ${\widehat{\alpha }}_{S-N}=-0.18\pm 0.08$ and −0.10 ± 0.08, respectively. And H19 and H21 have the same level of opposite contrast ${\widehat{\alpha }}_{S-N}=0.24\pm 0.04$ and 0.25 ± 0.05 respectively. Therefore, despite the opposite brighter side comparing J and H bands, ${\widehat{\alpha }}_{S-N}$ remains consistent at each given wavelength, without systematic trend over a time scale of a few years.

4. Discussion

4.1. Time variability

Our multi-epoch polarimetric data reveal notable integrated polarized reflectivity variability and azimuthal brightness variability in the polarized scattered light from the PDS 70 outer disk. This finding is based on calibrated imaging polarimetry taking the PSF smearing effects carefully into account. This variability can be linked with the recently reported light curve variability by Gaidos et al. (2024), which evolves between periodic and stochastic behavior, hinting at transitions between stable and unstable accretion regimes (Gaidos et al. 2024).

In the stable accretion regime, the matter flows along the magnetic field as funnel streams that are lifted above the equatorial plane. This configuration creates broad and quasi-periodic shadows on the outer disk, as the elevated streams intercept a significant portion of the stellar light before it reaches the outer regions. Consequently, the outer disk appears fainter in scattered light due to the reduced illumination. While the periodic shadows would result in predictable patterns in the azimuthal brightness variability, our observations with the shortest intervals of 6 months, are not sensitive to the quasi-periodic shadowing of 3 days (Thanathibodee et al. 2020).

On the other hand, during the unstable accretion regime, the system is dominated by Rayleigh-Taylor instabilities that drive matter through equatorial tongues (Romanova et al. 2012). Since these tongues are closer to the disk plane, they intercept less stellar light than the elevated funnels and more light reaches the outer disk, resulting in a brighter outer disk. The chaotic accretion and the thin tongue structures cause random and narrow shadows on the outer disk. This equatorial accretion alone is unlikely to obscure the star enough to explain the observed variability. Instead, a magnetized wind, as suggested by Gaidos et al. (2024), may also contribute to the dimming and variability.

The switching is not only supported by the optical light curve, but also by significant IR variations observed over different time intervals. For instance, near-IR measurements have shown a variation of 25% when comparing data taken years apart, indicating dynamic changes in the inner disk due to variable accretion (Christiaens et al. 2024). Additionally, mid-IR variability comparisons between JWST and Spitzer also support this dynamic inner disk geometry (Perotti et al. 2023). However, mapping our observations to the variability is not possible due to the temporal gaps and the lack of high-cadence photometry during our observations.

In the meantime, stellar flux variability also contributes to the observed integrated flux variation thought not directly recognizable in our observation. Repeated stellar flux measurements taken during our single observing epochs with a duration of about one hour remain stable at σ(I_⋆)/I_⋆ ≈ 1% and exhibit no systematic flux changes. This uncertainty falls short of resolving the periodic variation typically at 0.5% per hour. Although we can exclude the occurrence of deep dimming events (Δg> 0.3^m in 1 day) during our observation, shallow dimming may still be present. Gaidos et al. (2024) also reported long-term behavior in the optical, showing an upper envelope of brightness oscillating on a yearly timescale by a few tenths of a magnitude, which could add variability.

Additionally, part of the azimuthal variance can be attributed to data reduction effects. For example, Ma et al. (2023) demonstrated that for RX J1604, applying different intrinsic polarization subtraction can shift the broad shadowing features but does not contribute the integrated flux. Indeed, comparisons of our profiles with that of Keppler et al. (2018) reveal slight discrepancies (see Appendix A.2). However, the general agreement across most position angles suggests that data reduction effects are minor for PDS 70 and the azimuthal brightness variability is primarily influenced by shadowing effects. Lastly, the impact of multiple scattering on the large-scale azimuthal dependence is also minimal (see the U_φ images in Appendix A.2).

The time variability complicates the disk color measurement. Considering the higher fluxes as the unshadowed value at visible and H band, PDS 70 shows a gray color from 0.6 to 2.2 μm. However, using NR data would indicate a red color, and H19 data would introduce an abrupt flux drop. The variability also complicates the extraction of scattering phase functions (Ginski et al. 2023). Time-dependent changes in the phase functions could lead to incorrect interpretation of dust grain size, composition, or distribution if not properly accounted for. Therefore, the measurements of reflectivity colors and scattering phase functions must be interpreted with caution for PDS 70 and other systems showing significant brightness variation for the reflected light from the disk.

To explore the time variability and structural changes in the inner disk more thoroughly, continuous photometric monitoring in both visible and IR wavelengths, paired with simultaneous high-contrast imaging and/or interferometry, will be essential. These approaches will allow us to determine the accretion regime, identify shadowing effects, and capture the structure of the inner disk, providing deeper insights into the mechanisms driving the observed variability.

4.2. Wavelength dependence of brightness contrast

Despite variations in integrated flux and brightness profile, the brightness contrasts – both between near and far and between brightness maxima at given wavelengths – appear relatively stable over time based on the available observations. Here, we discuss some possibilities causing the observed tentative chromatic dependence from optical to near-IR.

The first possible cause for the wavelength dependence of the near-to-far side asymmetry is the polarized intensity scattering phase function of dust grains in the disk surface. Polarized intensity is a product of the degree of polarization and total intensity; thus, the polarized intensity data alone cannot inform us of which specific factor is creating this wavelength dependence. Assuming the degree of polarization is wavelength-independent, grains have to be more efficient backscatterers at optical compared to near-IR. For example, by assuming the Heyney-Greenstein phase function and a plane-parallel radiative transfer model developed in Ma & Schmid (2022), we estimated the asymmetry parameter of anisotropic scattering (Bohren et al. 1983) to be around g ≈ −0.25 in the optical to g ≈ 0.1 in the Ks band (see Appendix C for details). Backward-scattering particles is also suggested by Wahhaj et al. (2024) from PDS 70’s YJH total intensity images. Such a peculiar backward scattering behavior is not typically expected for small dielectric particles (e.g., Rayleigh scattering; Bohren et al. 1983) or large fluffy dust aggregates composed of sub-micron monomers (Volten et al. 2007; Tazaki et al. 2016). One possible explanation is the presence of small metallic dust grains, which exhibit strong backward scattering (van de Hulst 1957; Jackson 1975). However, this would predict a similar degree of backscattering at longer wavelengths, which is not observed, making small metallic particles unlikely. Another intriguing possibility is the presence of large compact aggregates in the disk surface, which might also account for the enhanced backward-scattering behavior. Min et al. (2016) demonstrated that compact dust aggregates of 4 μm in size could exhibit a negative g value if the observer only has access to intermediate scattering angles. In their simulations, a negative g value does not appear for aggregates smaller than 2 μm (see Figs. 4 and 5). This increase in backward scattering is thought to result from the diffusive reflection on the rough dust surfaces, such as the lunar phase (Mukai et al. 1982; Min et al. 2010). Therefore, both large size and some degree of compactness of aggregates seem crucial for this phenomenon. However, quantitative conditions for the dust properties remain unclear and require future studies. If this explanation accounts for the observed negative g value, the surface region of PDS 70 must be dominated by large compact grains or aggregates at least a few microns in size. On the other hand, the degree of polarization may change drastically with wavelength (Tazaki & Dominik 2022). A detailed investigation of such dust properties is beyond the scope of this paper, but observations of the degree of polarization are crucial to resolving parameter degeneracy in polarized intensity.

The second possibility is obscuration of the near-side scattered light by the inner edge of the outer disk. The scattering surface is expected to be at a progressively higher altitude for shorter wavelengths due to increased opacity, particularly notable for small grains, which can eventually cause self-obscuration of the near-side scattered light (Ma & Schmid 2022). Based on the plane-parallel radiative transfer model, we found that the wall slope needs to be > 90 ° −i = 40° (see Appendix C). Notably, this steep slope is consistent with the vertical gas structure constraints from recent CO isotopologue observations, supporting the feasibility of this kind of geometry (Law et al. 2024).

The third possibility is non-uniform outer disk illumination due to a misaligned inner disk. Misaligned disks have been suggested in various systems based on scattered light observations, (e.g., Marino et al. 2015; Stolker et al. 2016a; Benisty et al. 2017; Pinilla et al. 2018; Muro-Arena et al. 2020. PDS70 hosts an inner disk (Keppler et al. 2019), which potentially casts a shadow onto the outer disk. This shadowing effect is likely more efficient at shorter wavelengths (optical) due to higher opacity, while being partially transparent in the near-IR. However, in this case, the inner disk’s tilt has to be fine-tuned so that the shadow covers exactly the near side of the outer disk. The stellar light passing through the atmospheric layer of the inner disk before reaching the outer disk could impose an additional polarization signal. The current models of PDS 70 typically assume an aligned or even absent inner disk and often overlook time variations in scattered light. These simplifications could lead to an incomplete understanding of the disk’s geometry and the scattering properties of dust. Accounting for a misaligned inner disk and the time variability would be beneficial for developing more realistic models and interpretations of observed scattered light.

Stability is also observed in the brightness contrast between two maxima, ${\widehat{\alpha }}_{S-N}$. Consistent with our measurements, PDS 70’s HiCIAO PDI data from 2012 reported α_S − N = 0.13 ± 0.06 for H-band (Hashimoto et al. 2015). Meanwhile, SMA and ALMA observations of dust continuum emission indicate an azimuthal brightness enhancement on the N side at 0.87 mm (Long et al. 2018; Keppler et al. 2019) and 1.3 mm (Hashimoto et al. 2015). If the asymmetric α_S − N(λ) is a time-independent feature, Hashimoto et al. (2015) proposed a long-lived dust clump induced by accreting planets to explain the opposite asymmetry between H-band scattered light and ALMA dust emission. However, the consistency of J band scattered light asymmetry with ALMA contradicts the explanation. These tentative time-independent features need further observations to confirm, especially simultaneous multi-wavelength observations to mitigate the impact of temporal variability and properly interpret the wavelength dependence in scattered light observations. Simultaneous total intensity observations will also help distinguish between shadowing effects and intrinsic dust property variations in the faint regions.

5. Conclusions

We present the first polarized reflected image of PDS 70 at the N_R band, which is also the first high quality AO-data for λ < 1 μm taken under good seeing conditions. Combining the archival polarimetric observations across 0.6 to 2.2 μm, we report the variability in both integrated polarized reflectivity and the azimuthal brightness profiles. We found evidence of significant time variability in the polarized scattered light at optical wavelengths (NR and VBB), while the azimuthal brightness profiles remain stable between the two epochs (Figs. 2 and 3). In the near-IR wavelengths (J to Ks), we observed a hint of variability in the H band (H19 vs. H21); however, the current inhomogeneous datasets prevent us from drawing robust conclusions.

Additionally, we found a consistent trend of a brighter near side relative to the far side with wavelength regardless of observing epochs (Fig. 4, top). Interestingly, this near side disappears entirely at visible wavelengths. Furthermore, an excess brightness bump in the southeast ansa in the H band data is absent at other wavelengths (Fig. 4, bottom).

Given the variability of PDS70 and its significant influence on the outer disk scattered light, it is crucial to conduct simultaneous observations across the optical to near-IR wavelengths to extract the disk and dust properties of PDS70. Our results highlight a new complexity involved in analyzing scattered light observations made at different times around variable stars.

Acknowledgments

We thank the referee for the constructive comments, which have significantly improved the quality and clarity of this letter. We also thank Eric Gaidos for the discussion on photometric variability. JM thanks the Swiss National Science Foundation for financial support under grant number P500PT_222298. RT and FM acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation program (grant agreement No. 101053020, project Dust2Planets.)

References

Appendix A: Details on data selection and calibrations

A.1. Data selection

Table A.1 summarizes the selected dataset and the corresponding observing conditions. The reasons for discarding some science cycles are given in the footnotes.

A.2. Polarimetric calibration

ZIMPOL pipeline does not correct the telescope polarization and IRDAP corrects the instrumental polarization and cross-talks (van Holstein et al. 2020). For ZIMPOL observations we also estimated and corrected for the telescope polarization which depends on the parallactic angle Schmid et al. (2018). The interstellar polarization and the intrinsic polarization contribution estimation and removal are optional for both pipelines by normalizing the integrated Q and U in the given region. Given that our datasets cover different seeing conditions and include coronagraphic observations, the normalization region is not easy to define uniformly. Therefore, we took the non-normalized data and minimized the integrated absolute U_φ in the disk region to all datasets to correct for the interstellar and intrinsic polarization of the central object. For each cycle of polarimetric observation, we calculated Q′(x, y) = Q(x, y)−q_⋆ ⋅ I_Q(x, y), U′(x, y) = U(x, y)−u_⋆ ⋅ I_U(x, y) and computed U′_φ(x, y) with Eq. 2. We searched for the q_⋆ and u_⋆ to minimize the integrated abs(U′_φ(x, y)) in the disk region. The average and the standard deviation among all cycles for each epoch are summarized in Table A.2.

The mean U_φ images for each epoch are shown in Fig. A.1 corresponding to the Q_φ images in Fig. 1. The same color scale is used to show that the U_φ signals are generally weak for the outer disk. The rather strong positive and negative U_φ signal are expected for the inner disk because of PSF convolution effects and the contamination by the intrinsic polarization signal from the unresolved central source. The U_φ signals have also been reported by van Holstein et al. (2021) and explained by multiple scattering in the outer disk. However, the expected U_φ from multiple scattering (Canovas et al. 2015) is weaker than what we observe for the inner disk and what was reported by van Holstein et al. (2021).

Fig. A.1. Stokes Uφ images. All images are normalized in the same way as Fig. 1. The labels and color scales are also the same. The dashed lines, initially used to calculate the azimuthal brightness variation in the H19 panel of Fig. 1, serve to highlight the corresponding area of the disk here.

Our estimations of stellar polarization are in agreement with the previous studies in Keppler et al. (2018), Ma et al. (2024), although those works normalized the data within the central region. Pecaut & Mamajek (2016) measured low optical extinction for PDS 70, so the p_⋆ is mostly due to intrinsic polarization instead of interstellar polarization. We noticed that p_⋆ for H19 and H21 differs, hinting at variable polarization from the unresolved inner disk. Although different normalization step is applied, the effect of data reduction on the azimuthal profile for PDS 70 is minor. To facilitate a clearer comparison, we extracted the curves from Fig 3 in Keppler et al. (2018), shifted the starting position angle to match our convention, and overplotted our profiles for epoch J17 and VBB in Fig. A.2 As shown in the figure, the profiles agree with each other at most position angle within our errorbars. We used the standard deviation among all cycles as the uncertainty in our measurements, reflecting the variability observed across cycles. The typical uncertainty is 10-20% on each point.

Fig. A.2. Comparison of azimuthal profiles Qφ(φ)between this work and Keppler et al. (2018). The J17 profiles are offset by 1. The vertical dashed lines follow the same convention as in Fig. 3

A.3. Photometric calibration

In this subsection, we detail the integration of stellar flux, I_⋆, the polarized disk flux, Q_φ, and their ratio. For each non-saturated total intensity image I(x, y), we estimated and subtracted the background flux, so that the radial profile of the flux falls as a power-law I(r)∝r^α similar to the PSF halo. The stellar flux I_⋆ is integrated within the 3″ diameter round aperture. The uncertainty is the standard deviation among all frames of a given observing run. The I_⋆ for NR, VBB, and J17 are simultaneously obtained with the polarimetric data, allowing monitoring of I_⋆ during the whole observation and showing very stable σ(I_⋆)/⟨I_⋆⟩ = 0.6%,1.1% and 0.3%, respectively. For the H21 data, eight I_Q images show a scatter of σ(I_⋆)/⟨I_⋆⟩ = 5% due to noise at the large separations. No systematic flux increase or decrease with time is seen in the monitoring. The J16 and Ks data have only two flux frames, and as the two measurements agree with each other, we adopt σ(I_⋆)/⟨I_⋆⟩ = 3%. For H19 which has only one flux frame, we adopt σ(I_⋆)/⟨I_⋆⟩ = 5%. We note that the I_⋆ value includes disk contribution and, therefore overestimates the stellar intensity. The hot dust emission contributes up to about 10% at K-band (Wang et al. 2021), but less at shorter wavelengths.

The integrated polarized flux Q_φ are integrated from the polarized intensity images Q_φ(x, y) cycle by cycle within the elliptical aperture shown in Fig. 1 J16 panel. The uncertainty is the standard deviation of all cycles. Typical uncertainty σ(Q_φ)/⟨Q_φ⟩ = 7% does not show systematic trends. High uncertainty reaches 15% at N_R and J, mainly due to the low signal-to-noise ratio (S/N). For NR, VBB, and J17, simultaneous Q_φ and I_⋆ measurements are available, the Q_φ/I_⋆ and the uncertainty are the mean and standard deviation of all cycles. For non-simultaneous I_⋆, Q_φ/I_⋆ is the division of ⟨Q_φ⟩ and ⟨I_⋆⟩-values considering the uncertainty based on error propagation.

The same process is applied for the corrected integrated polarized flux ${\widehat{Q}}_{\phi }$. The correction follows the method described in Ma et al. (2024). The method does not require detailed parametric modeling, but only (quasi-) simultaneous PSFs to correct the smearing effect on the integrated flux. The VBB band is the most affected by the PSF smearing due to the poor seeing condition and the observed polarized flux is only 44% of the intrinsic polarized flux (see Table A.3). When applying A-type correction, ${\widehat{Q}}_{\phi }/I_{\star }$ is shown to have a smaller uncertainty than Q_φ/I_⋆ for NR and VBB because simultaneous PSF corrects the smearing effect better. For non-simultaneous I_⋆ and disk polarimetry data the smearing correction only retrieves the lost flux but does not improve the accuracy (B-type correction in Ma et al. (2024)). All the integrated fluxes are summarized in Table A.3.

Appendix B: Correction for azimuthal profile

The correction method developed by Ma et al. (2024) corrects for the smearing effect on the integrated flux but doesn’t address the smearing along the azimuthal direction. To investigate the PSF smearing effects on the azimuthal signal distribution we need a model fitting procedure similar to Ma et al. (2023) which takes the intrinsic azimuthal dependence of the polarization into account. We simulate the asymmetric asymmetric brightness along the disk with analytical formulation A(φ):

$$\begin{array}{c}\hfill A(\phi )=A_0+A_1·exp(-\frac{{(\phi -{\phi }_1)}^2}{2{\sigma }_1^2})+A_2·exp(-\frac{{(\phi -{\phi }_2)}^2}{2{\sigma }_2^2}),\end{array}$$(B.1)

where φ is the position angle measured from north to east, with A₁, A₂, φ₁, φ₂, σ₁, and σ₂ representing the intensities, azimuthal positions, and widths of two brightness peaks, respectively, and A₀ indicating the baseline brightness level.

In the radial direction, the disk is characterized by a double power law function:

$$\begin{array}{c}\hfill Q_{\phi }/I_{\star }(r,\phi )=A(\phi )·{({\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{\mathrm{in}}}+{\left(\frac{r}{r_0}\right)}^{-2{\alpha }_{\mathrm{out}}})}^{-1/2}.\end{array}$$(B.2)

The r value is the deprojected radius using diskmap, r₀ is the approximate peak radius, and the two power-law indices α_in > 0 and α_out < 0 define the inner and outer disk slope, respectively.

The inner edge is not spatially resolved and we fixed α_in= 15 for all wavelengths. The remaining seven parameters r₀, α_out, A₀, A₁, A₂, φ₁, andφ₂ are variable. The analytical model images are convolved with their corresponding PSFs and fitted to the averaged Q_φ/I_⋆(x, y) image for all selected epochs in the deprojected 35-70 au aperture as shown in Fig. 1. The fitted parameters are summarized in Table B.1. We note that the model aims not for a perfect azimuthal profile match, but to assess the relative change of the azimuthal profile caused by PSF convolution.

To assess the relative change, we convolved the best-fit analytical model with the corresponding PSFs and obtained the convolved model images ${\stackrel{\sim }{Q}}_{\phi }(x,y)$. We follow the same procedure to calculate the convolved azimuthal profile ${\stackrel{\sim }{Q}}_{\phi }(\phi )$, depicted in Fig. 3 as lines. The profiles ${\stackrel{\sim }{Q}}_{\phi }(\phi )$ and ${\stackrel{\sim }{Q}}_{\phi }(\phi )$ are normalized to their peak values. The relative difference $\mathrm{\Delta }{Q}_{\phi }(\phi )={\widehat{Q}}_{\phi }(\phi )-{\stackrel{\sim }{Q}}_{\phi }(\phi )$ illustrates the effect of PSF smearing on the azimuthal profile. Figure B.1 displays ΔQ_φ(φ) for each epoch, indicating a trend where brightness along the major axis gets weaker and that along the minor axis gets stronger, consistent with the "smearing" concept. The near side loses slightly more flux than the far side due to the inclined disk geometry, which makes the near side appear closer to the center. In general, the PSF smearing alters brightness by less than 10% across all position angles. Exceptions are the VBB and J band, where the relative difference reaches 20% near the major and minor axis respectively. Thus, the PSF smearing effect on PDS 70 is not substantial and preserves the overall profile shape before and after the convolution.

To correct the smearing effect on the azimuthal profiles, we multiplied the factor $c(\phi )={\widehat{Q}}_{\phi }(\phi )/{\stackrel{\sim }{Q}}_{\phi }(\phi )$ with the observed Q_φ(φ). Then we measured the brightness contrast as described in Sect. 3.2. The corrected brightness contrasts are summarized in Table B.2. Here, ${\widehat{\alpha }}_{S-N}$ doesn’t change much after correction, while ${\widehat{\alpha }}_{W-E}$ is slightly higher after the correction because the near-side (W) is closer to the center and therefore loses more flux.

Fig. B.1. Relative difference, $\mathrm{\Delta }{Q}_{\phi }(\phi )={\widehat{Q}}_{\phi }(\phi )-{\stackrel{\sim }{Q}}_{\phi }(\phi ),$ between the intrinsic azimuthal profile and the convolved azimuthal profile.

Appendix C: Flat disk modeling

This section aims to present flat surface disk results based on the model grid calculation from Ma & Schmid (2022), which could qualitatively reproduce the observed flux and near-to-far contrast as a function of wavelength.

For the disk geometry, we fixed the inclination to i = 47.5°, which is the closest value that the plane-parallel model can accommodate to the i = 49.7° observed in reflected disk images. We also fixed the disk opening angle α = arctan(h/r) = 8°, which follows from the estimated h = 7au and r = 50au derived by ellipse fitting (Takami et al. 2014; Keppler et al. 2018). For the scattering parameters, we fixed the single scattering albedo ω = 0.5 and the maximum polarization p_max = 0.5, as they do not greatly affect the brightness contrast (Ma & Schmid 2022). In the model grid, we varied the asymmetry parameter, g ∈ [ − 0.2, 0.5], with a step of 0.1 and disk wall slope of χ ∈ [10 ° ,50 ° ],  with a step of 5°, creating a coarse grid of disk models. For each model, we calculated the brightness contrast ${\widehat{\alpha }}_{W-E}$, following the same process described in Sect. 3, and the results ${\widehat{\alpha }}_{W-E}(g,\chi )$ are plotted in Fig. C.1.

Fig. C.1. Modelling result assuming plane-parallel surface disk. Left: Intrinsic brightness contrast ${\widehat{\alpha }}_{W-E}$ as function of the asymmetric parameter g and the disk wall slope χ. Levels are indicated with solid lines with values. Right: Examples of simulated polarized disk images, with parameters g = 0.1 and χ = 45° in the upper panel and g = 0.1 and χ = 20° in the lower panel.

To reproduce ${\widehat{\alpha }}_{W-E}<-0.25$ seen in the visible, backward scattering dust at the level g < −0.2 is needed for a flat disk wall χ ≈ 10°; otherwise, a high disk wall χ > 35° is required for small dust g ≈ 0. Even higher χ is needed for forward scattering dust. Figure C.1 presents two examples of flat models, which could resemble the observations. It should be noted that these simple disk models assume a plane-parallel disk surface structure. Therefore, the simulations produce no signal from the disk near side if χ > 90 − i = 40° for PDS 70. For a reliable estimate of the dust asymmetry parameter of a moderately inclined disk like PDS70, we should consider more detailed models with curved disk surface geometry and a vertical stratification for the scattering dust particles.

All Tables

All Figures

Fig. 1. Mean polarized intensity images Qφ(x, y) for N_R, VBB, J, H, and Ks band. All images are normalized to the stellar intensity I⋆ and divided by their pixel areas (3.6 × 3.6 mas2 for ZIMPOL and 12.25 × 12.25 mas2 for IRDIS) to provide the surface brightness contrast in mas−2, as indicated by the color bar. The first and the second rows present coronagraphic and non-coronagraphic observations respectively. The shared color bar is indicated on the right. The epochs are indicated on the top left in each panel. The grey dashed line in the J-band panel indicates the integrated region for disk signals. The black dotted line in the H panel indicates the deprojected region used to calculate the azimuthal profile. Up is to the north and left is to the east.

In the text

	Fig. 2. Corrected integrated polarized reflectivity ${\widehat{Q}}_{\phi }/I_{\star }$ as a function of wavelengths. The solid and hollow markers represent the values derived from non-coronagraphic and coronagraphic observations, respectively. The horizontal errorbars indicate the width of the band filters. The observing epochs are labeled beside each scatter.
In the text

Fig. 3. Azimuthal profile of the Qφ(φ) images averaged over an annulus of 35–70 au. The colored dots are calculated from the observations, with the uncertainties being the standard deviation over all selected cycles. The overplotted curves are calculated from the convolved fitted model. The vertical dashed lines indicate the position angles of the minor axis, smin(W) and smin(E), in the west and east directions, and the major axis, smaj(N) and smaj(S), in the north and south directions.

In the text

	Fig. 4. Brightness contrast of PDS 70. Upper: Brightness contrast between the W (near side) and the E (far side) along the minor axis. Lower: Brightness contrast between the brightness maximas on the S and the N close to the major axis. The errorbar on wavelength is omitted.
In the text

	Fig. A.1. Stokes Uφ images. All images are normalized in the same way as Fig. 1. The labels and color scales are also the same. The dashed lines, initially used to calculate the azimuthal brightness variation in the H19 panel of Fig. 1, serve to highlight the corresponding area of the disk here.
In the text

	Fig. A.2. Comparison of azimuthal profiles Qφ(φ)between this work and Keppler et al. (2018). The J17 profiles are offset by 1. The vertical dashed lines follow the same convention as in Fig. 3
In the text

	Fig. B.1. Relative difference, $\mathrm{\Delta }{Q}_{\phi }(\phi )={\widehat{Q}}_{\phi }(\phi )-{\stackrel{\sim }{Q}}_{\phi }(\phi ),$ between the intrinsic azimuthal profile and the convolved azimuthal profile.
In the text

	Fig. C.1. Modelling result assuming plane-parallel surface disk. Left: Intrinsic brightness contrast ${\widehat{\alpha }}_{W-E}$ as function of the asymmetric parameter g and the disk wall slope χ. Levels are indicated with solid lines with values. Right: Examples of simulated polarized disk images, with parameters g = 0.1 and χ = 45° in the upper panel and g = 0.1 and χ = 20° in the lower panel.
In the text