© The Authors 2026

1. Introduction

Most embedded protoplanet candidates known to date, such as PDS 70b,c, AB Aur b, WISPIT 2b, and 2MASS J16120668-3010270 b, are Hα sources and thought to be undergoing mass accretion (Keppler et al. 2018; Haffert et al. 2019; Wang et al. 2021; Christiaens et al. 2024; Hammond et al. 2025; Currie et al. 2025; van Capelleveen et al. 2025; Li et al. 2025; Close et al. 2025b). Among these candidates, only PDS 70c is known to have a radio counterpart (Isella et al. 2019; Benisty et al. 2021; Fasano et al. 2025; Domínguez-Jamett et al. 2025). The sub-millimetre continuum signal from PDS 70c is particularly important as it constrains the planetary accretion rate and conditions in its environment, provided the radio emission mechanisms are understood. By contrast, Hα observations (Zhou et al. 2025; Close et al. 2025a) can be affected by variable extinction (Cugno et al. 2025).

Based on nearly-coeval multi-frequency ALMA measurements, Domínguez-Jamett et al. (2025) suggest that, contrary to the standard hypothesis, the sub-millimetre signal is probably not due to thermal dust emission, as its optically thick spectral index over 97 GHz, 145 GHz, and 343 GHz (Band 7, hereafter IB23) is difficult to reconcile with the partially thick emission from a circumplanetary disk or envelope. It may instead be due to thermal plasma radiation, either with residual ionisation, as in a cool star (i.e. due to the H⁻ or H${}_2^{-}$ free-free continua), or fully ionised, as in an H II region. A 671 GHz drop suggests that the most likely explanation for the signal from PDS 70c is H I free-free from a thin H II region at the shocked surface of a circumplanetary disk (CPD) undergoing mass loading (as described by Aoyama et al. 2018; Szulágyi & Ercolano 2020).

The spectral analysis of Domínguez-Jamett et al. (2025) rests on the hypothesis that variability is negligible, which is justified by the short time span encompassing the multi-frequency observations, of two months at most. Indeed, the IB23 dataset, acquired in nine execution blocks (EBs) spanning six months, suggests that the Band 7 flux density of PDS 70c (hereafter F_B7) is constant within a standard deviation of 10%. In fact, F_B7 appears constant (Table 1) when comparing the discovery observation in the 2017 dataset (Isella et al. 2019), hereafter IB17 (following the nomenclature of Benisty et al. 2021), with measurements in 2019 (hereafter LB19, Benisty et al. 2021), in 2021 (Fasano et al. 2025), and in 2023 (i.e. in IB23, Domínguez-Jamett et al. 2025).

However, motivated by questions on the uncertainties in restored images and in the concatenation of observations, Casassus & Cárcamo (2022) and Casassus et al. (2023) chose PDS 70 as an example source to test an original scheme for multi-epoch visibility alignment (the VisAlign package) and to apply automatic self-calibration (with the snow package). Surprisingly, their revisit of IB17 yielded the non-detection of PDS 70c, despite sufficient signal-to-noise ratio (S/N) to reproduce the result of Isella et al. (2019). This non-detection led to the conclusion for variability of at least ∼40% over ∼2 yr.

We revised the analysis of Casassus & Cárcamo (2022) with the goal of resolving the discrepancy with all other measurements of F_B7. It turns out that IB17 is composed of three EBs, dated 2, 3, and 6 December. However, perhaps as a result of a bug in the handling of CASA tasks (which we attribute to intensive calls to the mstransform task, combining both time and channel averaging), the EB from 6 December was missing. As shown below, the inclusion of 6 December brings F_B7 into agreement with the original report by Isella et al. (2019). The variability reported by Casassus & Cárcamo (2022) therefore applies to timescales of days or less, rather than years.

Here, we report on a technique to estimate the flux density of unresolved yet variable embedded sources, and apply it to ALMA Band 7 datasets of PDS 70 (Sect. 2, the datasets are summarised in Table 1). We then discuss the expected timescales of the radio variability of PDS 70c for relevant emission mechanisms (Sect. 3). We conclude for a tentative detection of hourly variability, which favours an accretion shock at the CPD surface as the origin of the signal (Sect. 4).

