© The Authors 2025

1 Introduction

PDS 70 hosts two protoplanets detected in direct imaging, and probably even three (Keppler et al. 2018; Haffert et al. 2019; Wang et al. 2021; Christiaens et al. 2024; Hammond et al. 2025). Atacama Large Millimeter Array (ALMA) observations of PDS 70 resulted in the detection of a possible circumplanetary disc (CPD) around PDS 70c (Isella et al. 2019; Benisty et al. 2021), in the form of a point source in continuum emission at millimetre wavelengths in Band 7 (351 GHz), with a linear diameter of less than 1.2 au, or about one third the size of the Hill radius of the accreting protoplanet. This CPD is likely to play a crucial role in regulating the accretion process of PDS 70c. Its sub-millimetre flux density is variable by ∼40% on timescales shorter than 2 yr, and could be explained by the free-free continuum concomitant to the unresolved Hα signal (Casassus & Cárcamo 2022) reported by Haffert et al. (2019). For every proton-electron encounter leading to the recombination cascade, others will lead to free-free transitions.

The global hydrodynamic simulations of embedded thermalmass bodies by Szulágyi et al. (2018), as well as the parametric viscous disc models by Zhu et al. (2018), suggest that higher frequencies, such as those in ALMA Bands 9 and 10 (∼671 GHz and ∼850 GHz), are more suitable for CPD observations in the dust continuum. Circumplanetary discs are predicted to be hotter than the surrounding circumstellar disc, and their (unresolved) emergent flux densities, F_ν, in frequency, are thought to be mainly due to partially optically thick thermal dust with a spectral index of α ≈ 3 at ALMA frequencies, with α defined by F_ν ∝ ν^α.

However, if the radio signal from PDS 70c is mainly due to free-free emission, we should expect α ≤ 2. A rough approximation of partially thick free-free emission can be obtained with a 1/r² electron-density distribution, as in winds from early-type stars. The emergent free-free flux densities then follow a power-law spectrum, F_ν ∝ ν^α, with α ∼ 0.6 (Wright & Barlow 1975).

The sub-millimetre detection of PDS 70c thus raises the questions of what the nature of the emission is, and why PDS 70c is detected in the sub-millimetre continuum, while PDS 70b is not. If the emission is due to thermal dust emission, a multiwavelength analysis could shed light on the nature of the circumplanetary dust population. In turn, the free-free continuum component will inform us directly about the planetary accretion process, without the large uncertainties related to the levels of extinction that affect Hα (e.g., Hashimoto et al. 2020).

In this study, we present new ALMA Band 4 (145 GHz), Band 7 (343.5 GHz), and Band 9 (671 GHz) observations of PDS 70, and compare them to archival data in Band 3 data (97.5 GHz) from Doi et al. (2024), focussing on PDS 70c. Section 2 describes data acquisition, synthesis imaging, and the extraction of the spectral energy distribution (SED) of PDS 70c. Sect. 3 analyses our results in terms of parametric disc models. Sect. 4 discusses which models can or cannot fit the data. Sect. 5 concludes.

2 Observations

2.1 Data acquisition

New ALMA observations of PDS 70 in Bands 4, 7, and 9 were obtained as part of ALMA project IDs 2022.1.00592.S and 2022.1.01477.S. Each band was acquired in an extended array configuration, labelled TM1, and a compact configuration, labelled TM2. An observation log is provided in Table 1. We also include in our analysis the archival datasets for Bands 3, from ALMA projects 2022.1.00893.S (Doi et al. 2024) and 2022.1.01477.S (Liu et al. 2024).

2.2 Imaging, self-calibration and alignment

The same imaging strategy as for the 2017 and 2019 Band 7 datasets (Casassus & Cárcamo 2022; Casassus et al. 2023) was extended to the other bands. For each band and array configuration, automatic (user-independent) self-calibration was performed using the snow package. This package replaces the CASA tclean imager (CASA Team 2022) with gpuvmem (Cárcamo et al. 2018) in CASA self-calibration (i.e. in iterative calls to the gaincal and applycal CASA tasks). Image restoration was performed with a Briggs robustness parameter of r = 2. The gaincal solution interval was reduced by half in each iteration, starting with the scan length, and stopping when the peak signal-to-noise ratio (PS/N) decreases. The solution table with the highest PS/N was kept. After the phase calibration stage, a cycle of phase and amplitude self-calibration was performed to assess any PS/N improvements. During the antenna-gain calculations, using CASA task gaincal, all spectral windows were combined. Alignment was performed using the VisAlign package (Casassus & Cárcamo 2022; Casassus et al. 2023), selecting the datasets with the largest baselines as an astrometric reference (i.e. TM1), and a final round of self-calibration with snow was applied to the concatenated dataset. The resulting images are summarised in Fig. 1.

For Band 3, we concatenated the two extended (TM1) configurations from projects 2022.1.00893.S, 2022.1.01477.S. The PS/Ns of these data are limited by thermal noise, and selfcalibration did not produce any improvements. The datasets were acquired with different phase centres, and were aligned with VisAlign before concatenation. The reference was taken as 2022.1.00893.S. The resulting Band 3 continuum dataset yielded a PS/N of 13.6 with r = 2, limited by thermal noise. As is shown in Fig. 1a, a bright point source is detected at the phase centre, which is coincident with the position of PDS 70 given by the GAIA ICRS co-ordinates (Gaia Collaboration 2023). We therefore assume that this point source is due to accretion free-free emission on the central star, following Doi et al. (2024).

For Band 4, each array configuration (TM1 and TM2) was self-calibrated separately prior to concatenation. The PS/N for TM2 improved from 224 to 243 after one round of phase self-calibration and one round of combined phase and amplitude self-calibration. In contrast, TM1 is limited by thermal noise, with a PS/N of 49. An application of VisAlign, with TM1 as a reference, showed that the pointing is very consistent in these datasets, with a shift smaller than 1 mas, and consistent intensities within 1%. Such remarkable consistency probably reflects the favourable observing conditions. Joint selfcalibration after concatenation did not yield improvements in the PS/N.

The new Band 7 data is coarser than in Benisty et al. (2021), with a restored beam in natural weights of 0″.102 × 0″.082 in TM1. We avoided concatenation with TM2 to preserve angular resolution. The PS/N of TM1 improves from 104 to 257 after two rounds of phase self-calibration, and reaches 308 after amplitude and phase self-calibration. However, we selected the phase-only self-calibration for our final dataset, as the amplitude gains result in a drop in flux density for PDS 70 by 20%, compared to both TM2 and the phase self-calibrated TM1.

In Band 9, the PS/N for TM2 increases from 29.2 to 33.3 after two rounds of phase calibration, and to 111 after one round of combined amplitude and phase self-calibration. TM1 is affected by a low-spatial frequency modulation, probably caused by poor weather and resulting in noisy phases. We opted to discard baselines shorter than 0.4 Mλ. The uvrange-clipped dataset yielded a PS/N of 24, which did not improve with selfcalibration. An application of VisAlign to align TM2 to TM1 results in a shift of ∆α = 25.6 mas along the right ascension, and ∆δ = 0.13 mas along the declination, with a correction in the intensity scale of 78%. While seemingly large, the positional shift is within the pointing uncertainty of the TM2 dataset, which is about one tenth of the beam or 15 mas (at 1 σ). Likewise, the absolute flux calibration accuracy in Band 9 is around 10%. Failing to apply the positional shift to TM2 results in images of the concatenated dataset with a spurious excursion of the ring into the cavity. We applied the correction factor of 1/78% to TM1. Self-calibration of the concatenated dataset did not yield improvements.

While the star PDS 70 is readily picked up in Band 3, it bears no clear counterparts in the other bands. There may be a point-source associated with the star in Band 4, where the central emission appears to be split into four blobs in the deconvolved image of Fig. 1f, one of which coincides with the expected stellar position with 3 σ flux. The other blobs are offset by ~50 mas, which is too large to account for by a pointing error in Band 4, given by ~1/10 of the Band 4 TM1 beam, or ~8 mas. We fitted elliptical Gaussians to extract the stellar flux densities. Comparing 24 ± 5 μJy in Band 3 with 18 ± 4 μJy in Band 4 corresponds to a spectral index of α = ~0.7 ± 0.8, which is too noisy to distinguish partially thick free-free emission (α ~ 0.6, Wright & Barlow 1975), variability, or a synchrotron component. We applied VisAlign to check on the inter-band alignment, using Band 3 as reference given the detection of the star. Band 4 is shifted by ∆α = ~4.9 ± 1.1 mas, ∆δ = ~13 ± 1.1 mas, which is within the pointing errors in Band 4.

There is no counterpart of the star in Bands 7 or 9 with which to cross-check the inter-band alignment, and the ring appears dramatically different in radius and width compared to Band 3, hampering the application of VisAlign. We simply rely on the absolute pointing accuracy, of ~10 mas in Band 7 and ~5 mas in Band 9.

