© ESO 2022

1 Introduction

The late stages in the formation of gas giant planets are characterised by the presence of a disk made of dust and gas which surrounds and feeds the forming planet. Such a structure receives the name of circumplanetary disk (CPD) (Lubow et al. 1999; Kley 1999), and it is where the formation of large and regular moons is thought to take place (Canup & Ward 2002). Hydrodynamical simulations have shown that the emergence of these structures is common both in the core accretion and in the gravitational instability models of planetary formation (Szulágyi et al. 2017). Furthermore, the architectures of large moon systems, such as those observed around Jupiter and Saturn, are consistent with the constraints imposed by satellite formation via CPDs (Sasaki et al. 2010). However, direct and conclusive observational evidence of CPDs remains elusive. The sub-millimetre emission from the surroundings of planet-forming regions is an important window to scrutinise the complexity of planetary formation. Dust particles in the CPD are heated by the visible and near-infrared (NIR) radiation received from the star and the planet. This radiation is then re-emitted at longer wavelengths, making these disks shine in the far-infrared and sub-millimetre regimes (Draine 2010; Armitage 2010).

The first report of a potential CPD candidate was done by Welch et al. (2004). Observations made at 1.4 mm of the T Tauri star HL Tau showed a localised enhancement of the surface density at a distance of 70 AU from the central star. The feature was interpreted as dust emission associated with a protoplanet candidate, but independent measurements by Carrasco-González et al. (2009) at λ = 7 mm showed that no emission above 3σ is coming from the location of the potential candidate. They suggested that the signal was probably due to free-free rather than thermal emission making the real nature of the overdensity arguable.

Kraus & Ireland (2012) studied the LkCa 15 system known for hosting a transition disk with a gap of width ~55 au. They reported the direct imaging of a potential companion candidate close to the centre of the gap which shows an extended structure when observed in L′ and K′ bands, which is not consistent with the expected footprint of an unresolved point source, such as a planet. The observed structure was interpreted as a forming gas giant planet surrounded by circumplanetary material. Very large array (VLA) observations of the 7 mm continuum emission from the same system were carried out by Isella et al. (2014), but those did not detect any emission coming from the position of the potential companion. Similarly, more recent observations carried out with the Atacama large millimeter/submillimeter array (ALMA) could not detect with enough statistical significance any localised sub-millimetre signal within the cavity of LkCa 15 (Facchini et al. 2020).

These non-confirmed CPD candidates around HL Tau and LkCa 15 illustrate how challenging it is to achieve a non-ambiguous detection of CPDs embedded in dust-depleted cavities of protoplanetary disks. Modelling and high resolution observations at multiple wavelengths are required to understand the environment around forming giant planets.

An interesting candidate for the detection of CPDs is PDS 70 and its protoplanetary disk. PDS 70 is a 0.76 M_⊙ T Tauri star (Müller et al. 2018) located at a 113.5 pc distance in the Centaurus constellation (Gaia Collaboration 2016, 2018). The first observation made in scattered light in the K_s-band was able to mask the stellar brightness, keeping an angular resolution high enough to provide evidence for the presence of a protoplanetary disk around the star (Riaud et al. 2006). Observations in bands H and L′ showed that this disk has a cavity with a radial extension of ~70 au (Hashimoto et al. 2012). Based on the observational data, and after a process of radiative transfer modelling, the authors proposed that the cavity could be the result of the interaction between the disk and a number of companions of a sub-stellar or planetary nature. Their results also led them to propose a first estimate for the upper mass limit of the companion bodies which was about 50 M_Jup. This was the first time the existence of at least one forming giant planet around PDS 70 was hypothesised.

Keppler et al. (2018) reported the discovery of a planetary mass companion inside the cavity at a distance of ~ 22 au based on near-infrared observations. In addition, this work was the first to report polarised light coming from an inner disk. Haffert et al. (2019) provided an independent confirmation of PDS 70 b via H α observations and also announced the presence of an additional Hα source inside the cavity, labelled as PDS 70 c. The latter observation was confirmed at longer wavelengths by Isella et al. (2019) who claimed that the observed sub-millimetre flux could be explained by the presence of circumplanetary material around the forming planet PDS 70 c.

The signal coming from the PDS 70 c location is spatially unresolved and its closeness to the outer protoplanetary dust ring (Keppler et al. 2019; Mesa et al. 2019) makes it difficult to isolate the observation of the planet-CPD system from the protoplanetary disk. The ability to separate the individual contributions from the planet, CPD, and protoplanetary disk to the total sub-millimetre flux is essential. Indeed, it should be noted that the CPD brightness can be greater than the planet’s, at least by one order of magnitude (Szulágyi et al. 2019). Therefore, using the combined luminosity to estimate, for example, the planetary mass can lead to an overestimate of the real value by a similar extent. The recent high resolution continuum observations of PDS 70 reported in Benisty et al. (2021) allow one to separate the CPD emission from the outer disk and placed a robust constraint on the CPD size.

To conclude, although there is robust evidence to support the presence of two gap-carving planets around PDS 70, efforts are still required to fully characterise the architecture of the system. Recent works have put important constraints on both the orbital and bulk properties of the two planets, especially for PDS 70 b which is located well inside the gap (e.g. Müller et al. 2018; Wang et al. 2020, 2021). For PDS 70 c, however, the determination of those properties is more challenging and further analysis of the circumplanetary material and its impact on the observations is required. This work aims to shed light on these questions.

This paper is organised as follows. In Sect. 2 we present a two-dimensional radiative transfer model to explain self-consistently the sub-millimetre and NIR continuum emission from the PDS 70 disk detected at 855 μm and at 1.25 μm, respectively.This reference model does not include any planetary mass companions, but paves the way to evaluate the effect that an additional source mimicking PDS 70 c has on the observations. In Sect. 3 we extend the model to a three-dimensional case by including a planetary companion at the predicted location of PDS 70 c and we show the effects this additional source would have in the synthesised ALMA image when the planet is simulated with and without a CPD around it. In Sect. 4 we derive constraints for the amount of dust in the CPD and discuss our results in light of the work of Benisty et al. (2021). We devote special attention to the thermal structure of the CPD to constrain the relative influence of the star andthe planet in the heating of the dust grains. We summarise our findings in Sect. 5.

2 Multiwavelength radiative transfer modelling

In this section we describe the process followed to find a suitable model to reproduce three observables of the PDS 70 system: the ALMA image observed at 855 μm, the VLT/SPHERE image at 1.25 μm in polarised and scattered light and the photometry of the source. To our knowledge, this is the first attempt to develop a self-consistent radiative transfer model for these particular data sets. For a description of the observation and data reduction process, we refer the reader to the works of Isella et al. (2019), Keppler et al. (2018) and references therein.