2. Analysis and results

2.1. Individual EB imaging of IB17

Details on the alignment and self-calibration procedure are given in Appendix A. Imaging of the whole IB17 dataset (Fig. 1d, restored with Briggs parameter r = 0.3) confirms the detection of PDS 70c by Isella et al. (2019). Its flux density, given by the peak intensity at the position of the point source in restored maps, is ⟨F_B7⟩ = 172 ± 30 μJy, where the average represents a joint measurement over all three EBs. As detailed in Appendix A alternative measurements with Gaussian fits yield similar values within the errors (Table A.1). We note that ⟨F_B7⟩ is consistent at 1.8σ with the value reported by Isella et al. (2019).

Fig. 1. Restored images for IB17. The insets zoom in on PDS 70c (to the north-east), and symmetric positions along disk axes. The image noise and beams are provided in Table A.1. Peak intensities are all equal when brought to a common beam.

However, imaging of individual EBs, shown in Fig. 1, revealed no signal under PDS 70c on 2 and 3 December. Since the source is known to exist, rather than assigning zero flux density, we recorded the specific intensity at the expected position (Table A.1). This strategy for point-source photometry incorporates prior knowledge of the source position (as discussed by Domínguez-Jamett et al. 2025, their Sect. 2.3), and is sometimes referred to as ‘forced photometry’ (Samland et al. 2017).

The source is, however, quite conspicuous on 6 December, at 370 ±  65 μJy or ∼5.7σ (as given by the peak specific intensity). The F_B7 measurement on 6 December is roughly twice the average value of ⟨F_B7⟩, and the difference between 3 and 6 December is 276 ± 82.6 μJy and significant at 3.3σ. Fixing the average value of F_B7 as given by Domínguez-Jamett et al. (2025, Table 1), this rise corresponds to 228%±69%.

2.2. Time-differential photometry

The preceding estimates for F_B7 are, for simplicity, directly read out as the peak intensity under PDS 70c in restored images, where the source is separated from surrounding signal. But this approach requires artificially more weight in the longer baselines, overriding error propagation at the expense of sensitivity. In other words, with natural weights, the uncertainties would be lower, but the source would no longer be isolated from the outer ring (Fig. 1), thus introducing systematic errors.

We also tested the time variable nature of the radio signal from PDS 70c by measuring the relative variations in F_B7 in the uv plane. Our method was to fit a point source fixed at the location of PDS 70c (x_°, which was obtained with an elliptical Gaussian fit in 6 December), in the visibility residuals V^R after subtraction of the model visibilities calculated from the model image that best fits the concatenated dataset. Least-squares minimisation yields (notation follows from Casassus & Cárcamo 2022)

$$\begin{array}{c}\hfill F_{\nu }=\frac{{\sum }_{l=1}^L{\omega }_l\mathfrak{R}\left[V_l^R{e}^{-2\pi i{\mathit{u}}_l·{\mathit{x}}_{\mathbf{°}}}\right]}{\mathcal{A}({\mathit{x}}_{\mathbf{°}}){\displaystyle {\sum }_{l=1}^L}{\omega }_l},\end{array}$$(1)

where u_l is the uv-plane position for datum V_l^R and 𝒜(x_°) is the primary-beam attenuation. The associated errors are

$$\begin{array}{c}\hfill \sigma (F_{\nu })=\frac{1}{\mathcal{A}({\mathit{x}}_{\mathbf{°}})\sqrt{{\displaystyle \sum _{l=1}^L}{\omega }_l}}·\end{array}$$(2)

A similar measurement of F_B7 can also be obtained from the gridded visibilities. As shown in Appendix C, the gridded approach bypasses the model image, but it collapses the spectral information and depends on the gridding scheme. Since the model image is used in the self-calibration procedure anyways, we opted for photometric extraction on the raw data. Both approaches yield similar results for a flat spectral index (Fig. B.1).