Fig. 1 Multi-frequency imaging of PDS 70. The black insets at the top left of each image zoom on the central region (with tick-marks separated by 0".1), while the orange insets at the bottom right zoom on the expected position of PDS 70c for circular Keplerian rotation around a 0.97 M⊙ star (with tick-marks at 0″.05). The red plus symbol in the central inset marks the nominal position of the star. In all insets the linear grey scale stretches over the range of intensities in each region, and the contour levels start at 3 σ and are incremented in units of σ. Intensities within each contour level are colour-coded differently. The two circles are centred on the positions of PDS 70b and PDS 70d (no radio counterparts are detected). Images in a) to d), along the top row, correspond to restorations of the corresponding non-parametric model images (that fit the visibility data), in e) to h), along the bottom row. A beam ellipse in a) to d) is shown in blue, on the top right. The pixel size in the model images is fixed at 4 × 4 mas2. For each band the Briggs parameter used for the restoration, the clean bean and the noise can be found in Table 2.

2.3 Point source measurements of PDS 70c

We searched for radio counterparts to PDS 70 b, c, and d by estimating their positions at the epochs of observations, under the assumption that their orbit is circular around a 0.97 M_⊙ star, and in the plane of the disc, with position angle PA = 160 deg and inclination i = 130 deg (e.g. Casassus & Cάrcamo 2022). The two circles in each of the images in Fig. 1 are centred on the positions of b and d, using the co-ordinates provided by Wang et al. (2021, from UT2020-02-10) and Christiaens et al. (2024).

No counterparts were detected for PDS 70b or d, at any frequency. However, we did find counterparts to PDS 70c in Bands 3 and 4. In Band 3 a point-source signal at 2.6σ in the restored map is found at ∆α = 9 ± 8 mas, ∆δ = ~1 ± 8 mas, from the expected position of PDS 70c, as was estimated with an elliptical Gaussian fit to the Band 7 source (see Table 3). Given the pointing uncertainties in Band 7 (of ~5 mas), the uncertainties in the orbit, and the faint Band 3 counterpart that can be modulated by thermal noise, a shift of 9 mas is well within the errors. We recorded the peak intensity in Jy beam⁻¹ as a measurement of the point-source flux, and set its uncertainty to the noise in the residual image.

In Band 4, a point source at 5.2 σ is found at ∆α = 7 ± 7 mas, ∆δ = 2 ± 7 mas from the expected position of PDS 70c. The centroid of the point source was estimated with an elliptical Gaussian fit in the deconvolved image. Our choice of a Briggs parameter of r = 1 in Sect. 2.2 for image restoration is a compromise in terms of resolution to separate PDS 70c from the outer ring and PS/N.

PDS 70c is very conspicuous in the new Band 7 data, especially in the deconvolved image of Fig. 1g. The shift from its expected position is ∆α = 13 ± 12 mas, ∆δ = −6 ± 12 mas. However, its flux density is difficult to measure because of confusion with the cavity wall. The peak on PDS 70c, in the restored map of Fig. 1c, is 111 ± 14 μJy, and a Gaussian fit to the deconvolved map of Fig. 1g yields 157 ± 27 μJy. We report the weighted average, 121 ± 13 μJy.

Intriguingly, no counterparts of PDS 70c were found in Band 9. and the flux was obtained with the intensity of the peak in a region centred on its expected position, and with a radius given by 3 times the root-mean-square (rms) pointing accuracy (of ~5 mas). We did not attempt to correct for the diffuse emission near the edge of the cavity, which would lower the flux density reported in Table 3. As the absence of PDS 70c in Band 9 is particularly surprising, we performed point-source injections tests in the uv plane, using optically thick spectra tied to Band 7. As is shown in Fig. 2, such point sources should indeed be detectable, at least in the form of a protrusion of the cavity wall.

The Band 3 and Band 9 flux density measurements of PDS 70c presented here do not correspond to detections. However, their exclusion from the resulting SED would represent a severe loss of information. In turn, incorporating these nondetections as 2 or 3 σ upper limits injects a bias in attempts to fit this SED, by setting the expectation value for the flux densities to zero (or any other arbitrary value). The prior knowledge of the existence of PDS 70c, and of its orbit, is incorporated by reporting the peak flux density, within a beam, at the expected location of PDS 70c. This approach was referred to as ‘forced photometry’ by Samland et al. (2017), for the incorporation of non-detections in the SED of 51 Eri b.

With the flux measured in different bands, we use the spectral index, α, as a useful indicator of the emission mechanisms at radio frequencies, $$\alpha (\nu )=\frac{d\ln (F_{\nu })}{d\ln (\nu )}.$$(1)

We estimated α in PDS 70c using the measurements from Table 3, by comparing frequency pairs, $${\alpha }_{{\nu }_1}^{{\nu }_2}=\frac{\ln (F_{{\nu }_2}/F_{{\nu }_1})}{\ln ({\nu }_2/{\nu }_1)},$$(2)

and by fitting a power law to more than two frequency points, $$F_{\nu }=F_{{\nu }_{\circ }}{\left(\frac{\nu }{{\nu }_{\circ }}\right)}^{\alpha },$$(3)

with ν° fixed to 145 GHz.

A summary of different spectral index estimates is given in Table 4. Measurements involving Bands 3, 4, 7, and 9 in 2023 are spread over two months (except for the November 2023 execution blocks in Band 7 and the December 2023 TM2 execution block in Band 9), and are the least biased by variability compared to the 2017 and 2019 data.

Fig. 2 Point-source injections tests in Band 9. (a) Same as Fig. 1, with additional insets. (b) Same as a) but with four point sources, injected at the centre of each inset, and with a flux density corresponding to an optically thick extrapolation of the Band 7 flux density for PDS 70c.

Fig. 3 Dirty maps in natural weights of the residual Band 7 visibilities, after subtraction of the disc. The grey scale stretches over the full intensity scale, and for each of seven epochs. The inset zooms into PDS 70c, with contours at 90 and 120 μJy beam−1, fixed for ease of comparison across epochs. As in Fig. 1, the circles are centred on PDS 70b and d, and are drawn here to ease the search for any variable counterparts.

2.4 Variability study

The new Band 7 data were acquired in nine execution blocks spanning seven epochs, and each of these execution blocks reached thermal noise levels of ~20-30 μJy beam⁻¹, in natural weights, which should be deep enough to detect PDS 70c. However, the natural weight beams, of ~0″.12, are too coarse to separate PDS 70c from the bright ring. We therefore proceeded to subtract the disc emission from the visibility data. In this case we opted to work on the visibility data after amplitude and phase self-calibration, which ensures a consistent flux scale across all epochs. This was ascertained by aligning all epochs to 23 May 2023 with VisAlign, resulting in positional shifts of less than 1 mas, and amplitude shifts of less than 1%. We prepared a disc image with the gpuvmem model of the resulting visibilities (same as in Fig. 1g), and subtracted PDS 70c with an elliptical Gaussian fit. Simulated visibilities on this model image were then computed with the pyralysis package (Casassus & Cárcamo 2022) and subtracted from the observations. The dirty maps of the residuals are shown in Fig. 3, in natural weights (and also computed with pyralysis). The residual intensities are not perfectly thermal, as the model image fits the averaged dataset, not individual executions blocks. Nevertheless, PDS 70c stands out as the brightest residual in all epochs.

The technique of uv-plane subtraction of the disc signal allowed us to build visibility data containing mostly the signal from PDS 70c, with a consistent flux scale. We could thus proceed with point-source fits in the uv plane, fixed at x°, the position of the elliptical Gaussian centroid used to subtract PDS 70c, but with the best-fit amplitude (obtained by weighted least squares minimisation): $$F_{\nu }=\frac{\sum _k{\omega }_k\mathrm{\Re }\left[V_k^{\circ }\exp (-2\pi i\overrightarrow{u}\cdot \overrightarrow{x_{\circ }})\right]}{{\mathcal{A}}_{\circ }\sum _k{\omega }_k},$$(4)