Our strategy consists in first modelling the ALMA data. This allows us to find the modulated profile for the surface density of a homogeneous distribution of dust particles. This profile is then used as an input parameter to model the Q_ϕ component of the polarised-scattered emission. We show that this step requires to abandon the idea of a homogeneous distribution of dust particles and change it for a radial segregation of particle sizes throughout the disk that also shows a differential stratification in the vertical direction. A final tuning of some geometrical parameters is also required in order to fit the modelled spectral energy distribution with the observed photometry. We used the three-dimensional version of the continuum radiative transfer code MCMax (Min et al. 2009).

2.1 Modelling the dust continuum emission

Thermal emission is dominated by the large dust particles which tend to settle down into the midplane of the disk. The radial surface density profile is a good proxy to the thermal dust emission as it traces how the particles are distributed along the disk. We started our modelling with a simple prescription for the surface density profile which is assumed axially symmetric: $${\Sigma }_{\text{dust}}(r)={\Sigma }_0{(\frac{R_{\text{tap}}}{r})}^{ϵ}\exp [-{(\frac{r}{R_{\text{tap}}})}^{2-\gamma }],$$(1)

where R_tap is the tapering-off radius indicating the distance where the transition from a linear to an exponential decay of the surface density takes place. We use ϵ = γ = 1. The constant Σ₀ is determined by the total dust mass through the relation: $$M_{\text{dust}}={\displaystyle {\int }_{R_{{\text{in}}^{\prime }}}^{R_{{\text{out}}^{\prime }}}{\displaystyle {\int }_0^{2\pi }r}}{\Sigma }_{\text{dust}}(r)\text{d}\varphi \text{d}r,$$(2)

with $R_{{\text{in}}^{\prime }}$ and $R_{{\text{out}}^{\prime }}$ being the inner and outer limits of the protoplanetary disk, respectively. We adopt $R_{{\text{in}}^{\prime }}=0.04\text{au}$ following Keppler et al. (2018). For the outer radius we assume $R_{{\text{out}}^{\prime }}=130\text{au}$ since the value of the azimuthally averaged flux at 855 μm oscillates around zero Jy for r > 120 au. A typical value for R_tap of 100 au is taken from the T Tauri reference model introduced by Woitke et al. (2016). The initial surface density profile is shown in Fig. 1. The disk has a flaring exponent β and a scale height H(r) with a characteristic value H₀ at a reference distance r₀ = R_tap: $$H(r)=H_0{(\frac{r}{r_0})}^{\beta }.$$(3)

We assume an initial distribution of dust particles with sizes between 0.05 μm and 3 mm with a size distribution n(a) such that the number of particles with sizes between a and a + da follows a power law: $$n(a)\text{d}a\propto a^{-p}\text{d}a,$$(4)

with p being a positive exponent. The dust grains are a mixture of astronomical silicates and amorphous carbon with a porosity accounting for certain fraction of the total volume of the particle. This part of the modelling procedure used the DIANA standard dust opacities¹ (Woitke et al. 2016).

In MCMax3D, the treatment of dust settling follows the prescription given by Dubrulle et al. (1995). The settling is parameterised through the turbulence parameter α (Shakura & Sunyaev 1973; Min et al. 2012); small values of α imply a high level of settling. For a fixed value of α, the level of settling is higher for large and compact dust grains, although smaller particles can also follow the same trend but with a different efficiency. Therefore, the settling in our model is consistently implemented using a single parameter that affects both the small and the large grains. In Sect. 4.2, we discuss in more detail the parametrisation of the vertical dust distribution.

For the vertical density structure, we use the following prescription: $$\rho (r,z)=\frac{\Sigma (r)}{\sqrt{2\pi }\cdot H(r)}\exp (-\frac{z^2}{2H{(r)}^2}),$$(5)

where Σ(r) = 100 × Σ_dust(r) is the gas surface density and H(r) is given by Eq. (3).

We first ran a Monte Carlo simulation with the initial density profile shown in Fig. 1. The obtained disk structure was ray-traced at λ = 855 μm and the resulting image was convolved with a Gaussian kernel that represents the ALMA beam of the observed data set. The beam size is 67 × 50 mas (~ 8 × 6 au) and has a position angle of 61.5° east of north. Then, we performed aperture photometry both in the observed and synthesised images to extract and compare the respective azimuthally averaged radial brightness profiles. In order to do this, we used a series of elliptical annular apertures with constant width (equal to 1∕4 of the beam’s major axis) and increasing values of the inner semi-major axis to cover the whole extension of the disk, from $R_{{\text{in}}^{\prime }}$ to $R_{{\text{out}}^{\prime }}$. The eccentricity of the annulus is given by the sine of the disk inclination, measured with respect to the plane of the sky 51.7°.

To fit the main characteristics of the ALMA observation, the initial surface density profile had to be modified. This modification was done through the following iterative procedure which has been used in previous works (e.g. Pinte et al. 2016; Muro-Arena et al. 2018; Rab et al. 2020). The density profile was multiplied point-by-point by the ratio of the observed to modelled azimuthally averaged radial brightness profiles at each radial point. With the modified profile, we repeat the same procedure: run the Monte Carlo simulation, ray-trace the system and convolve the resulting image, extract a new azimuthally averaged radial brightness profile and compare with the observed one. After five iterations, our synthesised radial profile converged to a good extent towards the observed one and the iterations were stopped. The iteratively modified surface density is shown in Fig. 1 together with the profile reported by Keppler et al. (2018) for comparison.

The top panels in Fig. 2 show a comparison between the observed image and that obtained with our model. Using Eq. (2) we derive a dust mass of 4.4 × 10⁻⁵ M_⊙ in the protoplanetary disk, which is slightly higher than the value found by Keppler et al. (2018) via modelling of the near-infrared emission

The ALMA observation shows an asymmetry in the azimuthal distribution of the brightness in the outer disk which is especially visible along the major axis. The existence of this feature was pointed out by Long et al. (2018), confirmed by Keppler et al. (2019) and recently by Benisty et al. (2021) with observations at high angular resolution. To quantify this asymmetry, we measured the flux density as a function of the position angle inside a 2 au wide annulus centred at 75 au. We found the ratio of the maximum to minimum values to be 1.5 in the observation but only 1.2 in the simulated image. Our model does not aim to reproduce the observed asymmetry in the outer disk and any asymmetry is smeared out during the process we implemented to find the modified density profile. However, given the quality of the fit that our model provides to the photometric points at (sub-)millimetre wavelengths, and our focus to investigate the CPD, we did not modify the density profile any further and leave the detailed analysis of the asymmetry as a goal for a future study. It is also important to note from Fig. 1 that we restricted the column density to a minimum of 10⁻⁵ g cm⁻². This value is in essence determined by the sensitivity of the ALMA observation.

Finally, our approach to simulate the ALMA observation by convolving the ray-traced image with a Gaussian beam misses complex details that might arise in the reduction of interferometric data, such as spatial filtering and the presence of correlated noise in the final image. According to the array configuration used to image the source, the maximum recoverable scale is of the order of 1000 au (Isella et al. 2019). Therefore, this does not affect the analysis of compact emission sources, such as the inner disk or the CPD, given the small angular scales covered by them. On the other hand, we show in Sect. 4.3 that even the fingerprint of the lowest mass CPD considered in this study can be marginally detected above the rms noise level of the observation. The inclusion of these effects, although important for a detailed comparison with the observations, will not modify the fundamental conclusions of our work given the angular resolution and sensitivity of the analysed data set.