Such time-differential photometry is deceivingly simple as it assumes perfect calibration. Any imperfection yields non-thermal residuals. For instance, omitting EB alignment with VisAlign, or omitting the amplitude self-calibration, yields such residuals, as documented in Appendix B and Fig. B.1. However, alignment and self-calibration allow reaching close to thermal residuals, even in natural weights.

Another difficulty in differential photometry is that the uncertainties σ(F_ν) derive from the visibility weights (Eq. (2)), whose absolute calibration is uncertain. This can be overcome with the replacement of the original weights by the root-mean-square dispersion of each visibility datum taken over all individual integrations, as is performed with the CASA task statwt. An additional imaging run is in order after the application of statwt, to check that it did not appreciably deteriorate the dynamic range.

A complementary method to estimate the uncertainty on uv-plane point-source photometry is to take statistics on flux densities extracted at randomly distributed positions. We refer to this estimate as the ‘spatial scatter’. At a given epoch, we applied Eq. (1) at N = 120 random positions, and took the dispersion of the resulting N flux densities as a measure of its uncertainty. In the presence of an extended and strong signal, this uncertainty depends on position since imperfect antenna gains are multiplicative. Given that the source under consideration is roughly a ring, we chose to place the random positions along a projected circular orbit matching that of PDS 70c (but excluding the position of PDS 70c within one natural-weight beam).

The close to thermal residuals justify the use of Eq. (1) to extract point-source flux densities in the visibility residuals, and measure ΔF_B7(t) = F_B7(t)−⟨F_B7⟩. The result is summarised in Fig. 2, where F_B7 is consistent between 2 and 3 December, but rises by 206 ± 56 μJy from 3 to 6 December, which is significant at 3.7σ. Fixing F_B7 to the value for IB23 in Table 1, this rise corresponds to 170%±46%. The quoted errors were derived from the spatial scatter, and they are very close to the errors given by the visibility weights (after the application of statwt).

Fig. 2. Time-differential photometry of PDS 70c. The blue ±1σ error bars (2σ in total) record ΔFB7 = FB7 − ⟨FB7⟩ as a function of Julian day for each EB. The same quantity extracted from the spatial scatter (see text) is shown with thick orange error bars. Each type of marker corresponds to three synthetic and constant point sources injected at the positions of the insets in Fig. 1.

As a test on the significance of the variability of PDS 70c, we injected three synthetic and constant point sources in the concatenated dataset at symmetric positions from PDS 70c relative to the disk minor and major axis. These positions match the centres of the insets in Fig. 1. Using the uncertainties inferred from the spatial scatter, as for PDS 70c, the largest variation in ΔF_B7(t) for the synthetic point sources is 122 ± 56 μJy or 2.2σ.

An application of differential photometry to other Band 7 datasets is documented in Appendix D. Neither LB19 nor IB23 present detectable variability among EBs.

It is intriguing that the daily variability of F_B7 is picked up in IB17, but not in IB23. The deep IB23 dataset spanned six months in nine EBs with integrations of over 2 h, while the IB17 EBs are 1 h long. Yet a search for variability, using a related approach to the time-differential photometry presented here (but subtracting a Gaussian fit to PDS 70c in the model image, so the source stands out in the residuals), yielded a constant value for F_B7, with a scatter of ∼10% among EBs (Domínguez-Jamett et al. 2025). Perhaps the source underwent a peculiar glitch right on 6 December. However, given that the IB23 EBs are more than twice as long in IB17, an interesting possibility is that the variability of PDS 70c occurs on timescales of minutes or hours rather than days, and that it is averaged out in longer integrations.

The difficulty in observing fine-grained time variability is that splitting into shorter integrations results in noisier measurements. As documented in Appendix D, we nonetheless picked up some scatter in ΔF_B7(t) in ∼20 min intervals of the IB23 dataset, corresponding to a reduced χ² of 1.7 with 41 degrees of freedom, or significant at 2.6 σ. The intrinsic rms variability of PDS 70c is σ_int = 59 ± 25 μJy.