where $\mathrm{\Re }[\cdot$ denotes the real part, and with rms uncertainties $$\sigma \left(F_{\nu }\right)=\frac{1}{{\mathcal{A}}_{\circ }\sqrt{\sum _k{\omega }_k}},$$(5)

where the sums extend over all observed visibilities, $V_k^{\circ }$, at the frequency, ν, with weights, ωk, and with a primary beam attenuation of ${\mathcal{A}}_{\circ }=\mathcal{A}({\overrightarrow{x}}_{\circ })$ at the position, x°. We stress that while this technique is suitable for a variability study, it cannot be used to extract the flux density of the point source in the concatenated dataset, since the model image fits this point source.

In practice, we collapsed all channels in each spectral window, yielding four flux-density measurements. The weighted average flux density at each epoch is recorded in Fig. 4. The weighted average is 130 ± 4 μJy, with a weighted scatter of 9.3%. Therefore the flux density of PDS 70c is constant across all seven epochs within ~10%. A non-detectable variability of the Band 7 flux densities the PDS 70c in 2023 is also consistent, given the uncertainties, with the conclusion from Fasano et al. (2025), who compared measurements in 2019, 2021, and 2023. The Band 4–Band 7 spectral index derived here is consistent with the Band 6–Band 7 index reported by Fasano et al. (2025).

Fig. 4 Flux density of PDS 70c at 343.5 GHz (Band 7), as a function of Julian day. The scatter is consistent with the uncertainties.

3 Analysis

From the spectral indices in Table 4, we note that optically thick emission, that is α ≳ 2, is consistent with the data between Bands 3 and 7 but only marginally consistent when we add Band 9. This is very difficult to reconcile with dust emission, where the emergent total (unresolved) flux densities mix optically thick and thin regions, resulting in spectral indices larger than 2.5 and closer to 3, depending on the choice of grain population (Zhu et al. 2018; Shibaike et al. 2025). In fact, a blackbody fit to Bands 3, 4, and 7 predicts 465 ± 42 µJy in Band 9, and overshoots the observation by 2.6 σ.

The non-detection of PDS 70c in Band 9 is intriguing, as CPD models predict a bright signal at such short wavelengths, and favour Band 9 for CPD searches (Szulágyi et al. 2018; Zhu et al. 2018). It seems unlikely that this non-detection is caused by variability, given Sect. 2.4, and since PDS 70c is readily detected in Bands 4 and 7 at ∼5 σ and ∼9 σ, respectively. Instead, the multi-frequency data suggest that free-free emission, rather than the thermal dust component, may be the source of the PDS 70c radio signal. In this section we analyse our results in terms of simple parametric models of the radio signal from CPDs, which we use to reproduce the measurements from Bands 3, 4, 7, and 9.

3.1 PDS 70c SED models

The lack of CPD detection in Band 9 may be indicative of free-free radiation. To test this hypothesis, we extended the parametric CPD model proposed by Zhu et al. (2018) to include the free-free emission expected from the CPD itself (here we consider the radio free-free contributions from all species). For completeness we reproduce here some of the formulae from Zhu et al. (2018), before incurring some approximations and including our extensions.

The model of Zhu et al. (2018) incorporates H I free-free emission in the form of a jet. For this component, we used the formula from Anglada et al. (2018) to estimate the emission of a circumstellar disc at 1 cm: $$\left(\frac{F_{1\mathrm{c}\mathrm{m},\mathrm{j}\mathrm{e}\mathrm{t}}{d}^2}{\mathrm{m}\mathrm{J}\mathrm{y}\mathrm{k}\mathrm{p}{\mathrm{c}}^2}\right)={10}^{-1.9}{\left(\frac{L_{\mathrm{a}\mathrm{c}\mathrm{c}}}{L_{\odot }}\right)}^{0.6},$$(6)

with $$L_{\mathrm{a}\mathrm{c}\mathrm{c}}=\frac{G{M}_{\mathrm{p}}{\dot{M}}_{\mathrm{p}}}{2R_{\mathrm{i}\mathrm{n}}},$$(7)

where M_p is the mass of the planet, Ṁ_p the accretion rate onto it, and R_in the inner edge of the CPD.

The 1 cm flux density of the jet may be extrapolated to millimetre wavelengths under the hypothesis of a power-law spectrum. If the jet is partially optically thick, an upper limit is set on the spectral index with α_jet = 0.4 (Carrasco-González et al. 2016). In the case of an optically thin jet, α_jet = −0.1. It should be borne in mind that these formulae were developed for stellar jets, rather than protoplanets, and that they are being extrapolated into the planetary mass regime.

The parametric CPD model from Zhu et al. (2018) calculates the radio emission based on the disc accretion rate, Ṁ_p, and viscosity parameter, α. Prior to applying this model to PDS 70c, we made sure to reproduce the tables from Zhu et al. (2018) for the boundary layer case (this assumes that the planet accretion luminosity dominates over the central brightness of the planet itself), and in the Rayleigh-Jeans approximation for consistency.

We opted for a planetary mass of M_p = 4 M_Jup, given the upper limit of 4.9 M_Jup (Doi et al. 2024, based on a circumstellar disc viscosity parameter of α = 10⁻³) and the restriction of 4 M_Jup < M_p < 12 M_Jup (Haffert et al. 2019). PDS 70 is at a distance of 112.3 pc (Gaia Collaboration 2023), and we adopted a planetary orbital radius for PDS 70c of 34 au (Haffert et al. 2019), implying an outer CPD radius of R_out = 1/3 R_Hill ~ 1.23 au (Quillen & Trilling 1998), corresponding to an angular size of Θ = 0″.011. The radius of the planet, R_p, and the inner radius, R_in, were both set to 1 R_Jup following Zhu et al. (2018). This choice does not affect our results significantly since, as we shall see for the CPD models, most of the flux does not come from the inner region.

The posterior distributions on the free parameters of our SED models were sampled with the Nautilus package (Lange 2023), which implements nested sampling aided by neuronal networks. We assumed flat priors in logarithm, with very wide boundaries.

We considered four models for the observed SED. Firstly, we considered only the dust component from the CPD. Secondly, we included an optically thin jet. Thirdly, we included an optically thick jet. Finally, we included free-free emission from the disc itself, using an α-constant disc and a magnetic disc. We also guided our analysis using a uniform slab model, adding the contribution of free electrons from metals in the magnetic disc.

3.2 Dust component of the CPD

In the absence of a jet, a strong degeneracy was observed between the parameters for smaller values of the accretion rate, Ṁ_p, and larger values of the viscosity parameters, α. For example, using synthetic data obtained with the parametric model of Zhu et al. (2018), for Bands 3, 4, 7, and 9 with Ṁ_p = 10⁻⁹ M_jup yr⁻¹ and α = 10⁻³, we obtained the orange corner plot shown in Fig. 5, revealing that the important parameter is the ratio between these two parameters. By contrast, in the case of a larger accretion rate, Ṁ_p = 10⁻⁵ M_jup yr⁻¹, and the same viscosity parameter, the flux depends solely on the accretion rate, as can be seen in the grey corner plot shown in Fig. 5 (as noted by Zhu et al. 2018).

The degeneracy between Ṁ_p and α can be understood by writing the emergent flux density, $$F_{\lambda }=\frac{1}{d^2}{\int }_{R_{\mathrm{i}\mathrm{n}}}^{R_{out}}{B}_{\nu }(T_{\mathrm{b}})2\pi R\mathrm{d}R,$$(8)

where B_ν is the Planck function (we use the full blackbody function rather than the Rayleigh–Jeans approximation because the disc could reach low temperatures), d is the distance to PDS 70, and the brightness temperature is given by Zhu et al. (2018), $$T_{\mathrm{b}}=\{\begin{array}{ll}{(\frac{3}{8}\frac{{\kappa }_{\mathrm{R}}}{{\kappa }_{\mathrm{m}\mathrm{m}}}{T}_{\mathrm{e}\mathrm{f}\mathrm{f}}^4+T_{\mathrm{e}\mathrm{x}\mathrm{t}}^4)}^{1/4}& \mathrm{i}\mathrm{f}{\tau }_{\mathrm{m}\mathrm{m}}>0.5,\\ \\ 2T_{\mathrm{c}}{\tau }_{\mathrm{m}\mathrm{m}}& \mathrm{i}\mathrm{f}{\tau }_{\mathrm{m}\mathrm{m}}<\mathrm{0.5.}\end{array}$$(9)

The Rosseland mean opacity was chosen to be constant, at the nominal value κ_R = 10 cm²/g (as in Zhu et al. 2018). The dust opacity is ${\kappa }_{\mathrm{m}\mathrm{m}}=3.4\zeta \times \frac{0.87\mathrm{m}\mathrm{m}}{\lambda }{\mathrm{c}\mathrm{m}}^2/\mathrm{g}$ (Andrews et al. 2012, with ζ the dust-to-gas ratio usually set as 0.01), and the optical depth τ_mm = κ_mmΣ/2. The effective temperature, T_eff, is the temperature that the surface of the viscous disc would radiate if it were a blackbody (Hartmann 1998; Zhu et al. 2018), $$T_{\mathrm{e}\mathrm{f}\mathrm{f}}^4=\frac{3G{M}_{\mathrm{p}}{\dot{M}}_{\mathrm{p}}}{8\pi \sigma R^3}[1-{\left(\frac{R_{\mathrm{i}\mathrm{n}}}{R}\right)}^{1/2}],$$(10)

and the external temperature is $T_{\mathrm{e}\mathrm{x}\mathrm{t}}^4=T_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{o}\mathrm{r}}^4+T_{\mathrm{i}\mathrm{r}\mathrm{r}}^4$, where T_floor is the background temperature in the environment of the CPD and T_irr is the temperature-equivalent irradiation flux from the planetary photospheric luminosity.

The background temperature in the circumstellar disc environment surrounding the CPD, T_floor, is a crucial parameter as it sets the floor temperature in the CPD. We used T_floor = 27K, which was obtained by extrapolating the dust physical conditions in the outer ring obtained by multi-frequency SED fits (these dust diagnostics will be documented in a forthcoming article). The azimuthally averaged temperature profile of 21K at a de-projected radius of 0.5″ (56 au) was extrapolated to the position of PDS 70c at 0.3″ (i.e. 34 au) under the assumption of a temperature profile of R^−1/2. In practice, the uniform-slab dust diagnostic suggest a much higher temperature, as a powerlaw fit to the dust temperatures in the outer ring yields 49.9 K at 34 au. However, we adopted the value of 27 K in a conservative approach, as a lower value widens the range of possible dust temperatures in the CPD. For comparison, the peak brightness temperature in Band 9, of 14.7 K at 67 au, would translate to 21 K at 34 au with a R^−1/2 temperature profile. Additionally, the mid-plane temperature estimates by Law et al. (2024) suggest a temperature of at least 30K at the distance of PDS 70c.

When the planetary accretion luminosity dominates over the central brightness of the planet we find that the irradiation temperature is $$T_{\mathrm{i}\mathrm{r}\mathrm{r}}={\left(\frac{L_{\mathrm{a}\mathrm{c}\mathrm{c}}}{\sigma 40\pi R^2}\right)}^{1/4},$$(11)

if the disc receives one tenth of the irradiation luminosity in the direction perpendicular to the disc (i.e. as in Zhu et al. 2018).

The optically thick and thin regimes in τ_mm can be distinguished from Eq. (9). In the optically thick regime, since $T_{\mathrm{e}\mathrm{f}\mathrm{f}}\propto {\dot{M}}_{\mathrm{p}}^{1/4}$ and $T_{\mathrm{e}\mathrm{x}\mathrm{t}}\propto {\dot{M}}_{\mathrm{p}}^{1/4}$ this implies that $T_{\mathrm{b}}\propto {\dot{M}}_{\mathrm{p}}^{1/4}$ (except for large radii where T_ext ~ T_floor). This results in a millimetre-flux density depending solely on the accretion rate, $F_{\nu }\propto {\dot{M}}_{\mathrm{p}}^{1/4}$.

In the optically thin regime, $T_{\mathrm{b}}\propto T_{\mathrm{c}}\mathrm{\Sigma }$. The mass column density, Σ, was calculated following Zhu et al. (2018), as the minimum between the equations $$\mathrm{\Sigma }=\frac{2^{7/5}}{3^{6/5}}{\left(\frac{\sigma G{M}_{\mathrm{p}}{\dot{M}}_{\mathrm{p}}^3}{{\alpha }^4{\pi }^3{\kappa }_{\mathrm{R}}{R}^3}\right)}^{1/5}\left(\frac{\mu }{\mathfrak{R}}\right){[1-{\left(\frac{R_{\mathrm{i}\mathrm{n}}}{R}\right)}^{1/2}]}^{3/5}\propto {\left(\frac{{\dot{M}}_{\mathrm{p}}^3}{{\alpha }^4}\right)}^{1/5},$$(12)

and $$\mathrm{\Sigma }=\frac{{\dot{M}}_{\mathrm{p}}\mu \mathrm{\Omega }}{3\pi \alpha T_{\mathrm{e}\mathrm{x}\mathrm{t}}}\propto {\left(\frac{{\dot{M}}_{\mathrm{p}}^3}{{\alpha }^4}\right)}^{1/4},$$(13)

where R is the ideal gas constant and μ is the mean molecular weight set as μ = 1. The mid-plane temperature is given by (Zhu et al. 2018): $$T_{\mathrm{c}}={(\frac{9G{M}_{\mathrm{p}}{\dot{M}}_{\mathrm{p}}\mathrm{\Sigma }{\kappa }_{\mathrm{R}}}{128\pi \sigma R^3}[1-{\left(\frac{R_{\mathrm{i}\mathrm{n}}}{R}\right)}^{1/2}]+T_{\mathrm{e}\mathrm{x}\mathrm{t}}^4)}^{1/4}.$$(14)

If we neglect the contribution of T_ext, which is for a purely viscously heated disc, then we keep Σ from Eq. (12), and Eq. (14) yields $$T_{\mathrm{c}}\mathrm{\Sigma }=\sqrt{\frac{G{M}_{\mathrm{p}}}{{\pi }^2{R}^3}}\frac{{\dot{M}}_{\mathrm{p}}}{\alpha }{\left(\frac{\mu }{\mathfrak{R}}\right)}^{5/4}[1-{\left(\frac{R_{\mathrm{i}\mathrm{n}}}{R}\right)}^{1/2}].$$(15)

Finally from Eq. (8) if the disc becomes optically thin at a sufficient small radii, which is the case with τ_mm < 0.5 of Eq. (9), T_b can be calculated using Eq. (15) since most of the contribution will come from larger radii (optically thin case).

The following simplifications can be made if the accretion rate is small enough. First, neglect the planetary accretion luminosity relative to the viscous heating mechanism with T_ext ~ 0. Second, approximate that the CPD is mostly optically thin by lowering the surface density, $\mathrm{\Sigma }\propto {\left({\dot{M}}_{\mathrm{p}}^3/{\alpha }^4\right)}^{1/5}$. A CPD heated by viscous dissipation allows for a simplified model of the SED, with optically thin millimetre emission that depends on a single parameter, $$\gamma =\frac{{\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r})}{\alpha }.$$(16)