Fig. 1 Surface density profiles before and after the iterative procedure described in Sect. 2.1. The modified profile provides the best fit to the ALMA data. The grey area depicts the spatially unresolved inner disk due to the ALMA beam size.

Fig. 2 Comparison between the ray-traced images retrieved from our model (right column) and the ALMA-VLT/SPHERE observations (left column). Top panels: thermal emission at λ = 855 μm. Bottom panels: Qϕ component of the polarised and scattered emission at λ = 1.25 μm. North is up,and east is to the left.

2.2 Modelling the polarised and scattered light emission

Using the results obtained after modelling the ALMA image, we performed an additional ray-tracing at 1.25 μm leading to the simulated Q and U Stokes frames for linear polarisation. As the polarised intensity, given by $\text{PI}=\sqrt{Q^2+U^2}$, tends to be a small fraction of the total intensity (~10%), the absolutevalue of the independent Stokes frames tends to be close to the respective noise level. As a consequence, measurements based on the polarised intensity commonly come with large errors of systematic nature (Schmid et al. 2006). To circumvent this problem, we worked with the polar-coordinates Stokes parameters Q_ϕ and U_ϕ defined as: $$\begin{array}{ll}{Q}_{\varphi }\hfill & =+Q\cos (2\varphi )+U\sin (2\varphi ),\hfill \\ U_{\varphi }\hfill & =-Q\sin (2\varphi )+U\cos (2\varphi );\hfill \end{array}$$

where $$\varphi =\arctan (\frac{x-x_0}{y-y_0})+{\varphi }_0,$$(8)

is the angle between the position vector of a pixel with coordinates (x, y) in the Q and U frames, and a unit vector pointing to the north direction, with the origin of both vectors placed at the position of the star (x₀, y₀). The term ϕ₀ accounts for extra polarisation effects such as those induced by the instruments and was set equal to zero in our models. In this parametrisation, a positive signal in Q_ϕ accounts for scattered light polarised in the azimuthal direction whereas Q_ϕ negative implies radial polarisation. The parameter U_ϕ contains the light polarised with a direction ± 45° with respect to the radially polarised light.

Unlike the disk observed at ALMA wavelengths, the Q_ϕ frame derived from the observations is clearly non-azimuthally symmetric (see left figure in lower panel of Fig. 2). Instead of trying to fit an azimuthally averaged profile, we aimed to reproduce the Q_ϕ signal within a conical aperture centred along the major axis which is ± 10° wide in radial steps of ~3 au. This corresponds roughly to one-half of the resolution of the IRDIS polarimetric J-band observation (non-coronagraphic case) presented by Keppler et al. (2018). We did not consider stellar distances smaller than the image resolution of 49 mas (~5.6 au). The uncertainties were computed as the standard deviation of the signal within the correspondent radial bin. The resulting profile retrieved from the observation is shown in black markers in the central panel of Fig. 3.

A first comparison with the modelled profile failed at capturing the main features of the Q_ϕ frame which suggested to modify other parameters of our model. Scattered light is a tracer for the surface layers of protoplanetary disks which are mainly populated by small dust grains whose physical properties determine the polarisation state of the scattered wave. This suggests that a better representation of the data set could be obtained by modifying the physical properties of the dust grains and the geometry of the disk. We motivate the modification of the dust properties upon the assumption that some mechanisms of dust growth and radial drift are taking place in PDS 70.

With MCMax3D we can divide the protoplanetary disk in different zones. Each zone can have a distinct geometry respect to the others, and two zones can overlap to mimic, for example, the presence of an embedded CPD. We divided the protoplanetary disk in two zones: the inner from 0.04 to 40 au and the outer from 40 to 130 au. With these limits,the dust cavity and the inner disk are within the inner zone whereas the main dust ring is in the outer one. Our numerical experiments suggest that this configuration can capture many of the features observed in scattered light. Each zone has its own dustsize distribution (see Table 1). Our underlying assumption for the dust segregation is motivated by the results of Bae et al. (2019), where two-dimensional hydrodynamical simulations of the PDS 70 system, led them to demonstratethat large grains are filtrated at the edge of the gap whereas the inner protoplanetary disk is only populated by micron-sized grains.

With this information at hand, we populated the inner zone with grains of maximum size 100 μm whereas the outer zone contains larger grainswith a maximum size of 3 mm. Other bulk properties such as the power law index of the size distribution, the amount of carbon content and the porosityhave all equal values in the two zones but slightly different to those used in Sect. 2.1 to fit the ALMA observation. The modification of these properties was also needed to improve the fit of the SED (see Fig. 3 right panel). The turbulence level, controlled by the settling parameter α, decreases from the inner to the outer zone. The implemented values are all within the ranges explored by diverse authors for other systems (e.g. Muro-Arena et al. 2018). Table 1 summarises the whole set of properties of our radiative transfer model. In Fig. A.2 we show the extinction, absorption and scattering opacities for each zone.

Our fit to the Q_ϕ observation can still be improved. This might impose stronger constraints on the population of small size grains and their properties. We note that we cannot reproduce the azimuthal asymmetry in the brightness distribution of the outer disk, where the NE side appears brighter than the NW side (see Fig. 2 in this work and also Fig. 3 in Keppler et al. 2018). This effect might be connected to the scattering properties of the particles populating the surface of the outer disk, since the degree of forward scattering increases when those layers consist of larger grains. Geometrically, the NW side of the disk also represents the near side of the disk and therefore, an enhancement in the number of forward scattered photons due to the larger grains would make this side to appear brighter. However, since the main goal of the subsequent sections is related to the sub-millimetre regime, we do not aim to pursue a better fit to the infrared data.

A couple of differences in the assumptions made by Keppler et al. (2018) and this work stand out. Firstly, since our model assumes an inherent inhomogeneity in the radial distribution of grain sizes, our inner and outer disks have different opacities, as can be seen in Fig. A.2. At visual wavelengths, the opacity of the inner disk is is a factor of a few larger compared to the outer disk opacity, which is explained by the presence of smaller grains populating the inner disk. This impacts the modelling of structures in the outer disk since we have to compensate for the fewer photons reaching the outer disk by, for example, increasing the outer disk mass. In contrast, the dust opacities in Keppler’s model are expected to be the constant throughout the disk and the absolute values are also expected to differ from ours because of the differences in grain composition and size distribution slope². Secondly, the settling of dust in our model is controlled by the Shakura-Sunyaev α parameter through the steady-state solution of the diffusion equation for the dust density (see Sect. 4.1). This approach differs from Keppler et al. (2018) where the disk was modelled as a two-layer structure with different scale heights to account for the difference in the vertical distribution of large and small grains. Although both approaches provide a good fit to the data, we note that in order to fine-tune the value of the alpha parameter, a close look to the SED is required since this parameter determines the behaviour of the flux in the multicolour and sub-millimetre part of the spectrum. An SED fitting leaves a degeneracy between several parameters, namely the dust mass in the outer disk, the scale height, and the turbulence parameter. That said, we base our subsequent analysis on the parameters listed in Table 1 but we point out that these will not be unique