3. Discussion

The constraints on the radio variability of PDS 70c bring information on the emission mechanisms in the circumplanetary environment. For example, radiative heating from accretion bursts could provide external forcing to the emitting medium. The radio variability would result from the intrinsic variability of accretion smoothed by the timescale associated with thermal diffusion.

The interpretation that F_B7 stems from an accretion shock at a CPD surface around PDS 70c (Domínguez-Jamett et al. 2025) naturally accounts for fast variability. Since the accretion flows within the Hill sphere can be approximated as free fall (e.g. Tanigawa et al. 2012; Fung & Chiang 2016), the free-free emission responds directly to fluctuation in the inflow from the protoplanetary disk, that can vary by factors of a few along the planetary orbit in turbulent disks (Gressel et al. 2013) or for eccentric planets (Li et al. 2023). However, in their hydro-dynamical models of the CPD surface shock, Takasao et al. (2021) find that the interaction between the infalling gas and the CPD triggers high-frequency and stochastic variability through thermal instability of the standing shock. The net mass accretion rate on the CPD surface, and eventually on the protoplanet itself, result from the interaction of post-shock cooling and flow dynamics in the CPD. This accretion rate fluctuates by up to an order of magnitude on timescales much shorter than the orbital period, and down to hours. By contrast, accretion occurs on a quasi-steady state if averaged over days, as in the proposal of Sect. 2.2 to explain the lack of detectable variability on longer timescales.

In the case of a surface shock, the H I free-free emission reflects the temperature profile across the shock (Aoyama et al. 2018). From the shock profiles shown in Domínguez-Jamett et al. (2025) for two representative pre-shock conditions, in the region where τ(300 GHz)∼1 and T ∼ 3 × 10⁴ K to 3 × 10⁵ K, at depths of ∼10⁴ cm and ∼10⁵ cm, respectively. Given that the main source of opacity is optically thin, except at radio frequencies, energy transfer occurs over a light-crossing time, which is essentially instantaneous.

Alternatively, the radio signal may stem from the bulk mass in the planetary environment under external forcing. Assuming energy transfer in an optically thick medium via Rosseland diffusion, the typical photon mean-free-path is l_R = 1/(κ_R ρ), where κ_R is the Rosseland opacity and ρ is the mass density. The diffusion timescale is τ_d = N_scat l_R/c, with N_scat = (L/l_R)², and where L is the typical size of the emitting region.

We distributed the mass in the planetary environment M_env, in a cubic box with side L, representing the linear size of a CPD or of an envelope. The outer radius of the CPD could be a third of the Hill radius, or ∼1 au for PDS 70c, and a typical aspect ratio is ∼0.1, so L ∼ 0.1 au. Then the diffusion timescale is τ_d = M_env κ_R/(c L). Reducing the box side, or increasing κ_R, results in larger diffusion timescales for a given mass.

The floor κ_R in a dust-deprived medium, with residual ionisation, is due to the H⁻ or H₂⁻ free-free continua and to atomic or molecular lines (e.g. Malygin et al. 2014). In a CPD accounting for the PDS 70c radio spectrum (Domínguez-Jamett et al. 2025), the annulus with the peak contribution to the net flux density, at R = 0.05 au, has a mid-plane ρ  ∼  10⁻⁸ g cm⁻³ and T  ∼  2000 K, and κ_R ≳ 0.1 (for solar metallicity, Malygin et al. 2014, their Fig. 3). The contribution from dust can hardly exceed this floor value, at least for large dust grains which have κ_R  ∼  0.1 cm² g⁻¹ for a standard gas-to-dust mass ratio of 100 (e.g. Chrenko et al. 2025).

An environment mass of M_env ∼ 0.07 M_Jup for a planet-to-CPD mass ratio of 10⁻² (considering a planet mass of 7 M_Jup, Wang et al. 2021) is, in practice, an upper limit. The corresponding dust mass of 0.2 M_⊕ is larger than the total mass of the Galilean moons, which is 0.07 M_⊕. It is also much larger than the dust mass inferred from modelling F_B7 assuming thermal dust emission, or ∼1.4 × 10⁻² M_⊕ (Shibaike & Mordasini 2024).