Fig. 5 Corner plot for synthetic data. The orange posteriors show the synthetic optically thin case, calculated with Ṁp = 10−9 MJup yr−1 and α = 10−3. The outcome illustrates the degeneracy between Ṁp and α inherent to the parametric CPD model of Zhu et al. (2018). The grey posteriors show the synthetic optically thick case, calculated with Ṁp = 10−5 MJupyr−1 (and the same α = 10−3). Both synthetic SEDs were obtained using only the dust component of the CPD model following Zhu et al. (2018), with errors set to 10% of the flux.

3.3 HI free-free component of the CPD

Another component of the SED from PDS 70c is free-free emission from the disc itself. This scenario is supported by mid-plane temperatures exceeding 10⁴ K in the steady-state CPD models of Zhu et al. (2018). Also, the hydrogen recombination leading to the Hα signal from PDS 70c (Haffert et al. 2019) implies that the radio signal is, at some level, partly due to H I free-free continuum emission.

Here we consider the H I free-free luminosity of the parametric disc models from Zhu et al. (2018). The surface density, Σ(R) (Eqs. (12) and (13)), gives information on the gas density. Only the regions hot enough to ionise hydrogen (T ≳ 10 000 K) will contribute to the flux, and we neglect free electrons contributed by other elements.

Assuming that the gas follows a hydro-static and Gaussian vertical profile, $$\rho (R,z)=\frac{\mathrm{\Sigma }(R)}{\sqrt{2\pi }H(R)}\exp (-\frac{z^2}{2H(R{)}^2}),$$(17)

with the scale height H(R) = c_s∕Ω(R). We solved the radiative transfer equation with vertical stratification of the temperature (following Hubeny 1990), $$T(z{)}^4=\frac{3}{8}{T}_{\mathrm{e}\mathrm{f}\mathrm{f}}^4{\tau }_{\mathrm{R}}+T_{\mathrm{e}\mathrm{x}\mathrm{t}}^4,$$(18)

where τ_R is the Rosseland mean optical depth at z. We can calculate the number density of electrons using the Saha equation, $$\frac{n(\mathrm{H}\mathrm{I}\mathrm{I})}{n(\mathrm{H}\mathrm{I})}=\frac{2{\mathrm{Q}}_{\mathrm{I}\mathrm{I}}}{n_{\mathrm{e}}{\mathrm{Q}}_{\mathrm{I}}}{\left(\frac{2\pi m_{\mathrm{e}}{k}_{\mathrm{B}}T}{h^2}\right)}^{3/2}\exp (-\frac{\chi }{k_{\mathrm{B}}T}),$$(19)

where m_e is the electron mass, k_B is the Boltzmann constant, χ is the collisional ionisation potential of hydrogen (13.6 eV), h is the Planck constant, and the partition function was calculated as $${\mathrm{Q}}_{\mathrm{I}}=\sum _{n=1}^{3}2n^2\exp (-\frac{\chi (1-1/n^2)}{k_{\mathrm{B}}T}),$$(20)

for neutral hydrogen, and Q_II = 1 for ionized hydrogen. The ionisation fraction of an initially neutral gas is given by f = n_e∕n_H, where n_H = ρ(R, z)∕m_H (we assume the gas consists of hydrogen only). We can then solve for n_e, $$\begin{array}{r}{n}_{\mathrm{e}}(R,z)=\frac{{\mathrm{Q}}_{\mathrm{I}\mathrm{I}}}{{\mathrm{Q}}_{\mathrm{I}}}{\left(\frac{2\pi m_{\mathrm{e}}{k}_{\mathrm{B}}T}{h^2}\right)}^{3/2}\exp \left(\frac{-\chi }{k_{\mathrm{B}}T}\right)\\ (-1+\sqrt{1+2n_{\mathrm{H}}\frac{{\mathrm{Q}}_{\mathrm{I}}}{{\mathrm{Q}}_{\mathrm{I}\mathrm{I}}}{\left(\frac{2\pi m_{\mathrm{e}}{k}_{\mathrm{B}}T}{h^2}\right)}^{-3/2}\exp \left(\frac{\chi }{k_{\mathrm{B}}T}\right)}).\end{array}$$(21)

With the number density of electrons we calculated the emission measure (EM), $$\mathrm{E}\mathrm{M}(R,z)={\int }_z^{\mathrm{\infty }}{n}_{\mathrm{e}}(R,z)n_{\mathrm{i}}(R,z)dz,$$(22)

where n_i is number density of positive ions (for this initially neutral gas n_i = n_e). We also calculated the optical depth (Wright & Barlow 1975; Allen 1973, including stimulated emission), $${\tau }_{\nu }(R,z)={\int }_z^{\mathrm{\infty }}{K}_{\nu }(T)n_{\mathrm{e}}(R,z)n_{\mathrm{i}}(R,z)dz,$$(23)

with $$\begin{array}{r}{K}_{\nu }(T)=\frac{4\pi }{3\sqrt{3}}\frac{Z^2{e}^6}{hc{m}_{\mathrm{e}}^2\upsilon }\frac{g_{\nu }}{{\nu }^3}(1-\exp (-\frac{h\nu }{k_{\mathrm{B}}T}))\\ =\frac{3.692\times {10}^8}{\mathrm{c}{\mathrm{m}}^{-5}{\mathrm{K}}^{-1/2}\mathrm{H}{\mathrm{z}}^{-3}}(1-\exp (-\frac{h\nu }{k_{\mathrm{B}}T}))Z^2{g}_{\nu }{T}^{-1/2}{\nu }^{-3},\end{array}$$(24)

where the electron velocity υ was replaced by the mean velocity $\overline{\upsilon }=(\pi kT/2m_{\mathrm{e}}{)}^{1/2}$, Z is the ionic charge (in this case Z = 1), e is the electron charge, and c is the velocity of light in vacuum. The gaunt factor, g_ν, was approximated using the equation from Oster (1961), $$\begin{array}{r}{g}_{\nu }(T)=\ln \left({\left(\frac{2k_{\mathrm{B}}T}{\gamma m_{\mathrm{e}}}\right)}^{3/2}\left(\frac{2m_{\mathrm{e}}}{\gamma Z{e}^2\omega }\right)\right)\\ =\ln (4.955\times {10}^7{\left(\frac{\nu }{\mathrm{H}\mathrm{z}}\right)}^{-1})+1.5\ln \left(\frac{T}{\mathrm{K}}\right),\end{array}$$(25)

with $\gamma =e^{{\gamma }^{\ast }}\approx 1.781$ the exponential of Euler’s constant γ^* and ω = 2πν the angular frequency. When g_ν(T) < 1, we set g_ν(T) = 1, which is the limit at higher frequencies¹.

Finally, we calculated the free-free emission from the CPD by integrating over the disc, $$F_{\nu }={\int }_{R_{\mathrm{i}\mathrm{n}}}^{R_{\mathrm{o}\mathrm{u}\mathrm{t}}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-({\tau }_{\nu }-{\tau }_{\nu }^{\prime })}{B}_{\nu }(T(z^{\prime }))K_{\nu }(T(z^{\prime }))n_{\mathrm{e}}^2(R,z^{\prime })d{z}^{\prime }\frac{2\pi R}{d^2}dR,$$(26)

where we have already replaced n_i = n_e, B_ν(T) is the Planck function, τ_ν = τ_ν(R, z → -∞), is the total optical depth from Eq. (23), and ${\tau }_{\nu }^{\prime }={\tau }_{\nu }(\mathrm{R},{\mathrm{z}}^{\prime })$, is the optical depth calculated using Eq. (23), integrating from z′ to ∞. The distance to PDS 70 is d.

3.4 Uniform slab model

In parallel to the calculation of the free-free component from the analytic viscous disc model, we also implemented a simpler treatment in order to inform ourselves about the EMs and disc sizes that would match the observed SED. We used a pill-box model consisting of a uniform slab with a constant temperature and EM, and a cross-section of $A=\pi R_{\mathrm{m}\mathrm{a}\mathrm{x}}^2$. The optical depth given by Eq. (23) reduces to $${\tau }_{\nu }^{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{l}}=K_{\nu }(T)\times \mathrm{E}\mathrm{M},$$(22)

then from Eq. (26), since T and EM are constants, we solved both integrals analytically using R_in = 0 and R_out = R_max and obtained $$F_{\nu }^{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{l}}=\pi B_{\nu }(T)(1-e^{-{\tau }_{\nu }^{\mathrm{p}\mathrm{i}\mathrm{l}\mathrm{l}}}){\left(\frac{R_{\mathrm{m}\mathrm{a}\mathrm{x}}}{d}\right)}^2.$$(28)

3.5 Magnetic disc

The previous models did not take into account the free electrons coming from metals, since the viscous disc models in consideration reach hydrogen-ionising temperatures. However, the very high effective viscosity in a strongly magnetised disc would result in lower peak temperatures, as in the parametric disc models from Hasegawa et al. (2021). Such ‘magnetic disc’ models assume that the planet is undergoing magnetospheric accretion, which means that the inner region of the CPD has enough free electrons to couple the gas to the magnetic field of the planet.

We followed Hasegawa et al. (2021), with a magnetic field (B_p) described as a dipole, $$B_{\mathrm{p}}=B_{\mathrm{p}\mathrm{s}}{\left(\frac{R_{\mathrm{p}}}{R}\right)}^3,$$(29)

where B_ps is the field strength at the surface of the planet. We used an alpha-viscosity (α_m) inspired by the results of the ideal magnetohydrodynamics (MHD) simulations of Salvesen et al. (2016), $${\alpha }_{\mathrm{m}}=11{\beta }^{-0.53},$$(30)