Fig. 3 Left panel: azimuthally averaged radial brightness profile from the ALMA image (black dots) and from our model (blue line). Middle panel: radial cut along the major axis of the Qϕ frames (negative radial direction is along the red-shifted semi-major axis, which is nearly coincident with the south-east direction). Right panel: sepectral energy distribution towards PDS 70; black dots correspond to the source’s photometry and the blue line is our prediction.

3 The effect of a planetary mass companion

3.1 The effect of a planetary companion on the sub-millimetre flux

We explore the fingerprint of a planetary companion assuming no circumplanetary material around it. The implementation of an additional photon source requires defining its luminosity and effective temperature. Keppler et al. (2019) pointed out the existence of a bridge feature located at a position angle of 288° emerging from the outer ring and possibly connecting it with the inner disk. The closeness of this spur to the location of PDS 70 c alongside the inverted P-cygni profile observed in Hα, may indicate that the planet is still accreting material from the protoplanetary disk. However, the accretion rate ~ 10⁻² M_Jup Myr⁻¹ (Haffert et al. 2019) makes it unlikely that the planet is going through a runaway phase of growth.

The model for giant planet formation presented in Ginzburg & Chiang (2019) is applicable to a post-runaway scenario. Ginzburg & Chiang (2019) have shown that just after the runaway phase, during the runaway cooling, a planet with a two-layered atmosphere (a convective layer surrounded by a radiative one) will contract and increase its mass with nearly similar characteristic times. As both processes are regulated by the cooling efficiency and by the accretion rate respectively,the similarity between these characteristic timescales allows us to formulate the following equivalence: $$4\pi \sigma R_{\text{p}}^2{T}_{\text{eff}}^4=\frac{G{M}_{\text{p}}{\dot{M}}_{\text{p}}}{R_{\text{p}}},$$(9)

where M_p, Ṁ_p, R_p and T_eff represent the mass, mass accretion rate, radius and effective temperature of the planet; σ and G are the Stefan-Boltzmann and gravitational constants, respectively. As demonstrated by Ginzburg & Chiang (2019), the common assumption of taking R_p ≈ 2 R_Jup during the last Myrs of the planet formation (e.g. Eisner 2015) is well justified as long as the opacity at the boundary of the atmospheric layers adopts low values of the order of ~ 10⁻² cm² g⁻¹. We assumed this and worked with a fixed planetary radius of R_p = 2 R_Jup. Therefore, an estimate for the effective temperature can be derived using Eq. (9) once a value for the planet mass is assumed.

Several works have constrained the mass of PDS 70 c, finding values ranging from ~ 0.5 to 12 M_Jup (e.g. Haffert et al. 2019; Mesa et al. 2019; Bae et al. 2019; Wang et al. 2020) (see Table 2). We explored a low and a high mass regime given by these values and a third intermediate case with a mass of 5 M_Jup. The spectral energy distribution of the planet is modelled as a blackbody. The location of the planet is taken to match the position angle and the mean orbital radius of the dynamically stable orbital solution for PDS 70 c reported in Wang et al. (2020) and in Wang et al. (2021), namely, r_c = 34.3 au and PA_c = 280.4°.

We ran one radiative transfer calculation for each value of the companion’s mass. We used 5 × 10⁸ photon packages and the resulting disk structure was ray-traced at 855 μm. Analysing the temperature structure of the protoplanetary disk, we find that the planet does heat the dust material in its close vicinity, even though the surface density is just 10⁻⁵ g cm⁻² (see Fig. 1). Indeed, a locally enhanced temperature of ~40 K is measured inside the computational cell where the planet is located. This is notoriously different to a temperature of ~ 20 K that is measured at the diametrically opposed computational cell. These values refer to the midplane temperatureand correspond to the simulation that includes a planetary companion with an effective temperature of 1054 K (Wang et al. 2021). The ray-traced images were convolved with the same beam properties described in Sect. 2.1. After the convolution, we find that the signal from the planet is smeared out to an extent that it is not possible to differentiate it from the simulation without the planet. Therefore, our simulations suggest that although the luminosity of PDS 70 c can modify the local temperature of the circumstellar material, the observational evidence of such an effect could not be retrieved with ALMA in the continuum. With this result at hand, we proceed to assess if the sub-millimetre signal from the vicinity of the planet can be explained by the presence of circumplanetary material.

3.2 The effect of a circumplanetary disk around PDS 70 c on the sub-millimetre flux

In light of the recent work by Benisty et al. (2021) that constrained the CPD size and the works of Wang et al. (2020, 2021) that provideestimates for the temperature and luminosity of PDS 70 c, we investigate the observational fingerprint of a CPD using two different approaches. First, we fix the planet luminosity and temperature and the CPD size to their recent estimates and vary the amount of dust in the CPD to reproduce the observation. In the second approach, we fix the CPD mass found in the previous step and investigate the impact of having a less luminous planet. We defer this second part to the discussion in Sect. 4.3.

3.2.1 Constraining the dust content in the CPD

We modelled the CPD as a flared disk centred at the location of the planet. The outer radius of the CPD was set to 1.2 au according to the upper limit set by Benisty et al. (2021). Assuming this value is equal to the truncation radius of the CPD, which has been shown to scale as one-third of the planet Hill’s radius³ (e.g. Quillen & Trilling 1998) we derive a rough estimate for the planet mass of 2.8 M_Jup. We set the inner radius of the CPD to 0.007 au⁴ and we populated it with dust grains ranging from 0.05 μm to 3 mm in size with the same content of amorphous carbon and porosity as the grains in the protoplanetary disk. Opacities are again computed using the DIANA standard (see Fig. A.2). The rest of the physical and geometrical properties for the CPD are listed in Table 1. A simulation with 5 × 10⁸ photon packages was ran for three different CPD masses, namely, 0.007 M_⊕, 0.07 M_⊕ and 0.7 M_⊕. The results of these simulations are shown in Fig. 4 and discussed in detail in Sect. 4.3. Unlike the simulations presented in Sect. 3.1 where only the planet was included, it is now evident how the CPD modifies the ray-traced image. The photons emitted by the planet are absorbed and reprocessed by the dust grains in the CPD and their subsequent re-emission boosts the sub-millimetre flux.