If M_env = 0.07 M_Jup, if L = 0.1 au, and if κ_R = 0.1 cm² g⁻¹, then τ_d = 3.5 d. Fixing τ_d = 1 h would require a box with a side of L = κ_R × M_env/(c × τ_d)∼8 au, which is much larger than the constraint of L < 1 au (from the limit on the angular extent of the source). Alternatively, the environment mass would have to be reduced to M_env = 8 × 10⁻⁴ M_Jup to match an hourly variability.

4. Conclusions

Imaging individual (∼1 h) EBs in the 2017 Band 7 observations of PDS 70 shows that PDS 70c is conspicuous on 6 December 2017, but undetectable on both 2 and 3 December. Estimates of the flux density from PDS 70c in the restored images (F_B7) result in a rise from 3 to 6 December, which is significant at 3.3σ.

Time-differential photometry allows the extraction of the relative variations of F_B7 in the uv plane, at a higher significance than in the restored images, where the source is confused with the outer ring. The residual Band 7 flux density of PDS 70c (ΔF_B7, after subtraction of the time-averaged extended emission) also rose from 3 to 6 December, by 3.7σ. Since the residuals in single EBs are not perfectly thermal, this estimate is tentative.

Tests on the variability of PDS 70c in Band 7 observations from 2023, with EBs twice as long than in 2017, yield a limit of ∼15%, as given by the dispersion of ΔF_B7 among seven EBs, which is consistent with a constant. However, splitting these 2 h EBs in 20 min intervals results in some intrinsic variability, significant at 2.6σ, and with an intrinsic rms scatter of 59 ± 25 μJy. These results tentatively suggest that the radio signal from PDS 70c is variable by ∼50%±21% on timescales of hours, but that its variability is averaged out on longer timescales.

As shown recently, existing models of surface accretion onto a CPD account for the PDS 70c radio spectrum (Domínguez-Jamett et al. 2025). Such models also predict high-frequency and stochastic planetary accretion variability, on timescales of hours, and averaged to a steady state on longer timescales (Takasao et al. 2021). If the radio signal stems from the surface shock itself, then the observed radio variability directly follows the accretion variability. In turn, if the radio signal stems from bulk mass in the planetary environment, then the environment mass must be at least a factor ∼10⁻⁴ lower than the planetary mass to bring the diffusion timescale down to 1 h.

Acknowledgments

We thank the referee, Luis Peña-Moñino, for a constructive report, and also Myriam Benisty and Daniele Fasano for sharing their version of the IB17 dataset, which included Dec. 6. S.C. and M.C. acknowledge support from Agencia Nacional de Investigación y Desarrollo de Chile (ANID) given by FONDECYT Regular grant 1251456, and ANID project Data Observatory Foundation DO210001. G.-D.M. acknowledges the support of the Deutsche Forschungsgemeinschaft (DFG) through grant MA 9185/2-1. O.C. acknowledges the Czech Science Foundation (grant 25-16507S), the Charles University Research Centre program (No. UNCE/24/SCI/005), and the Ministry of Education, Youth and Sports through the e-INFRA CZ (ID:90254). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.A.00006.S, ADS/JAO.ALMA#2018.A.00030.S, ADS/JAO.ALMA#2022.1.01477.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

References

Appendix A: Alignment, self-calibration and image-plane F_B7 measurements in the IB17 dataset