where $\beta ={\rho }_{\circ }{c}_{\mathrm{s}}^2/(B_{\mathrm{p}}^2/8\pi )$ and ${\rho }_{\circ }=\mathrm{\Sigma }\mathrm{\Omega }/(\sqrt{2\pi }{c}_{\mathrm{s}})$ is the mid-plane gas volume density.

Since the dusty disc model from Sect. 3.2 already has the disc surface density, Σ(R), and the temperature profile, T_c(R), for a given α and Ṁ_p, we can solve Eq. (30) by imposing α = α_m for a fixed B_ps and Ṁ_p. Thus, the magnetic disc model has a variable viscous parameter α(R), with the rest of the physical conditions calculated analogously to Sect. 3.2. However, these viscosity estimates neglect other sources of ionisation (i.e. interstellar UV and cosmic rays), and α values can be unrealistically low at large radii. We therefore imposed a lower limit, α(R) > 10⁻⁵.

With the temperature and density profiles as a function of Ṁ_p and B_ps, the next step was to consider the free electrons coming from the metals. For this we used the abundances and collisional ionisation potentials of the most relevant elements in the solar photosphere, shown in Table 5 (Asplund et al. 2009; Lide 2008; Keith & Wardle 2014).

We used these abundances and solved the Saha equation for each element, imposing the neutrality condition $$n_{\mathrm{e}}=\sum _i{n}_{\mathrm{i}},$$(31)

where the sum runs over the five elements in consideration, and n_i is the number density of the i-ionised element. We assumed that the metals are not depleted into dust. However, solving this equation takes too long for an optimisation. We decided to use a simpler approach that takes advantage of the differences in the abundance of each element. We assume that at a given temperature a single element contributes most of the free electrons. We solved the Saha equation separately, setting n_e = n_i for each element, and used the abundances to determine which element provided most of the electrons at a given temperature (a similar approach can be found in Mohanty et al. 2018). This allows for an approximation of the real number density of electrons without requiring an excessively long computing time.

Following Hasegawa et al. (2021), we take into account that the CPD is undergoing magnetospheric accretion. This means that the radius at which the planetary magnetic field dominates over the ram pressure (called the truncation radius, R_T) is larger than the radius at which the disc has enough free electrons to couple the gas to the magnetic field. This condition is satisfied for T ~ 1000 K, which is the temperature threshold for K, and we used the condition R_T > R(T_c = 1000 K) to ensure that the model is consistent (Hasegawa et al. 2021). Finally, we calculated the optical depth due to metal free-free using the electron number density in Eq. (23).

3.6 Negative hydrogen ion contribution

The optical/IR opacity of the sun and cooler stars is dominated by the negative hydrogen ion H⁻ (Wildt 1939), and it is worth exploring how it contributes to the CPD radio opacity. Because of its simple and highly polarised structure, the hydrogen atom can hold a second electron with an ionization potential of 0.7542 eV (Frolov 2015).

We included Wildt’s photo-detachment, which is the bound-free mechanism, $$h\nu +{\mathrm{H}}^{-}\to \mathrm{H}+{\mathrm{e}}^{-},$$(32)

where any photon with λ < 16 440 Å had sufficient energy to ionise the H⁻. We also included Pannekoek’s free-free mechanism (Pannekoek 1931), $$h\nu +{\mathrm{e}}^{-}+\mathrm{H}\to \mathrm{H}+{\mathrm{e}}^{-},$$(33)

which contributes to the opacity, especially at longer wavelengths.

The presence of ionised metals is key to the abundance of H⁻, since we need extra electrons to encounter neutral hydrogen, and H⁻ also draws free electrons from the medium. This means we have to solve Saha’s equation, again adding the H⁻. We have its ionisation energy (0.7542 eV), partition function, ${\mathrm{Q}}_{{\mathrm{H}}^{-}}\sim 1$, Q_H ~ 2, and a new neutrality condition, $$n_{\mathrm{e}}+n_{{\mathrm{H}}^{-}}=\sum _i{n}_{\mathrm{i}},$$(34)

which follows from Eq. (31). The difference is that some of the free electrons from the ionised metals are now bound to the H⁻.

Here again we simplified the solution of the Saha equation by taking advantage of the different abundances of each element. In most cases, the abundance of H⁻ was negligible compared with the ions, and we solved Saha’s equation as shown in Sect. 3.5. We calculated the number density of H⁻ with the number density of electrons (n_e) obtained with metals only. In the case in which the number density of H⁻ becomes significant (10% of the ions or more), the free electrons are less abundant in the medium. The number density of free electrons decreases drastically and we set n_i = n_H⁻ = n_e = 0.

Provided with the number density of electrons, we calculated the bound-free and free-free opacities of H⁻ following Gray (2022, their Chap. 8). The final total optical depth is $${\tau }_{\nu }^{\mathrm{t}\mathrm{o}\mathrm{t}}={\tau }_{\nu }^{\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{s}}+{\tau }_{\nu }^{{{\mathrm{H}}^{-}}_{\mathrm{f}\mathrm{f}}}+{\tau }_{\nu }^{{{\mathrm{H}}^{-}}_{\mathrm{b}\mathrm{f}}},$$(35)

with ${\tau }_{\nu }^{\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{s}}$ the optical depth calculated as Sect. 3.5, and ${\tau }_{\nu }^{{{\mathrm{H}}^{-}}_{\mathrm{b}\mathrm{f}}}$ and ${\tau }_{\nu }^{{{\mathrm{H}}^{-}}_{\mathrm{f}\mathrm{f}}}$ the bound-free and free-free optical depth of H⁻, respectively. The total flux was then calculated as $$F_{\nu }={\int }_{R_{\mathrm{i}\mathrm{n}}}^{R_{\mathrm{o}\mathrm{u}\mathrm{t}}}{\int }_0^{{\tau }_{\nu }}{e}^{-({\tau }_{\nu }-{\tau }_{\nu }^{\prime })}{B}_{\nu }({\tau }_{\nu }^{\prime })d{\tau }_{\nu }^{\prime }\frac{2\pi R}{d^2}dR,$$(36)

where τ_ν is the ${\tau }_{\nu }^{\mathrm{t}\mathrm{o}\mathrm{t}}$ calculated using Eq. (35).

Fig. 6 Best fits for the data of Bands 3, 4, 7, and 9 using the dusty disc (Zhu et al. 2018) model and the simple dusty disc model (described in Sect. 3.2) are shown as the solid orange and light grey lines, respectively. For the dusty disc, we find that the best fit parameters of the model are $\log ({\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r}))=-{9.04}_{-2.75}^{+0.46}$ and $\log (\alpha )=-{5.21}_{-2.82}^{+0.53}$ at 1 σ. The resulting spectral index is 2.35 and has a minimum reduced χ2 of 6.50. For the simpler dusty disc model, we find that the best fit parameter is $\log (\gamma )=-{2.67}_{-0.05}^{+0.04}$. The resulting spectral index is 3 and has a minimum reduced χ2 of 14.61. Here the reduced χ2 is the χ2 divided by the number of points (data) minus the number of degrees of freedom in the model.

3.7 Fitting the data

We first fitted the nearly coeval data in Bands 3, 4, 7, and 9 using the dusty CPD model from Sect. 3.2, varying both parameters (α and Ṁ_p). We observed that the accretion rate and viscosity parameters were degenerate, and therefore applied the simpler model that only depends on γ (Eq. (16)). The best fit of both models is shown in Fig. 6, which poorly fits the data.

For the dusty disc model shown in Fig. 6, we observed that the bulk of the flux originates from the outer radius, where the mid-plane temperature is T_c = T_ext ≈ T_floor = 27 K. This indicates that the SED is highly dependent on the floor temperature, and our approximation using the single parameter, γ, was invalid (T_ext cannot be neglected). As is shown in Fig. 7, if the floor temperature were lower (i.e. T < 15 K), we could fit the data.

The addition of an optically thin and thick jet to the dusty disc model was performed. We fixed α = 10⁻² and also varied the dust/gas ratio, ζ, using uniform priors with log(γ) and log(ζ) between −15 and 15. However, in both cases the contribution of the jet after the fit was negligible and it just returned the dusty disc case without a jet.

We then turned to the new model that sums the free-free component (from Sect. 3.3) and the thermal dust emission from the disc. Since there is no certainty that the approximation made in Sect. 3.2 is valid (with the single parameter, γ), we used the complete model proposed by Zhu et al. (2018). We considered two different models: fixed α, for which we kept a fixed alpha value of α = 10 ² and allowed the accretion rate to vary, and fixed Ṁ_p, for which we kept the accretion rate fixed at Ṁ_p = 10⁻⁶ M_Jup yr⁻¹ and allowed the alpha value to vary. Fig. 8 shows the resulting best fits for each case. The models were able to produce spectral indices of $({\alpha }_{97.5}^{145},{\alpha }_{145}^{343.5},{\alpha }_{343.5}^{671})=(1.87,1.86,1.97)$ and $({\alpha }_{97.5}^{145},{\alpha }_{145}^{343.5},{\alpha }_{343.5}^{671})=(2.02,1.93,1.83)$, respectively, which are more consistent with Table 4 for Bands 3 to 7. However, these models only fit the data points for Bands 3, 4, and 7 and do not account for the non-detection in Band 9 with a flux, F₆₇₁, of 417 µJy and 389 µJy in each case.

In the model with fixed α we derive a well-constrained mass accretion rate, and vice versa. However, the results also indicate that free-free emission is the dominant component of the CPD signal. The inferred upper limits on the dust-to-gas ratio are 10^−4.80 and 10^−7.05 for the models with fixed α and fixed Ṁ_p, respectively. This suggests a strongly dust-depleted environment, which could point to reduced dust transport into the CPD, with larger grains being trapped in a pressure bump outside the planet’s orbit (e.g. Zhu et al. 2012). Weber et al. (2018) show that in the presence of viscosity, even with strong gap-opening planetary torques, there is a flow of gas and micron-sized grains through the planet-induced gap, while grains larger than a certain size are filtered out at the gap edge. In contrast, a study more focussed on dust transport to the CPD by Szulágyi et al. (2022) reports a contradictory result: the planet’s gravitational perturbation elevates millimetre-sized particles from the midplane and induces a vertical accretion flow from the circumstellar disc towards the CPD. Since Weber et al. (2018) considered a two-dimensional case over long evolutionary timescales – leaving the planet sufficient time to clear its orbit – whereas Szulágyi et al. (2022) conducted three-dimensional simulations evolved for only 200 orbits, it remains an open question how dust filtration and transport operate in more evolved, three-dimensional cases.