Numerical simulations show that the depletion timescales for the dust component in CPDs are smaller than the respective counterpart in protoplanetary disks. For example, Rab et al. (2019) show, using dust evolution models applied to CPDs around wide-orbit companions, that the dust-to-gas ratio can decrease from 10⁻² to 10⁻⁴ in just 1 Myr. This fast depletion timescale can be mitigated by invoking the presence of dust traps such as those proposed byDrążkowska & Szulágyi (2018). They demonstrated that the advection of gas directed outwards can create an efficient dust trap allowing dust grains to accumulate and grow. In order to constrain the gas-to-dust ratio in CPDs further observations at high angular resolution both in the continuum and in emission lines are required. Therefore, in the absence of observational constraints, we work with a canonical value of 100 for the gas-to-dust ratio throughout the CPD.

To try and interpret the result shown in Fig. 4 a bit further, at least for the low mass regime where some regions of the CPD are expected to be optically thin, we note that if the CPD is located at a distance d from an observer on Earth with an inclination i with respect to the plane of the sky, thenthe power at frequency ν received per unit area is (Woitke 2015): $$\nu F_{\nu }=\frac{2\pi \cos i}{d^2}{\displaystyle \underset{r_{\text{in}}}{\overset{r_{\text{out}}}{\int }}r}\nu B_{\nu }(T(r))[1-\exp (-\frac{{\tau }_{\nu }(r)}{\cos i})]\text{d}r,$$(10)

where the integration limits are determined by the CPD size, T is the dust temperature at a planetocentric position r, and τ_ν is the position-dependent vertical optical depth: ${\tau }_{\nu }(r)={\kappa }_{\nu }^{\text{abs}}\Sigma (r)$. Using the Rayleigh-Jeans formula, Eq. (10) becomes⁵: $$\nu F_{\nu }\approx \frac{2{\nu }^3{k}_{\text{B}}}{{(dc)}^2}{\kappa }_{\nu }^{\text{abs}}\langle T\rangle M_{\text{dust}},$$(11)

with c, the speed of light. It is clear then that the flux from an optically thin CPD scales linearly with its dust content. However, due to the relatively high values of the surface density, especially in the inner regions of the CPD (see Sect. 4.5), we expect optically thick emission from the CPD. In fact, analysing the density profile of the disk with a value for the absorption opacity of 8 cm² g⁻¹ (see Fig. A.2), we found that only the 0.007 M_⊕ disk has an optical depth lower than 1 at 855 μm for r > 0.01 au. The two more massive disks are entirely optically thick. A full numerical radiative transfer treatment for the study of these compact structures seems thus justified.

4 Results and discussion

4.1 Particle size segregation

Our modelling requires to segregate dust grains by size since the initial implementation of a radially homogeneous distribution of small and large grains did not provide a good fit to the scattered light observation and multiple photometric points were too far from the predicted SED. We do not attempt to have a consistent dust evolutionary model aimed to explain the dust grain segregation but we note that such behaviour has been implemented also in other works to explain the observations of protoplanetary disks (e.g. Muro-Arena et al. 2018; Villenave et al. 2019).

In general, small dust particles tend to be well coupled to the gas component which is moving at sub-Keplerian velocities as a result of radial pressure gradients in the disk. Although the motion of larger grains is not directly affected by the pressure gradient, the interaction with smaller grains and gas particles will make the former lose angular momentum and migrate inwards (Birnstiel et al. 2010). However, this migration can be halted via several mechanisms. Some of those invoke the presence of dead zones or giant planets. Dead zones are weakly ionised regions in a protoplanetary disk characterised by the suppression of gas flow (Pinilla et al. 2016). The accumulation of gas at the outer rim of the dead zone leads to an increment in the pressure facilitating the trapping of large grains (Villenave et al. 2019). On the other hand, a giant planet with a mass larger than 1 M_Jup can open a gap in the gas distribution leading again to the enhancement of the pressure at the outer border of the gap. In particular, the scenario involving a giant planet may apply to PDS 70 where at least two planets have been detected. In the case of a two-planet system, we note that Pinilla et al. (2015) carried out hydrodynamical and dust evolution simulations and find that, depending on the disk viscosity and fragmentation velocity of the grains, the presence of a second planet can modulate another overdensity region independent of the one created by the first planet. Regardless the operating scenario, only the smaller particles can penetrate the pressure wall and migrate inwards while the larger ones remain trapped (Fouchet et al. 2007).

Fig. 4 Synthetic images including a planetary companion surrounded by a CPD. The planet was placed at a distance of 34.3 au with a position angle of 280.4° measured east from north. The three images have similar pixel scales (4 mas px−1) and were convolved with the same beam properties of the ALMA observation (see Sect. 2.1). The colour scale is set equal to the square root of the surface density to boost the signal from the CPD. North is up, east is to the left.

4.2 Vertical distribution of dust grains

In our models, the vertical stratification of grains is implemented following the prescription of Dubrulle et al. (1995) where the settling is the result of a diffusion process mediated by the gas turbulence and the gravitational effect of the central body. The strength of the turbulence is proportional to the turbulent viscosity and can be parameterised by the dimensionless α parameter introduced by Shakura & Sunyaev (1973). Prior to the formal Monte Carlo iteration, the steady-state solution of the diffusion equation is computed to find the dust density as a function of the vertical coordinate for each bin of the size distribution.

The steady-state solution to the diffusion equation can be written in terms of the scale height H_d (r, a) of a dust grain with size a. This expression is simply (Dubrulle et al. 1995; Woitke et al. 2016): $${(\frac{H_{\text{d}}(r,a)}{H(r)})}^2=\frac{\alpha {(1+{\gamma }_0)}^{-1/2}}{a{\rho }_{\text{gr}}}{\rho }_{\text{mid}}(r)H(r),$$(12)

where, γ₀ ≈ 2 for compressible turbulence, ρ_gr is the mass density of a dust grain (ρ_gr = 2.18 g cm⁻³ in our simulations) and ρ_mid is the gas density in the midplane (see Eq. (5)). Therefore, for each cell in the disk we identify a particle size distribution of the form $f(a,r)\propto \text{exp}(-z^2)/[2H_{\text{d}}{(r,a)}^2]$, where the vector r denotes theposition of the cell in the three-dimensional grid and z = rcosθ is the height above the midplane for a point with co-latitude θ.

To visualise the effect of the vertical segregation we first define the mean grain size via the third moment of the size distribution: $$\overline{a}(r)\equiv {\langle a^3\rangle }^{1/3}(r)={[\frac{1}{\overline{n}}{\displaystyle {\int }_{a_{\text{min}}}^{a_{\text{max}}}f}(a,r)a^3\text{d}a]}^{1/3},$$(13)

where $\overline{n}$ is the total number density of representative equal-sized particles as used in Vasyunin et al. (2011). The resulting structure is computed for both zones in the protoplanetary disk and shown in Fig. 5. The shape of the contours is clearly related to the surface density profile which in turn modulates the shape of the gas density.