Description of the IB17 ALMA dataset, alignment of EBs and imaging follow from Casassus & Cárcamo (2022), Casassus et al. (2023), and references therein. In brief, the 3 continuum spectral windows were time-averaged to 6.06 s, and binned into 8 channels. Individual EBs from 3 and 6 December 2017 were then aligned to the first EB, from 2 December, using the VisAlign package, including a correction for visibility amplitudes α_R. The resulting offsets were α_R = 1.102 ± 0.002, $\mathrm{\Delta }\alpha =0\stackrel{″}{.}00845\pm 0\stackrel{″}{.}00034$, $\mathrm{\Delta }\delta =-0.{00045}_{-0.00044}^{+0.0004}$, for 3 December, and α_R = 1.225 ± 0.003, $\mathrm{\Delta }\alpha =0\stackrel{″}{.}0107\pm 0\stackrel{″}{.}0004$, $\mathrm{\Delta }\delta =0\stackrel{″}{.}0043\pm 0.0005$, for 6 December. An offset in flux scale of ∼20% is consistent with the nominal calibration uncertainty in ALMA Band 7, of ∼10% rms, while the positional offsets along R.A. and Dec. are all consistent with the nominal pointing accuracy, of ∼1/10 the naturally weighted beam of $0\stackrel{″}{.}09\times 0\stackrel{″}{.}08$.

Following alignment, we concatenated the data into a single measurement set, and applied the snow package for automatic self calibration (snow replaces the tclean imager by gpu-uvmem, Cárcamo et al. 2018, in otherwise standard CASA self-calibration). Phase-only self-calibration did not yield improvements. However, a single iteration in amplitude and phase self-calibration, with a solution interval limited to a single scan (roughly 50 s), improved the peak S/N from 55 to 115.

Once self-calibrated, each EB was split for separate imaging with gpu-uvmem, in pure-χ² reconstructions (i.e. without entropy regularisation, but truncating the optimisation at ten iterations). The (non-parametric) model images were restored using Briggs weights with a robustness parameter r = 0.3. The resulting images are shown in Fig. 1.

Table A.1 provides a summary of the image properties in the IB17 dataset. We also tabulate the values of F_B7 measured in each EB as well as in the concatenated dataset. The elliptical Gaussian fits were carried out following two strategies: 1- by fitting a beam, assuming that the source is unresolved, and 2- by fitting an elliptical Gaussian. Fitting a beam resulted in a bias for larger flux densities to account for extended low-level emission. Fitting an elliptical Gaussian was sensitive on the definition of the zoomed region around PDS 70c, and turned out to be impossible in the case of the concatenated dataset, where the source is connected to the ring.

Appendix B: Residual maps on single EBs and impact of alignment and self-calibration

Reaching thermal residuals in a dataset comprised of several EBs is routinely achieved, provided with sufficiently accurate calibration. However, because of the varying flux scale and pointing, the residuals under individual EBs will, in general, not be thermal. The application of both EB alignment (with VisAlign), along with amplitude and phase self-calibration, allowed reaching close to thermal residuals in the 2017 dataset, as illustrated in Fig. B.1. The omission of alignment or self-calibration leads to strong residuals, due to imperfect calibration of the visibility data.

Fig. B.1. Dirty map residuals in natural weights for each EB, after subtraction of the model visibilities from the concatenated dataset. Insets zoom on PDS 70c, in the same colour stretch as for the image. Each row, from top to bottom, shows: ‘no visalign, apcal’, without visalign alignment, but after amplitude and phase self-calibration (apcal); ‘visalign, no apcal’: with alignment, but no self-calibration; ‘visaligned, ap’: here with both alignment and self-calibration, and reaching close to thermal residuals; ‘gridded’: also with both alignment and self-calibration, but from the gridded approach (i.e. Idiff in Eq. C.3); ‘restored’: restored maps, obtained with r = 0.3 (as in Fig. 1).

Amplitude self-calibration of the IB17 dataset required short solution intervals, down to a single scan or 50 s. As a result, the derived antenna gains can be fairly noisy. We found the resulting gains vary strongly, up to ∼100% rms. These strongly variable gains modulate the point-source flux densities, in ways that cannot be corrected for with a single scaling factor (as with VisAlign). These noisy gains are probably at the origin of the persistent non-thermal residuals. Improvements in self-calibrations techniques, such as with ‘closure imaging’ (Chael et al. 2018), might allow purely thermal residuals in this dataset.

Appendix C: Gridded uv-plane time-differential photometry