The model with a fixed α = 10⁻² requires very high accretion rates, Ṁ_p ~ 10^−3.5 M_Jup yr⁻¹. However, by reducing the viscosity parameter, for example to α = 10⁻³, the flux ends up being overestimated by more than a factor of two. This means that reducing the viscosity parameter would require a lower accretion rate to fit the data. Using the same example, for α = 10⁻³ the accretion rate that fits the data is Ṁ_p ∼ 10⁻⁴ M_Jup yr⁻¹.

The properties of the disc are listed in Figs. 9 and 10 for the case with fixed α and fixed Ṁ_p, respectively. At small radii (R ≲ 0.02au) the fraction of ionisation is f ~ 1, and the resulting EMs are extremely large (EM ~ 10^26-30 cm⁻⁶ pc). The profile for the surface density (Σ) is in the viscous-disc regime (small T_ext), and in the case with fixed Ṁ_p and small α (i.e. α ~ 10⁻⁷), Σ is two orders of magnitude higher than for large values of α (i.e. α ~ 0.01). Another outstanding feature of the physical conditions in these analytical viscously accreting CPDs is the occurrence of extremely hot inner discs, for any α < 0.01. Such hydrogen-ionising temperatures are counter-intuitive for a CPD.

For the uniform slab model, we fixed a temperature of 10³ K and 10⁴ K and chose uniform priors for log (R_max/au) and log (EM/(cm⁻⁶ pc)) (with very wide boundaries). The resulting best fits are shown in Fig. 11. Now we could produce spectral indices of $({\alpha }_{97.5}^{145},{\alpha }_{145}^{343.5},{\alpha }_{343.5}^{671})=(2.00,1.91,0.59)$ and $({\alpha }_{97.5}^{145},{\alpha }_{145}^{343.5},{\alpha }_{343.5}^{671})=(2.00,1.88,0.81)$, respectively, which are more in line with Table 4 even for Band 9.

We notice that the values for the EM are smaller compared with the parametric model of Zhu et al. (2018), having the best fit EM = 1.19 × 10¹¹ cm⁻⁶ pc and EM = 1.04 × 10¹² cm⁻⁶ pc for each case. On the other hand, the radius of the pill seems to match the order of magnitude of the radius where we start having a smaller ionisation fraction in Figs. 9 and 10, with a best fit size of R_max = 0.078 au and R_max = 0.025 au for the pill in each case.

For the magnetic disc, we chose uniform priors for log(B_ps/G) between −4 and 4 and log $\log ({\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r}))$ between –10 and –2. We also added a penalty in the log likelihood for the magnetospheric accretion criterion by subtracting the quantity $x=10\times \frac{|R_{\mathrm{T}}-R(T_{\mathrm{c}}=1000\mathrm{K})|}{R(T_{\mathrm{c}}=1000\mathrm{K})}$ for the cases in which R_T < R(T_c = 1000 K). The resulting best fit is shown in Fig. 12. Now we could produce spectral indices of $({\alpha }_{97.5}^{145},{\alpha }_{145}^{343.5},{\alpha }_{343.5}^{671})=(1.85,1.73,1.76)$ that could marginally fit the data, and overestimate Band 9 with a flux F₆₇₁ of 346 µJy.

From Fig. 13 we can see that the model satisfies the magnetospheric accretion criterion (R_T > R(T_c = 1000 K)). The radius of the peak EM, with EM ~ 10^11.5 cm⁻⁶ pc, is R ~ 0.02 au, and both EM and R approximate those required by the uniform slab model. The size of the magnetic disc model is smaller than the 1/10 the Hill radius, which is reminiscent of the upper limit to the extension of the sub-millimetre signal observed in SR 12c (Wu et al. 2022).

Because of the residual ionisation stemming from collisional ionisation of metals, part of the atomic phase of the CPD can be seen as an atomic plasma. This state of matter is characterised by temperatures of T ~ 1000-2000 K. The free electrons released by the metals interact with both neutral hydrogen and ionised metals, and thus produce the radio signal.

Atomic plasma radiation in the magnetic disc model can marginally fit the data. However, an optically thick spectrum persists at higher frequencies and the drop in Band 9 cannot be reproduced. This is because, even though the magnetic disc reaches the necessary EMs for free-free emission from ionised metals, H⁻ free-free also contributes to the total optical depth. In Fig. 14 the optical depth of metals behaves exactly as we would like it to, going from optically thick to thin from 97.5 GHz to ∼500 GHz. However, the optical depth of H⁻ begins to dominate above 500 GHz, and the total optical depth remains τ ~ 1 without reaching the optically thin regime.

The atomic plasma model, however, depends heavily on the accuracy of the extrapolation from the H⁻ opacities, specially the free-free opacity that dominates in the sub-millimetre continuum. Gray (2022) was inspired by works that calculated opacities of up to 10 µm and in a range of temperatures from 2500 to 10 000 K. This is very important here, since we are calculating the opacity at millimetre wavelengths, and in a cold disc where the majority of the flux comes from temperatures of T ≈ 1000–2000 K.

Fig. 7 Best fit for the data of Bands 3, 4, 7, and 9 for the dusty disc (Zhu et al. 2018) using different Tfloor. In order, the solid orange, light grey, and grey lines correspond to Tfloor of 10, 15, and 20 K.

Fig. 8 Best fits to the data from Bands 3, 4, 7, and 9 using two free-free emission models (sum the dusty disc component from Sect. 3.2 and the free-free component from Sect. 3.3). The first has a fixed α = 10−2, shown as the solid orange line and has a minimum reduced χ2 of 3.18. We find that the best fit parameters of this model are $\log ({\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r}))=-{3.48}_{-0.09}^{+0.06}$ at 1 σ, with a strict upper limit on the dust-to-gas ratio of log(ζ) = −4.80. The second has a fixed Ṁp = 10−6 MJup yr−1 as the solid light grey line and has a minimum reduced χ2 of 2.54. We find that the best fit parameters for this model are $\log (\alpha )=-{6.86}_{-0.07}^{+0.08}$ at 1 σ, with a strict upper limit on the dust-to-gas ratio of log(ζ) = −7.05.

Fig. 9 Physical conditions including free-free in the best-fit CPD model with α fixed. Tc is the mid-plane temperature calculated as Eq. (14), Σ(R) is the surface density, EM is the emission measure, f is the ionisation fraction in the mid-plane, and Hhwhm is the height at which the density of free electrons is half the density of free electrons in the midplane (half-width at half-maximum).

Fig. 10 Physical conditions in the best-fit CPD model with Ṁp fixed. Annotations follow from Fig. 9.

Fig. 11 Best fit to the data from Bands 3, 4, 7, and 9 using the uniform slab model described in Sect. 3.4. On the upper right panel, the solid orange line corresponds to the case T = 103 K, has a minimum reduced χ2 of 0.18, and a p value of 83%. We find that the best fit parameters of this model are $\log (\mathrm{E}\mathrm{M}/(\mathrm{c}{\mathrm{m}}^{-6}\mathrm{p}\mathrm{c}))={11.08}_{-0.19}^{+0.53}$ and $\log \left(R_{\mathrm{m}\mathrm{a}\mathrm{x}}/\mathrm{a}\mathrm{u}\right)=-{1.11}_{-0.06}^{+0.03}$. The solid light grey line corresponds to the case T = 104 K, has a minimum reduced χ2 of 0.33, and a p value of 72%. We find that the best fit parameters of this model are $\log (\mathrm{E}\mathrm{M}/(\mathrm{c}{\mathrm{m}}^{-6}\mathrm{p}\mathrm{c}))={12.00}_{-0.18}^{+0.42}$ and $\log \left(R_{\mathrm{m}\mathrm{a}\mathrm{x}}/\mathrm{a}\mathrm{u}\right)=-{1.61}_{-0.04}^{+0.02}$. The bottom left panel shows the corner plots obtained for each model. The orange and grey posteriors also correspond to the cases of T = 103 K and T = 104 K, respectively.

4 Discussion

4.1 Which models fit the data

The best fits obtained for the dusty disc model, shown in Fig. 6, have spectral indices of α > 2.35 and do not fit the data, even between Bands 3 and 7. Adding a jet, either optically thin or thick, does not change this result. However, we observe that the bulk of the flux stems from the outer regions of the disc. If the temperature is low enough (T_floor < 15), then a dusty disc can fit the data, as is shown in Fig. 7. This can be understood from the Wien displacement law, which states that at such low temperatures the turnover of the black-body spectrum is close to the Band 9 frequency. However, these low temperatures are unlikely, as is mentioned in Sect. 3.2. Band 9 shows a lower temperature limit of 21 K at 34 au, and the presence of CO in the PDS 70 cavity suggests temperatures of at least 30 K (Law et al. 2024). Therefore, we conclude that it is not possible to fit the observed SED using only a dusty disc model with the current derivations of T_floor.

If viscous heating is efficient and reaches H-ionising temperatures, our best fit with H I free-free from the H II disc reaches EM ∼ 10²⁶⁻³⁰ cm⁻⁶ pc (see Figs. 9 and 10), which seems overwhelming and would approximate the radio spectrum of an early-type star. This results in an optically thick SED that overshoots the flux density in Band 9.

A uniform slab model can inform one about the parameters needed to produce the Band 9 drop, with an EM of ∼10¹¹⁻¹² cm⁻⁶ pc, temperatures of T ∼ 10³–10⁴ K, and sizes of ∼0.1 au. Comparing these EMs and temperatures with those from the H II disc (see the implicit dependence of EM and T in Figs. 9 and 10), we see that such EMs in an HII CPD would require low ionisation fractions, and hence fine-tuning of the temperature. Since f(T) is extremely sensitive to T, we conclude that the CPD is probably not hot enough to ionise hydrogen, as, unless the temperature field is ‘just right’, hydrogen-ionisation results in huge EM and a very optically thick SED (and inconsistent, at 2.6 σ, with the non-detection in Band 9).

Finally, if the planet has a significant magnetic field (B_ps > 200 G), magnetospheric accretion occurs, the temperature drops to ∼10³ K, and the signal stems from free-free emission from ionised metals at lower frequencies (ν ≲ 500 GHz) and from H⁻ at higher frequencies (ν > 500 GHz). This atomic plasma radiation model maintains a partially optically thick SED above 500 GHz and also overshoots Band 9.