Inside the gap, the level of settling is enhanced since the surface density drops to 10⁻⁵ g cm⁻² (see Fig. 1); the larger particles in the inner zone are settled to the midplane whereas for z values close to one scale height, the smaller particles prevail. On the other hand, the right panel of Fig. 5 shows how particles in the outer zone are vertically segregated. The main dust ring shows a low degree of settling, even though the value of α is slightly lower than in the inner zone (see Table 1). This is again related to the high value of the gas density around 80 au. In Fig. 5, it is observed an apparent lack of grains with size close to the a_max value of the respective zone. This can be understood considering the following arguments. First, this has to do with the relatively high value of the exponent of the size distribution used in the model (see Eq. (4)). Such a high exponent increases the steepness of the distribution function which in turn amplifies the statistical weight of the smaller grains. Second, it is also a numerical effect due to our definition of the mean particle size which is taken to be proportional to the cubic root of the third moment of the distribution. We note that other definitions are also suitable and they can make explicit the presence of larger particles in Fig. 5 (e.g. Vasyunin et al. 2011; Facchini et al. 2017). However, the conclusions about the vertical stratification will remain unaltered regardless the formalism used to define the mean particle size.

4.3 Constraining the dust environment around PDS 70 c

A visual inspection of the ALMA data set (Fig. 2) suggests the presence of a dim spot nearly matching the location of the dynamically stable orbital solution for PDS 70 c in Wang et al. (2021). Such a feature can be retrieved from our simulations when we include a CPD that is massive enough, as shown in Fig. 4. Therefore, we use our models to measure the brightness profile in the vicinity of PDS 70 c and aim to draw some conclusions about the expected mass ranges of the CPD. We extracted the radial brightness profile using a conical aperture with ± 10° width around the position angle of PDS 70 c (PA = 280.4° measured east-from-north) with radial bins equally spaced by 0.017 arcsec (~ 2 au), equivalent to one quarter of the ALMA beam’s major axis. This procedure was performed with the GoFish package (Teague 2019). Figure 6 shows the cuts obtained from the observation compared to the three CPD mass cases simulated in this work.

Although the width and height of the modelled curves fit the observation, all the curves show an offset of ~ 2 au with respectto the observation. It might be that our fitting procedure based on extracting azimuthal averages is not well suited to reproduce the observed flux along individual directions especially when they present evident asymmetries. We do not attempt to dive into this issue. Despite the mismatch, our results can be used to explore the dust content in the CPD. The blue points in the interval 26 ≲ r ≲ 40 au were extracted from the side of the disk opposite to the direction of PDS 70 c (i.e. PA ≈ 100°) where there is no evidence of emission other than the background continuum. Due to the level of dust depletion in the cavity, this background emission and its uncertainty allow us to conclude that a 0.007 M_⊕ CPD would be marginally detectable with the sensitivity of our data set. Figure 6 also depicts the 19 μJy beam⁻¹ rms noise level of the observation as reported by Isella et al. (2019). The photometric curves in Fig. 6 indicate that the model with a 0.7 M_⊕ CPD predicts a flux density that agrees with the observed values within the error bars.

Our estimate for the mass of PDS 70 c of ~2.8 M_Jup is very closeto the upper value estimated by Wang et al. (2020) via SED fitting (see Table 2). This is remarkable considering that we used a very simple argument, namely, that the CPD radius equals one-third of the Hill’s radius. However, as noted by Wang et al. (2020), the extent of their wavelength coverage opens the question whether or not the luminosity they derived properly approximates to the bolometric luminosity needed to estimate the mass of the planet with the evolutionary model of Ginzburg & Chiang (2019). Part of the bolometric luminosity could be captured by the CPD and then radiated away at wavelengths longer than ~ 4 μm, which is the limit of the wavelength coverage of Wang’s data set. If this is the case, an infrared SED fitting with a limited wavelength coverage will underestimate the total luminosity because the SED is expected to peak at longer wavelengths (Zhu 2015; Szulágyi et al. 2019). Consequently, the planet mass will be underestimated. In our models, the assumed mass of the planet affects the dust settling within the CPD and this could alter our predictions for the sub-millimetre flux. Additional observationsat longer wavelengths are needed to obtain more robust estimates to the mass of PDS 70 c via SED fitting.

The amplitude of the crest at 34 au in Fig. 6 depends not only on the CPD mass but also on the planet properties, particularly on its luminosity and mass. In order to test how sensitive our results are to these variables, we abandon the estimates provided in Wang et al. (2021) and use instead the theoretically expected values from the Ginzburg & Chiang (2019) model for the low, intermediate and high planet mass cases (see Table 2). In all cases we fixed the size of the CPD to 1.2 au in accordance to the most recent observations (Benisty et al. 2021). The mass of the CPD was fixed to 0.7 M_⊕ and the results are shown in Fig. 7. Evidently, the assumed properties for the planet impact the flux emitted by the optically thick CPD. Since the dust mass was kept fixed and even the highest planetary luminosity is still lower than Wang’s luminosity, it is not a surprise that the observations tend to be better explained by high mass planets. Given the results presented in this subsection, we claim that it is not possible to constrain individually the mass of the planet and the dust budget in the CPD.

We performed aperture photometry in the ray-traced images using a circular aperture 9 au in radius around the planet’s mean location. When the CPD mass is 0.07 M_⊕, the measured flux is 64% of the 0.7 M_⊕ case and when the CPD mass is 0.007 M_⊕ the flux represents just ~37%. This indicates that at least one of the assumptions made going from Eqs. (10) to (11) (which predicts a linear scaling) breaks and we reinforce the idea that the CPD cannot be treated as optically thin.

We can also compare the values of the peak intensities between our model and the ALMA observation. Isella et al. (2019) reported a CPD flux peak value of 106 ± 19 μJy beam⁻¹ whereas our model predicts: 150, 106 and 58 μJy beam⁻¹ for the 0.7, 0.07 and 0.007 M_⊕ CPD mass cases, respectively. Similar to the brightness profiles in Fig. 6, this result also suggests that the mass of the CPD is constrained to the range 0.07−0.7 M_⊕.

Fig. 5 Left: cut in the vertical distribution of mean particle sizes in the inner zone of the protoplanetary disk. Labels along the dashed contours show the mean particle size in μm. Right: same as the firgure on the left, but for the outer zone.

Fig. 6 Radial cuts extracted from Fig. 4. The conical aperture was centred at a position angle of 280.4° in line with the astrometric position of PDS 70 c. The aperture has a width of ± 10° around the mean position and the radial step size is 2 au.

Fig. 7 Effect on the CPD flux due to the adopted planet mass and luminosity. In this experiment, our reference CPD mass was fixed to 0.7 M⊕. Vertical dashed lines represents the limits of the circular aperture used to measure the flux emitted by the CPD.

4.4 Comparison with recent works

Benisty et al. (2021) observed the continuum emission at λ = 855 μm with an angular resolution of ~20 mas. The observation clearly shows a localised signal in the expected location of PDS 70 c and it is detached from the disk although still unresolved. The authors interpret this sub-millimetre emission as originating from dusty material in a CPD.