The photometric extraction of the variable signal can also be measured in gridded visibility data. We can compute the gridded visibilities, in natural weights, for both the concatenated dataset, ${\stackrel{\sim }{V}}_{\mathrm{conc}}$, with weights W_conc, and in each individual EB, ${\stackrel{\sim }{V}}_{\mathrm{EB}}$, with weights W_EB. In an extension of the scheme used in the visalign package (Casassus & Cárcamo 2022), the residuals are

$$\begin{array}{c}\hfill {\stackrel{\sim }{V}}_{\mathrm{diff}}={\stackrel{\sim }{V}}_{\mathrm{EB}}-{\stackrel{\sim }{V}}_{\mathrm{conc}}\end{array}$$(C.1)

with weights

$$\begin{array}{c}\hfill W_{\mathrm{diff}}=\frac{W_{\mathrm{conc}}\times W_{\mathrm{EB}}}{W_{\mathrm{conc}}+W_{\mathrm{EB}}}.\end{array}$$(C.2)

Such weights naturally match the uv sampling functions of both datasets in comparison.

The dirty images of the variable signal in each EB can be used as a consistency test between the photometry in the raw and gridded approaches. The point-source flux density in Eq. 1 is a function of position, and can also be viewed as the dirty map of the residual visibilities corrected for primary-beam attenuation. In the gridded approach, the dirty map is

$$\begin{array}{c}\hfill I_{\mathrm{diff}}(\mathit{x})=\frac{{\sum }_{i,j}^N{W}_{\mathrm{diff}}({\mathit{u}}_{i,j}){\stackrel{\sim }{V}}_{\mathrm{diff}}({\mathit{u}}_{i,j})e^{-2\pi i{\mathit{u}}_{i,j}·\mathit{x}}}{{\sum }_{i,j}^N{W}_{\mathrm{diff}}({\mathit{u}}_{i,j})},\end{array}$$(C.3)

in units of Jy beam⁻¹, and for a grid of side N.

Fig. B.1 compares the dirty maps of the best visibility residuals, i.e. after both alignment and self-calibration, for both approaches. The row labelled ‘visalign, apcal’ refers to the dirty maps of the raw residuals, after subtraction of the model visibilities that fit the concatenated dataset, while the row labelled ‘gridded’ refers to the direct subtraction of the concatenated dataset, i.e. ${\stackrel{\sim }{V}}_{\mathrm{difff}}$. Both version yield similar images, with differences below the noise level, which satisfies this consistency test.

Appendix D: Application to deep Band 7 data

We applied differential photometry to the LB19 Band 7 dataset (Benisty et al. 2021), which was acquired in 3 EBs, two from 27/7/2019, and one from 30/7/2019, and each ∼ 1h15m long. Differential photometry yields the residual flux densities ΔF_B7 shown in Fig. D.1, where no variability can be picked up among the three EBs.

Fig. D.1. Same as Fig. 2, but for the LB19 dataset. The blue error bars stem from the original visibility weights (in this case without the application of statwt), while the orange bars stem from the spatial scatter.

We also revisited the full reduction and alignment of the IB23 dataset, and applied time-differential photometry for each EB, as shown in Fig. D.2. The weighted-average residual flux density is ΔF_B7 = 0.2 ± 8 μJy, with a weighted scatter of 18 μJy, which represents ∼15% of the flux density reported by (Domínguez-Jamett et al. 2025, of 121 ± 13 μJy). The reduced χ² is 0.6 (with 6 degrees of freedom), and consistent with a constant distribution. We therefore reach the same conclusion as Domínguez-Jamett et al. (2025), in the sense that the observed scatter stems purely from the thermal uncertainties (which appear to be somewhat conservative).

Fig. D.2. Same as Fig. 2, but for the IB23 dataset (2023 Band 7 observations).

In an attempt to pick-up fine-grained variability, we split the seven EBs from IB23 in six equal intervals, each ∼20 min long, and aligned the resulting ‘sub-EBs’ using VisAlign. The reference EB was chosen as 2/6/2023, as it was acquired under the best conditions. The corrections in flux scales were all less than 2%, and the pointing corrections less than 16 mas, which reflects the excellent weather. We then concatenated the aligned sub-EBs for joint imaging, and applied differential photometry (i.e. we used the residuals for the photometry of PDS 70c in the uv-plane).