Fig. 12 Results using the magnetic disc model described in Sect. 3.5. This takes into account the free-free from metals and H− and the bound-free from H− described in Sect. 3.6. On the upper right, we have the best fit for the data of Bands 3, 4, 7, and 9 represented with their error bars using the fluxes of Table 3. The modelled flux is represented by the solid line, obtaining a partially thick free-free emission with a spectral index of 1.76. At the bottom left, we have the corner plot with the best fit parameters $\log ({\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r}))=-{5.37}_{-0.12}^{+0.06}$ and $B_{\mathrm{p}\mathrm{s}}={610}_{-375}^{+98}\mathrm{G}$ at 1 σ. The minimum reduced χ2 is 2.37.

Fig. 13 Physical conditions including metal free-free emission in the best-fit MHD model. α is the viscosity parameter, Tc is the mid-plane temperature calculated according to Eq. (14), Σ(R) is the surface density, EM is the emission measure, f is the ionisation fraction in the midplane, and Hhwhm is the height at which the density of free electrons is half the density of free electrons in the mid-plane (half-width at halfmaximum). The vertical dashed red line is the truncation radius, Rt = 1.69RJup, and the disc total mass is 1.93MJup.

Fig. 14 Optical depth of H− (including free-free and bound-free) and metals, in the annulus that makes the greatest contribution to the flux for the magnetic disc model (Sect. 3.5). The solid blue and orange lines correspond to the optical depth by H− and metals, respectively. The dashed line is the total optical depth.

4.2 Surface shock on the CPD

The small EMs from the uniform-slab models are reminiscent of shock models, as was studied by Aoyama et al. (2018) and Aoyama et al. (2020) for the CPD surface shock and Aoyama et al. (2020) for the photospheric shock. These shocks arise from the supersonic, nearly vertical gas flows that collide with the CPD or planetary surfaces, as was predicted in detailed hydro-dynamical simulations Tanigawa et al. (2012); Szulágyi & Ercolano (2020). The strong shock heats the gas and results in a thin, non-equilibrium layer of ionised hydrogen.

Because of the high temperatures (T ≳ 10⁴ K) in the postshock gas, we can neglect contributions from H⁻ to the net opacity. The bulk of the mass is at a relative constant temperature of ~10⁴ K. Since the most emission originates from the region with T ~ 10⁴ K, we refer to 10⁴ K as the effective temperature for free-free emission (Aoyama et al. 2018, 2020). As is seen in the uniform slab model of Sect. 3.4 for temperatures of 10⁴ K, an effective radius of $\sim {52}_{-5}^{+5}{R}_{\mathrm{J}\mathrm{u}\mathrm{p}}$ is required to fit the SED. This indicates that the flux in the radio should come from shocks on the CPD and neither from the photosphere nor accretion funnels (whose sizes are ~2R_Jup).

To firmly rule out photospheric shocks and accretion funnels, we estimated the temperature and EM needed for a fixed size of R_max = 2 R_Jup in the uniform slab model. The best fit parameters are $\log \left(T/\mathrm{K}\right)={6.83}_{-0.10}^{+0.04}$ and $\log (\mathrm{E}\mathrm{M}/(\mathrm{c}{\mathrm{m}}^{-6}\mathrm{p}\mathrm{c}))={15.75}_{-0.14}^{+10.11}$. Such temperatures are only achieved in the post-shock region of an object > 200 M_Jup with a 2 R_Jup radius, which is outside of the planetary regime. High temperatures, in excess of ~10⁶ K, also quench Hα radiation as recombination becomes less probable. Therefore a photospheric shock, or free-free emission from the accretion funnels, cannot be the source of the radio signal as such shocks would require excessively high temperatures, given the free-fall velocities in the planetary regime. Note that part of the Hα signal from PDS 70c could still stem from a T ~ 10⁴ K shock on the planetary surface, even if such a shock does not account for the radio signal.

To confirm that the post-shock emission on the CPD surface can indeed reproduce the observed radio emission, we performed two simulations of shock-heated gas following Aoyama et al. (2018), whose results are shown in Figs. 15 and 16. These correspond to pre-shock velocities and neutral hydrogen number densities of (υ₀, n₀) = (30kms⁻¹, 10¹⁵cm⁻³) and (50kms⁻¹, 10¹² cm⁻³), respectively. In both cases, the emission at 300 GHz is optically thick (τ > 1) but turns over into an optically thin spectrum at 600 GHz (τ ~ a few times 0.1). This is also the case for the SED of PDS 70c obtained with the uniform slab model, with a best fit of τ_{300 GHZ} ≈ 3 and τ_600GHZ ≈ 0.6.

In these simulations, the optical-depth weighted temperatures are ~7000 and ~7500 K, for υ₀ = 30kms⁻¹ and υ₀ = 50kms⁻¹, respectively. A uniform slab at 7000 K requires a radius of R_max = 62 R_Jup, which is in between the values obtained for the uniform slab models at 10³ K and 10⁴ K shown in Fig 11. We can estimate the planetary mass required to achieve a preshock velocity of υ₀ > 30kms⁻¹ at 62 R_jup, by assuming that υ₀ is the free-fall velocity, $${\upsilon }_0=\sqrt{\frac{2G{M}_{\mathrm{p}}}{r}}.$$(37)

A minimum planetary mass of 15.7 M_Jup is required, which is 2.1 M_jup above the 2σ upper limit of 13.6 M_jup estimated for the dynamical mass of PDS 70c (Trevascus et al. 2025). This minimum planetary mass could be reduced by considering the increasing υ₀ at smaller radii, and its resulting radio flux, but this is beyond the scope of this work.

Using a mean molecular weight of μ = 1.35, v₀ = 30kms⁻¹, and n₀ = 10¹⁵ cc⁻¹ implies a mass flux of μm_h × n₀ × υ₀ = 6.75 × 10⁻³ cm⁻² s⁻¹ g. If this flux is is spatially constant and falls on both sides of a circular slab with a radius of 62 R_Jup, the mass influx rate within this area is 0.1 M_Jup yr⁻¹. For υ₀ = 50 kms⁻¹ and n₀ = 10¹² cc⁻¹, the crude estimate of the mass influx rate is instead 2 × 10⁻⁴ M_Jup yr⁻¹, for simplicity using the same effective radius of 62 R_Jup. These rough estimates of the accretion rate may seem extremely high. However, not all of the mass falling on the CPD is necessarily accreted by the planet because of equatorial decretion and recycling of the CPD mass along the planetary wakes (Szulágyi et al. 2016; Batygin 2018).

Moreover, while we assumed a constant mass flux across the entire CPD, as in Aoyama et al. (2018), a concentration of mass flux towards the inner region could significantly lower the mass-accretion-rate estimate, while preserving the radio flux. The larger υ₀ and radio intensity at smaller radii could produce the same radio flux with a smaller mass accretion rate. The radial distribution of CPD mass-loading is, however, an open question, with some studies favouring the outer regions of CPDs (Tanigawa et al. 2012; Marleau et al. 2024).

If the observed radio signal from PDS 70c is entirely due to H I free-free emission from the CPD surface shock, the emission from the CPD itself must be faint. This places an upper limit on the planetary accretion rate. Using the magnetic disc model from Sect. 3.5 with the best fit value of B_ps = 610G, we obtain an upper limit of $\log ({\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r}))<-5.42$ using the 3 σ flux limit in Band 7. In the case of a smaller magnetic field (for example B_ps = 200 G), the upper limit for the accretion rate diminishes to $\log ({\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r}))<-5.77$.

The non-detection of PDS 70b can be understood in the context of the CPD-surface shock if PDS 70b is less massive than PDS 70c (e.g. Trevascus et al. 2025). In this case, PDS 70b would have the same υ₀ at a smaller radii. As an illustration, for a 3 M_Jup planet, the pre-shock velocity υ₀ = 30 km s⁻¹ is obtained at R_shock,b = 12 R_Jup. This means that the radii for the shocked region where free-free emission arises is smaller and consequently the total flux decreases. As a rough approximation, the radio flux from PDS 70b would be lower than PDS 70c by a factor corresponding to the ratio of the shocked surface area, or ~(R_max/R_shock,b) ~ (62/12)² ~ 30.

Fig. 15 Physical conditions of the post-shock region for υ0 = 30 km s−1 and n0 = 1015cc−1. The temperature, free electron local density (ne−), and HII local density (nH+) are the solid red, black, and blue lines, respectively. The optical depths at 300 and 600 GHz (τ300GHz and τ600GHz) are the dashed green and yellow lines, where the dotted black line corresponds to τ = 1.

Fig. 16 Physical conditions of the post-shock region for v0 = 50 km s−1 and n0 = 1012cc−1. Annotations follow from Fig. 15.

5 Conclusion

We report on new ALMA observations of PDS 70, focussing on the point-source signal coincident with PDS 70c, with new detections in Bands 4 and 7, and a marginal detection in Band 3. The spectrum across Bands 3, 4, and 7 is consistent with optically thick emission, with α = 2.0 ± 0.2. However, a nondetection in Band 9 breaks this trend, with a flux density lower than an optically thick extrapolation of lower frequencies by at least 2.6 σ. The Band 7 signal from PDS 70c appears to be constant within 10% rms, from May to December 2023. Since the multi-frequency measurements are coeval within two months, we assume that the SED collected here is constant.

The non-detection in Band 9 is particularly surprising, given the model predictions for bright signal at short sub-millimetre wavelengths, and suggest a very dust poor disc, whose radio signal stems from partially thick free-free radiation. We interpret the observed SED in an adaptation of parametric CPD models, including the free-free emission from the disc itself.

A purely dusty disc cannot reproduce the SED of PDS 70c, even with the inclusion of free-free radiation from a jet. However, the H I free-free continuum from the CPD yields an optically thick SED, and marginally accounts for the data (within 3 σ of Band 9). This free-free signal is intrinsic to the CPD and a product of the high temperatures reached by viscous heating in these analytic CPD models. Such an H II disc is at odds with the atomic or molecular environments expected in the context of planetary accretion.

Given an optically thick free-free spectrum, any steeper component is negligible, and the dust-to-gas mass ratio, ζ, must be very low, log(ζ) < ~4.80. Such a dust-deprived environment could be the result of strong dust filtering inside the planet-induced gap.