They detected emission from a resolved inner disk and argue that the presence of the planets can produce an effective mechanism for dust trapping, so only the small dust grains can flow towards the inner disk. This agrees with our findings, where the particles in the inner disk have sizes from 0.05 to 100 μm. Our model for the protoplanetary disk predicts a flux from the inner disk equal to 0.76 mJy, which is very close to the observed mean value of 0.85 mJy. We estimated the dust mass in the inner disk to be 0.02 M_⊕. This value is low in comparison to the 0.08 M_⊕ lower limit estimated in Benisty et al. (2021). This difference could be explained by the steepness of our size distribution (p = −3.9). This causes the abundance of large grains in the inner zone to be lower than in Benisty’s work (p = −3.5). Combining the result in Fig. 1 with an absorption opacity of 8 cm² g⁻¹ suggests that the inner disk is optically thin at λ = 855 μm.

Working under the assumption of an optically thin CPD populated with submicron and centimetre-sized grains, Benisty et al. (2021) estimated the dust mass in the CPD as a function of the maximum grain size. In particular, when the dust size exponent is assumed to be − 3.5, they found a mass of ~0.009 M_⊕. Our 0.007 M_⊕ CPD case does not provide a satisfactory fit to the observation. As our more massive CPDs become optically thick, we only comment on the comparison of these two low-mass cases. First of all, the semi-analytical estimate of Benisty et al. (2021) only depends on the absorption opacity. However, at mm wavelengths and assuming a size distribution with maximum size above 100 μm, the scattering opacity exceeds the absorption opacity (Birnstiel et al. 2018). A similar behaviour is observed in the DIANA opacities (see Fig. A.2 for the CPD zone). Therefore, the inclusion of scattering effects is expected to alter the emitted flux even in optically thin systems (Zhu et al. 2019). Our model takes into account those non-isotropically scattered photons that will not contribute to the dust heating. As a result, the same amount of dust will emit a lower millimetre flux when scattering is considered compared to the case when scattering is neglected.

Another point relates to the assumed CPD temperature. Benisty et al. (2021) estimated the dust temperaturefrom the combined effect of several heating sources: the viscous accretion and the planet and stellar irradiation at the location of PDS 70 c. With this assumptions, they found a temperature of ~ 26 K at 1 au. In contrast, our simulation predicts a temperature of 10 K at the same distance. Since the dust mass is inversely proportional to the temperature in an optically thin disk (see Eq. (11)), our lower estimate for the temperature requires a higher dust mass to reproduce the same flux.

Finally,we note that although most of the emission from the 0.007 M_⊕ CPD is expected to be optically thin, our model suggests that the innermost region could still be optically thick because of the dust density decaying exponentially with distance. This, of course, is just an assumption since spatially resolved observations of CPDs are still not available. However, it raises the question whether or not the heating of those grains in the outskirts is due to photons from the planet’s photosphere travelling freely through the medium. Infrared photons re-emitted by the optically thick layers might have an important role setting the equilibrium temperature of the outer CPD. We note that the potential effects of optically thick emission are accounted in our model.

4.5 The thermal environment of the CPD

Since the CPD was modelled as an additional zone placed on top of the protoplanetary disk, the thermal structure was self-consistently calculated during the same Monte-Carlo run that allowed us to find the temperature in the protoplanetary disk. In our simulations, the dust temperature in the CPD is set by the stellar photons reprocessed by the protoplanetary disk and by the planet luminosity⁶, after a prescription for the vertical structure of the CPD has been provided (see Table 1). We aim to disentangle the contribution to the total CPD temperature of both sources by switching on and off the effect of the planet luminosity in the simulation.

We use the reference model for the protoplanetary disk described in Sect. 2 and include a 2.8 M_Jup planet with the same estimates for the luminosity and effective temperature as given by Wang et al. (2021). We also include an optically thin CPD with a dust mass of 0.007 M_⊕ and a size of 1.2 au. Considering an absorption opacity of 8.5 cm² g⁻¹ and the density profile of Fig. 8, the CPD will be optically thin at sub-millimetre wavelengths down to 0.01 au; this ensures that a high enough number of photon packages is able to reach the disk’s midplane and a reliable dust temperature estimate can be then retrieved.

The whole PDS 70 system was simulated with and without the planet and the azimuthally averaged midplane temperaturewas extracted afterwards for both cases. Results are shown with red points in Fig. 9. The filled circles represent the temperature obtained when the planet was included and the open circles when the planet was removed from the simulation. To quantify the relative contribution of the star to the CPD heating, we computed the ratio between the simulations and identified the distance where it becomes larger than 0.5. Beyond this location, the heating due to stellar irradiation (direct or scattered by the protoplanetary disk) contributes more than 50% to the total CPD heating. This location is marked with the vertical red line. The blue curve in Fig. 9 is the temperature of an isolated CPD with the same properties but not embedded within a protoplanetary disk so the only heating source is the planet. The horizontal line in Fig. 9 at 19 K represents the predicted temperature of the protoplanetary disk at the planet location.

From Fig. 9, we can infer that the planet luminosity is the main heating source in the CPD out to 0.6 au but beyond this radius the temperature is dominated by the stellar heating, either in the form of direct irradiation or in the form of photons re-emitted and scattered by the protoplanetary disk.

We conclude that the temperature of dust grains in the midplane of a CPD that is embedded in the dust cavity of PDS 70 at 34 au is determined by two heating mechanisms which dominate at different spatial locations. Additionally, we can see that the temperature profile is monotonically decreasing when the planet+CPD system is isolated but when it is embedded in the natal protoplanetary disk a reversal of the temperature happens at ~ 0.8 au. As a result of this reversal, the predicted temperature at the outer edge of the CPD is 13 K in the embedded case but only 6 K in the isolated case.

The predicted temperature reversal in the outskirts demonstrates the relevance of a self-consistent radiative transfer modelling ofa CPD embedded in the circumstellar environment. This result motivates further studies about the position of icelines in CPDs and their impact on the formation of exomoons.

Fig. 8 Surface density profile for an optically thin CPD at sub-millimetre wavelengths of 1.2 au in size obtained from Eq. (1) using the parameters shown in the lower-left corner. The dust mass is 0.007 M⊕.

Fig. 9 Azimuthally averaged midplane temperature for an optically thin CPD at sub-millimetre wavelengths around a 2.8 MJup planet. Filled circles indicate the temperature retrieved from the simulation that includes to the planet and to the star as heating sources; open circles indicate the temperature measured when the planet was not included. Vertical line indicates the planetocentric distance beyond which the star alone provides more than 50% of the total CPD heating. The blue curve shows the case of the isolated planet+CPD system where the heating from the circumstellar material is not modelled. The horizontal line marks the characteristic temperature of the protoplanetary disk at the location of the planet.

5 Conclusions

In this work we present a three-dimensional radiative transfer model that jointly explains the observed continuum emission from the PDS 70 system at 855 μm and at 1.25 μm in the Q_ϕ component of the polarised and scattered light. To this model, we then added a second source of photons that mimics the planetary mass companion PDS 70 c and the circumplanetary disk (CPD) around it. Here, we summarise our results:

• 1.

  Our simulations show that dust grains are spatially segregated by size in the PDS 70 system. The inner disk is populated by small grains with a maximum size of 100 μm and the outer disk is populated by larger grains with a maximum size of 3 mm. Thevertical settling is controlled by turbulence parameters equal or lower than 2 × 10⁻². Other bulk parameters such as the carbon content and the porosity are considered to be constant throughout the disk. These properties for the dust grains, alongside a density profile modulated by the substructures observed with ALMA, allow us to reproduce the sub-millimetre and the near-infrared polarised-scattered emission from the PDS 70 system.

• 2.

  The inclusion of a planetary companion of 2.8 M_Jup with values of the effective temperature and luminosity constrained by the work of Wang et al. (2021) (i.e. T_p = 1054 K and L_p = 4.5 × 10⁻⁵ L_⊙) does not generate any detectable signal at sub-millimetre wavelengths if the model is convolved with a Gaussian beam of size 67 × 50 mas. Since the observation shows a clear substructure around the location of the planetary companion PDS 70 c, our models support the presence of an optically thick circumplanetary disk that surrounds the planet.

• 3.

  Given our assumptions on the planet luminosity and grain properties, and using the most recent estimate for the upper limit for the size of the CPD, we constrained the dust mass in the CPD to the range 0.07−0.7 M_⊕.

• 4.

  For planetocentric distances lower than 0.6 au, the thermal structure in the midplane of the CPD around PDS 70 c is primarily determined by the planet luminosity. For larger distances, photons emitted by the star and reprocessed by the protoplanetary disks are the dominant heating source.

Although we adopted the premise of a CPD orbiting PDS 70 c, we point out that other interesting scenarios have been proposed, such as the effect of accretion streams and their ability to reproduce some of the observed features without invoking the presence of a CPD (e.g. Toci et al. 2020). However, the recent high resolution continuum observations of PDS 70 (Benisty et al. 2021) seem to favour the CPD scenario. Spatially resolved observations of the gas dynamics in the vicinity of PDS 70 c will help to complete the puzzle of this interesting planet forming region.

Acknowledgements

We thank the anonymous referee for a constructive report that improved the quality of the paper. This work is partly supported by the Netherlands Research School for Astronomy (NOVA). C.H.R. acknowledges support from the DFG Research Unit “Transition Disks” project number 325594231 (FOR 2634/1 and FOR 2634/2). C.H.R. is grateful for support from the Max Planck Society.

Appendix A Physical structure of the modelled protoplanetary disk around PDS 70.

In this section we show the vertical cut in the three-dimensional structure of some relevant quantities of our model for the PDS 70 disk. In Fig.A.1 we report the cuts in the density structure and the equilibrium temperature. In the absence of a planetary companion or any other localised substructure, these fields are azimuthally symmetric. The dust opacities for each zone are shown in Fig. A.2.

Fig. A.1 Left: Vertical cut in the density structure of the protoplanetary disk. Labels along the contours represent the logarithm of the density. The patched area shows the locations in the disk where log ρdust < −30. Right: Vertical cut of the equilibrium temperature field. Labels along the dashed contours show the logarithm of the temperature.

Fig. A.2 Extinction, absorption, and scattering opacities computed for the dust grains populating each zone of the model. The wiggles observed at λ > 300 μm are a numerical artefact product of the implementation of a continuous size distribution in the form of discrete bins.

References

All Tables

All Figures

	Fig. 1 Surface density profiles before and after the iterative procedure described in Sect. 2.1. The modified profile provides the best fit to the ALMA data. The grey area depicts the spatially unresolved inner disk due to the ALMA beam size.
In the text

	Fig. 2 Comparison between the ray-traced images retrieved from our model (right column) and the ALMA-VLT/SPHERE observations (left column). Top panels: thermal emission at λ = 855 μm. Bottom panels: Qϕ component of the polarised and scattered emission at λ = 1.25 μm. North is up,and east is to the left.
In the text

Fig. 3 Left panel: azimuthally averaged radial brightness profile from the ALMA image (black dots) and from our model (blue line). Middle panel: radial cut along the major axis of the Qϕ frames (negative radial direction is along the red-shifted semi-major axis, which is nearly coincident with the south-east direction). Right panel: sepectral energy distribution towards PDS 70; black dots correspond to the source’s photometry and the blue line is our prediction.

In the text

Fig. 4 Synthetic images including a planetary companion surrounded by a CPD. The planet was placed at a distance of 34.3 au with a position angle of 280.4° measured east from north. The three images have similar pixel scales (4 mas px−1) and were convolved with the same beam properties of the ALMA observation (see Sect. 2.1). The colour scale is set equal to the square root of the surface density to boost the signal from the CPD. North is up, east is to the left.

In the text

	Fig. 5 Left: cut in the vertical distribution of mean particle sizes in the inner zone of the protoplanetary disk. Labels along the dashed contours show the mean particle size in μm. Right: same as the firgure on the left, but for the outer zone.
In the text

	Fig. 6 Radial cuts extracted from Fig. 4. The conical aperture was centred at a position angle of 280.4° in line with the astrometric position of PDS 70 c. The aperture has a width of ± 10° around the mean position and the radial step size is 2 au.
In the text

	Fig. 7 Effect on the CPD flux due to the adopted planet mass and luminosity. In this experiment, our reference CPD mass was fixed to 0.7 M⊕. Vertical dashed lines represents the limits of the circular aperture used to measure the flux emitted by the CPD.
In the text

	Fig. 8 Surface density profile for an optically thin CPD at sub-millimetre wavelengths of 1.2 au in size obtained from Eq. (1) using the parameters shown in the lower-left corner. The dust mass is 0.007 M⊕.
In the text

Fig. 9 Azimuthally averaged midplane temperature for an optically thin CPD at sub-millimetre wavelengths around a 2.8 MJup planet. Filled circles indicate the temperature retrieved from the simulation that includes to the planet and to the star as heating sources; open circles indicate the temperature measured when the planet was not included. Vertical line indicates the planetocentric distance beyond which the star alone provides more than 50% of the total CPD heating. The blue curve shows the case of the isolated planet+CPD system where the heating from the circumstellar material is not modelled. The horizontal line marks the characteristic temperature of the protoplanetary disk at the location of the planet.

In the text

	Fig. A.1 Left: Vertical cut in the density structure of the protoplanetary disk. Labels along the contours represent the logarithm of the density. The patched area shows the locations in the disk where log ρdust < −30. Right: Vertical cut of the equilibrium temperature field. Labels along the dashed contours show the logarithm of the temperature.
In the text

	Fig. A.2 Extinction, absorption, and scattering opacities computed for the dust grains populating each zone of the model. The wiggles observed at λ > 300 μm are a numerical artefact product of the implementation of a continuous size distribution in the form of discrete bins.
In the text