The fine-grained measurements of ΔF_B7 shown in Fig. D.3 seem to be indicative of variability, as the reduced χ² is 1.67 (with 41 degrees of freedom). This corresponds to a detection of variability at 99.6% confidence, or 2.62 σ. The weighted scatter is σ_s = 13 ± 2 μJy, with bootstrapped uncertainties. The intrinsic variability of PDS 70c, σ_int, can be estimated as an additional source of noise, added in quadrature to each measurement of ΔF_B7(t), such that the reduced χ² of the PDS 70c time-series is brought down to unity. In this case we find σ_int = 59 ± 25 μJy, with bootstrapped uncertainties. Fixing F_B7 to the value for IB23 in Table 1, this scatter corresponds to 49%±21%.

Fig. D.3. Same as Fig. 2, but for the fine-grained IB23 dataset, splitting each EB into 6 intervals. The blue error bars stem from the statwt visibility weights, while the orange bars stem from the spatial scatter. The coarser measurements shown in transparent red boxes, with total vertical extension of 2σ, are reproduced from Fig. D.2, and correspond to extractions over the full EB.

However, the limitations on self-calibration discussed in Sect. B suggest that the rapidly varying antenna gains represent a source of systematic noise that cannot be corrected through alignment to a global flux scale with VisAlign. Also, such non-thermal residuals would not be accounted for in the spatial scatter along the orbit of PDS 70c, since that region is approximately devoid of signal. A test on synthetic and constant point sources is required, with the same flux density as PDS 70c, and at the positions described in Sect. 2. These three control point-sources have ΔF_B7(t) time series which result in reduced χ² values of 1.35, 1.29 and 1.36, which are all consistent with a constant value (the corresponding significance of a non-constant value is ∼1.4σ).

All Tables

All Figures

	Fig. 1. Restored images for IB17. The insets zoom in on PDS 70c (to the north-east), and symmetric positions along disk axes. The image noise and beams are provided in Table A.1. Peak intensities are all equal when brought to a common beam.
In the text

	Fig. 2. Time-differential photometry of PDS 70c. The blue ±1σ error bars (2σ in total) record ΔFB7 = FB7 − ⟨FB7⟩ as a function of Julian day for each EB. The same quantity extracted from the spatial scatter (see text) is shown with thick orange error bars. Each type of marker corresponds to three synthetic and constant point sources injected at the positions of the insets in Fig. 1.
In the text

Fig. B.1. Dirty map residuals in natural weights for each EB, after subtraction of the model visibilities from the concatenated dataset. Insets zoom on PDS 70c, in the same colour stretch as for the image. Each row, from top to bottom, shows: ‘no visalign, apcal’, without visalign alignment, but after amplitude and phase self-calibration (apcal); ‘visalign, no apcal’: with alignment, but no self-calibration; ‘visaligned, ap’: here with both alignment and self-calibration, and reaching close to thermal residuals; ‘gridded’: also with both alignment and self-calibration, but from the gridded approach (i.e. Idiff in Eq. C.3); ‘restored’: restored maps, obtained with r = 0.3 (as in Fig. 1).

In the text

	Fig. D.1. Same as Fig. 2, but for the LB19 dataset. The blue error bars stem from the original visibility weights (in this case without the application of statwt), while the orange bars stem from the spatial scatter.
In the text

	Fig. D.2. Same as Fig. 2, but for the IB23 dataset (2023 Band 7 observations).
In the text

	Fig. D.3. Same as Fig. 2, but for the fine-grained IB23 dataset, splitting each EB into 6 intervals. The blue error bars stem from the statwt visibility weights, while the orange bars stem from the spatial scatter. The coarser measurements shown in transparent red boxes, with total vertical extension of 2σ, are reproduced from Fig. D.2, and correspond to extractions over the full EB.
In the text