Uniform slab models indicate the need for lower EMs than are obtained in constant-viscosity CPDs to account for the turnover to partially thin free-free emission above Band 7 (see Fig. 11). The free-free signal would then stem from an atomic plasma, with residual ionisation from metals. The required EMs and temperatures, of T ∼ 2000 K, are naturally produced in a magnetic disc models where large disc viscosity (up to α ∼ 0.1) results from a strong magnetic field.

Magnetic disc models result in similar radio SEDs as for the H II disc, but without reaching H-ionising temperatures, and with comparatively moderate EMs (EM ∼ 10^11.5 cm⁻⁶ pc). However, because of the opacity of H⁻ free-free, an optically thick SED persists despite the lower EMs, and still overshoots the Band 9 flux density.

The CPD surface shock model is the closest one to the uniform slab requirements, and is therefore the only viable interpretation of the observed SED, including the Band 9 drop. Such a surface shock requires inefficient planetary accretion with a strong mass flux onto the CPD, and correspondingly strong equatorial mass loss.

The radio SED of PDS 70c rules out a planetary surface shock, or accretion funnels, as the source of the signal. Such small emitting regions would require temperatures in excess of 10⁶ K to reach the observed flux densities, which cannot be produced in shock models in the planetary mass regime.

Acknowledgements

We thank the referee for constructive comments. O.D., S.C., M.C. and P.W. acknowledge support from Agencia Nacional de Investigación y Desarrollo de Chile (ANID) given by FONDECYT Regular grants 1211496, ANID MAGISTER NACIONAL BECAS CHILE/2025-22250525, ANID PFCHA/DOCTORADO BECAS CHILE/2018-72190574, ANID project Data Observatory Foundation DO210001, FONDECYT postdoctoral grant 3220399 and Millennium Science Initiative Program Center Code NCN2024_001, FONDECYT grant 3220399 and ANID – Millennium Science Initiative Program – Center Code NCN2024_001. The work of OC was supported by the Czech Science Foundation (grant 21-23067M), the Charles University Research Centre program (No. UNCE/24/SCI/005), and the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2022.1.00592.S, #2022.1.01477.S, 2022.1.00893.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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

All Tables

All Figures

Fig. 1 Multi-frequency imaging of PDS 70. The black insets at the top left of each image zoom on the central region (with tick-marks separated by 0".1), while the orange insets at the bottom right zoom on the expected position of PDS 70c for circular Keplerian rotation around a 0.97 M⊙ star (with tick-marks at 0″.05). The red plus symbol in the central inset marks the nominal position of the star. In all insets the linear grey scale stretches over the range of intensities in each region, and the contour levels start at 3 σ and are incremented in units of σ. Intensities within each contour level are colour-coded differently. The two circles are centred on the positions of PDS 70b and PDS 70d (no radio counterparts are detected). Images in a) to d), along the top row, correspond to restorations of the corresponding non-parametric model images (that fit the visibility data), in e) to h), along the bottom row. A beam ellipse in a) to d) is shown in blue, on the top right. The pixel size in the model images is fixed at 4 × 4 mas2. For each band the Briggs parameter used for the restoration, the clean bean and the noise can be found in Table 2.

In the text

	Fig. 2 Point-source injections tests in Band 9. (a) Same as Fig. 1, with additional insets. (b) Same as a) but with four point sources, injected at the centre of each inset, and with a flux density corresponding to an optically thick extrapolation of the Band 7 flux density for PDS 70c.
In the text

Fig. 3 Dirty maps in natural weights of the residual Band 7 visibilities, after subtraction of the disc. The grey scale stretches over the full intensity scale, and for each of seven epochs. The inset zooms into PDS 70c, with contours at 90 and 120 μJy beam−1, fixed for ease of comparison across epochs. As in Fig. 1, the circles are centred on PDS 70b and d, and are drawn here to ease the search for any variable counterparts.

In the text

	Fig. 4 Flux density of PDS 70c at 343.5 GHz (Band 7), as a function of Julian day. The scatter is consistent with the uncertainties.
In the text

Fig. 5 Corner plot for synthetic data. The orange posteriors show the synthetic optically thin case, calculated with Ṁp = 10−9 MJup yr−1 and α = 10−3. The outcome illustrates the degeneracy between Ṁp and α inherent to the parametric CPD model of Zhu et al. (2018). The grey posteriors show the synthetic optically thick case, calculated with Ṁp = 10−5 MJupyr−1 (and the same α = 10−3). Both synthetic SEDs were obtained using only the dust component of the CPD model following Zhu et al. (2018), with errors set to 10% of the flux.

In the text

Fig. 6 Best fits for the data of Bands 3, 4, 7, and 9 using the dusty disc (Zhu et al. 2018) model and the simple dusty disc model (described in Sect. 3.2) are shown as the solid orange and light grey lines, respectively. For the dusty disc, we find that the best fit parameters of the model are $\log ({\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r}))=-{9.04}_{-2.75}^{+0.46}$ and $\log (\alpha )=-{5.21}_{-2.82}^{+0.53}$ at 1 σ. The resulting spectral index is 2.35 and has a minimum reduced χ2 of 6.50. For the simpler dusty disc model, we find that the best fit parameter is $\log (\gamma )=-{2.67}_{-0.05}^{+0.04}$. The resulting spectral index is 3 and has a minimum reduced χ2 of 14.61. Here the reduced χ2 is the χ2 divided by the number of points (data) minus the number of degrees of freedom in the model.

In the text

	Fig. 7 Best fit for the data of Bands 3, 4, 7, and 9 for the dusty disc (Zhu et al. 2018) using different Tfloor. In order, the solid orange, light grey, and grey lines correspond to Tfloor of 10, 15, and 20 K.
In the text

Fig. 8 Best fits to the data from Bands 3, 4, 7, and 9 using two free-free emission models (sum the dusty disc component from Sect. 3.2 and the free-free component from Sect. 3.3). The first has a fixed α = 10−2, shown as the solid orange line and has a minimum reduced χ2 of 3.18. We find that the best fit parameters of this model are $\log ({\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r}))=-{3.48}_{-0.09}^{+0.06}$ at 1 σ, with a strict upper limit on the dust-to-gas ratio of log(ζ) = −4.80. The second has a fixed Ṁp = 10−6 MJup yr−1 as the solid light grey line and has a minimum reduced χ2 of 2.54. We find that the best fit parameters for this model are $\log (\alpha )=-{6.86}_{-0.07}^{+0.08}$ at 1 σ, with a strict upper limit on the dust-to-gas ratio of log(ζ) = −7.05.

In the text

	Fig. 9 Physical conditions including free-free in the best-fit CPD model with α fixed. Tc is the mid-plane temperature calculated as Eq. (14), Σ(R) is the surface density, EM is the emission measure, f is the ionisation fraction in the mid-plane, and Hhwhm is the height at which the density of free electrons is half the density of free electrons in the midplane (half-width at half-maximum).
In the text

	Fig. 10 Physical conditions in the best-fit CPD model with Ṁp fixed. Annotations follow from Fig. 9.
In the text

Fig. 11 Best fit to the data from Bands 3, 4, 7, and 9 using the uniform slab model described in Sect. 3.4. On the upper right panel, the solid orange line corresponds to the case T = 103 K, has a minimum reduced χ2 of 0.18, and a p value of 83%. We find that the best fit parameters of this model are $\log (\mathrm{E}\mathrm{M}/(\mathrm{c}{\mathrm{m}}^{-6}\mathrm{p}\mathrm{c}))={11.08}_{-0.19}^{+0.53}$ and $\log \left(R_{\mathrm{m}\mathrm{a}\mathrm{x}}/\mathrm{a}\mathrm{u}\right)=-{1.11}_{-0.06}^{+0.03}$. The solid light grey line corresponds to the case T = 104 K, has a minimum reduced χ2 of 0.33, and a p value of 72%. We find that the best fit parameters of this model are $\log (\mathrm{E}\mathrm{M}/(\mathrm{c}{\mathrm{m}}^{-6}\mathrm{p}\mathrm{c}))={12.00}_{-0.18}^{+0.42}$ and $\log \left(R_{\mathrm{m}\mathrm{a}\mathrm{x}}/\mathrm{a}\mathrm{u}\right)=-{1.61}_{-0.04}^{+0.02}$. The bottom left panel shows the corner plots obtained for each model. The orange and grey posteriors also correspond to the cases of T = 103 K and T = 104 K, respectively.

In the text

Fig. 12 Results using the magnetic disc model described in Sect. 3.5. This takes into account the free-free from metals and H− and the bound-free from H− described in Sect. 3.6. On the upper right, we have the best fit for the data of Bands 3, 4, 7, and 9 represented with their error bars using the fluxes of Table 3. The modelled flux is represented by the solid line, obtaining a partially thick free-free emission with a spectral index of 1.76. At the bottom left, we have the corner plot with the best fit parameters $\log ({\dot{M}}_{\mathrm{p}}/(M_{\mathrm{J}\mathrm{u}\mathrm{p}}/\mathrm{y}\mathrm{r}))=-{5.37}_{-0.12}^{+0.06}$ and $B_{\mathrm{p}\mathrm{s}}={610}_{-375}^{+98}\mathrm{G}$ at 1 σ. The minimum reduced χ2 is 2.37.

In the text

Fig. 13 Physical conditions including metal free-free emission in the best-fit MHD model. α is the viscosity parameter, Tc is the mid-plane temperature calculated according to Eq. (14), Σ(R) is the surface density, EM is the emission measure, f is the ionisation fraction in the midplane, and Hhwhm is the height at which the density of free electrons is half the density of free electrons in the mid-plane (half-width at halfmaximum). The vertical dashed red line is the truncation radius, Rt = 1.69RJup, and the disc total mass is 1.93MJup.

In the text

	Fig. 14 Optical depth of H− (including free-free and bound-free) and metals, in the annulus that makes the greatest contribution to the flux for the magnetic disc model (Sect. 3.5). The solid blue and orange lines correspond to the optical depth by H− and metals, respectively. The dashed line is the total optical depth.
In the text

	Fig. 15 Physical conditions of the post-shock region for υ0 = 30 km s−1 and n0 = 1015cc−1. The temperature, free electron local density (ne−), and HII local density (nH+) are the solid red, black, and blue lines, respectively. The optical depths at 300 and 600 GHz (τ300GHz and τ600GHz) are the dashed green and yellow lines, where the dotted black line corresponds to τ = 1.
In the text

	Fig. 16 Physical conditions of the post-shock region for v0 = 50 km s−1 and n0 = 1012cc−1. Annotations follow from Fig. 15.
In the text