© The Authors 2025

1 Introduction

Planet formation is a by-product of the star formation process. A significant mass fraction of the collapsing cloud that makes up a future star (Hueso & Guillot 2005; Pascucci et al. 2016; Manara et al. 2018; Tychoniec et al. 2018) collapses into a rotationally supported circumstellar disc (CSD) that forms planets (Najita & Kenyon 2014). The dust component that is incorporated into the CSD will initally bear resemblance to dust found in the interstellar medium (ISM). This ISM dust consists of amorphous olivine mineral and amorphous graphite components with a maximum size of ∼0.3 µm (Draine 2003), whereas the carbon- bearing species could also be organic in nature and grown from polycyclic aromatic hydrocarbons (PAHs) (Pollack et al. 1994). Beginning in the ISM, dust plays an important double role throughout the entirety of the planet formation process. On the one hand, it provides the mass for the larger building blocks of planets. On the other hand, it is responsible for cooling and compacting the gaseous envelopes of the forming planets.

The earliest stages of planet formation can be accessed via studies of the spectra of nearby star-forming regions in the Milky Way (Martin et al. 2012) and pre-stellar cores (Chacón-Tanarro et al. 2019). Those works make it clear that dust undergoes a significant and environment-dependent evolution as it passes through ever-denser mass reservoirs. Indirect evidence from class 0 (Harsono et al. 2018) and class II discs around brown dwarves (Pinilla et al. 2017) and detailed dust coagulation models (Birnstiel et al. 2012) furthermore suggest that the dust can grow to millimetres in size in a number of sources, although it remains unclear whether the growth occurs in the olivine, organics, or an icy component. Dust then can continue to form planetesimals and terrestrial planets and nucleate the cores of giant planets (Lambrechts et al. 2014; Morbidelli et al. 2015).

Cores of giant planets surround themselves with a significant amount of gas, called the protoplanetary envelope. This envelope can be thought of as an enormous extension of the protoplanetary atmosphere, which can still be physically connected to the parent CSD, particularly before gap-opening. Because the gas is unable to form a significantly overlapping spectral line continuum in the low-pressure environment of envelopes, the envelope cooling capability is dominated by the dust continuum opacities. Those opacities have to be low enough to allow proto-giant envelopes to radiate their compressional heat and reach runaway gas accretion before the CSD dissipates (Pollack et al. 1996; Piso & Youdin 2014).

The opacity of a dust mass is inversely proportional to the size of the individual particles, because larger dust distributes the same mass in a lower effective surface area and is therefore generally less opaque (Ossenkopf & Henning 1994). As a consequence, the opacity of ISM-like dust is too high to cool gaseous envelopes in time to form giants (Hubickyj et al. 2005). However, it was shown that the envelope provides the incorporated dust with the dense environments necessary to grow up to millimetre sizes; hence, the required low opacities to form gas giants can be reached self-consistently (Movshovitz et al. 2010; Mordasini et al. 2014). The results from the latter works indicate that the key cooling parameter, the Rosseland opacity, should be of the order of κ_R = 10⁻³ – 1 cm² g⁻¹, corresponding to single-grain sizes of a centimetre to 10 µm, with variations around this number allowed.

Computational models in 2D and 3D predict that the envelopes of giant planets form structures, in the form of circumplanetary discs (CPDs; Lubow et al. 1999; Kley 1999; Ayliffe & Bate 2009a; Tanigawa et al. 2012; Szulágyi et al. 2016; Maeda et al. 2022; Li et al. 2023), when simplified assumptions about thermodynamics, and in particular the cooling rates, are made. The CPD would furthermore form an environment for moon formation (Canup & Ward 2006; Szulágyi et al. 2016; Shibaike et al. 2019; Ronnet & Johansen 2020; Anderson et al. 2021). We recently were able to show in Schulik et al. (2020) that the formation of a CPD is possible at κ_R = 10⁻² cm² g⁻¹ , while it is suppressed at κ_R = 1 cm² g⁻¹ . Both our simulated giant planet envelopes accrete gas vigorously, but only for the low opacity does the accreted gas pass through the CPD in the midplane, leading to mass build-up in the CPD. The latter result has recently been understood to belong to a general class of cooling criteria (Krapp et al. 2024), while viscosity-based criteria are also being investigated (Lega et al. 2024).

We now have strong evidence for two planetary objects, PDS 70b and PDS 70c, being young, actively accreting protoplanets in the PDS 70 system (Keppler et al. 2018; Haffert et al. 2019; Aoyama & Ikoma 2019; Zhou et al. 2021). Recently, a third object, PDS 70d, of an as-yet-unknown nature, has been proposed (Christiaens et al. 2024). This system, which has an age of 3–5 Myrs and a possible CPD around PDS 70c (Isella et al. 2019; Benisty et al. 2021), will be used as an inspirational case in this work. Further analysis of the system was carried out in Bae et al. (2019), who found that the two planets are close to a 2:1 resonance and that the dust ring widths are consistent with filtration by the planets, emphasizing this system’s potential as a laboratory for planet formation theory. We have investigated the CPDs of simulated Jupiter-mass planets, whose masses are compatible with the extremes of the mass determination from Wang et al. (2020). We emphasize the notion of a ‘protoplanet’ as a planet that is observed during its ongoing formation process, which makes the PDS 70 planets interesting. In particular, the finding of a beam-resolved sub-millimetre source in PDS 70c now opens the possibility of studying the building material of protoplanets and CPDs up close, but this requires an understanding of the dust that produces the PDS 70c sub-millimetre emission.

With the understanding that CPDs need low opacities to form, the next step is to investigate possible dust components that are consistent with those values, and whether one could infer the presence of those components from spectral observations, such as those of PDS 70c. Additionally, it is of importance to ascertain the source of the dust emission flux, as there is a possibility that the observed spectrum of a protoplanet is a mix of planet and CPD contributions (Haffert et al. 2019; Szulágyi et al. 2019).

Therefore, in this work, we have used the densitytemperature structures from 3D radiation hydrodynamics simulations that were presented initially in Schulik et al. (2020). We then post-processed those simulations with a number of plausible dust populations that are consistent with low κ_R, in order to generate synthetic spectra of CPDs and assess their importance in relation to the planets. Our results indicate that if millimetresized particles dominate the particle distribution, they will act to significantly grey the spectrum of the planet and the CPD. Furthermore, spectra of objects that emit at a number of temperatures are always a sum of black-bodies. Therefore, depending on the exact geometry of emitting source temperatures, this summation can leave traces by flattening the far-infrared slope of a spectrum. This would be analogous to the process in the spectra of CSDs, which we see weak evidence of in our synthetic observations. Finally, we find that the contribution of the CPD to the total spectrum generally peaks between 30 µm and 500 µm. Our work, while not attempting to fit an exact model to PDS70c, nonetheless makes it plausible therefore that the broadness of this spectrum indicates dust growth.

This paper is organized as follows. In Sect. 2, we introduce briefly the density-temperature structures of our CPD simulations. Furthermore, we outline the provenance of the laboratory dust analogues that we use and discuss the size, fluffiness, composition, and ice coating parameters of the basic dust particles employed. The same section discusses the computational tools employed to arrive at the final data product, which are images and spectra ranging from the NIR to the submillimetre. In Sect. 3, we briefly outline the geometry of the emission surfaces for one example CPD at various wavelengths, and present the general analysis procedure. We later employ all dust opacity carriers of interest to discuss synthetic spectra and imaging. Section 4 summarizes our findings and discusses further work that is needed.

2 Methods

To generate synthetic spectra and observations, we used the radiative transfer code radmc3d (Dullemond et al. 2012) v0.41¹. This code needs density data, temperature data and wavelength- resolved opacities for all dust species involved. The latter need to be computed with separate tools for individual dust parameters from the laboratory optical constants, n and k. In the following paragraphs, we describe our choices of each of those components. Table 1 summarizes the combinations of hydrodynamic data and dust post-processing parameters used.

2.1 Optical constants, dust parameters, and radiation transport tools

Dust populations in the ISM consist of ∼0.3 µm amorphous olivines and amorphous carbonates, with additional smaller components of amorphous carbonates at particle sizes a_P ≈ 10 nm and PAHs of ≈2 nm sizes (Draine 2003; Compiègne et al. 2011; Planck Collaboration Int. XVII 2014). The PAH macromolecules, however, do not seem to be important for the spectra of protostellar discs and are absent in the presolar meteoritic record (Scott 2007), which is why we neglect this component in our work.

Moreover, interstellar dust does not seem to be a good proxy for planetary building blocks. The terrestrial meteoritic record shows that chondrules of up to cm sizes formed in the solar system continuously for 3 Myrs after CAIs (Connelly et al. 2012), and are mostly products of repeated heating and crystallization in the early solar system (Scott 2007). While their high-pressure, high-temperature processing is thought to have happend in extreme environments, such as nebular shocks of various kinds, their mass assembly must have happened prior to the thermal processing and hence implies a presumably fluffy, precursor population of dust particles similar to what is seen in cometary (protosolar) dust (Brownlee 2014; Güttler et al. 2019), only larger, at millimetre to centimetre sizes. Together with the previously mentioned indirect evidence for millimetre-sized particles in CPDs (Pinilla et al. 2017; Harsono et al. 2018) we think it is reasonable to assume that a population of millimetre-sized particles would dominate the dust mass budget in a several Myr old source orbiting a fully formed giant planet, like PDS 70c.

The exact number distribution and chemical makeup of those grown particle populations will be of paramount importance for spectral modelling. Motivated by analytical models (Birnstiel et al. 2012) and observations of the ISM (Draine 2003; Compiègne et al. 2011) we use three mixes of dust distributions in order to match the two different average Rosseland opacities from our hydrodynamic simualtions. The average Rosseland opacity κR is defined by $${\kappa }_{\text{R}}={\left(\frac{{\displaystyle {\int }_0^{\infty }d}\lambda {(\kappa (\lambda ))}^{-1}\frac{\partial B}{\partial T}}{{\displaystyle {\int }_0^{\infty }d}\lambda \frac{\partial B}{\partial T}}\right)}^{-1},$$(1)

where κ_λ = Σ_i κ_i(λ) are the non-averaged opacities with contributions from all dust species and B = B(λ, T) is the Planck function for a black-body of temperature T . As we use two simulations with hydrodynamic data for constants κ_R = (1,0.01) cm² g⁻¹, those values must be matched at least approximately by the dust populations in order to ensure self- consistent modelling.

To compute the κ_i (λ) one needs the real and imgainary optical constants, n(λ) and k(λ), of the material in question. We used amorphous olivine, organics, and water ice, with the optical constants taken from the St. Petersburg–Jena–Heidelberg database of optical constants (Henning et al. 1999; Jäger et al. 2003) (the ‘new’ dataset² with the ‘normal’ olivines³). We converted the optical constants into κ_i (λ) in the simplest case of spherical grains by specifying the particle size, a_P, and particle density, ρ_P, and computing their Mie-scattering properties with the Bohren & Huffman (1983) code.

The resulting opacity functions for the two key dust sizes in this work, are computed by the Bohren & Huffman (1983) code and can be seen in Fig. 1 (Left). The millimetre-sized grains display very grey spectral behaviour throughout the entire wavelength range; hence, we can expect the total spectrum of the CPD to be a convolution of black bodies. The small dust is nevertheless also fairly grey, as can very clearly be seen in the absence of the otherwise no table 10 µm feature for the olivine dust component. In Fig. 1 (Right), we show that this disappearance of the 10 µm feature is due to a dramatic transition of spectral flattening that happens between 1–10 µm of particle size, as would be expected from Mie theory in the interference regime, 2πa ≈ λ.

The dust mixes we used were constructed as follows. Each dust mix consists of amorphous olivine and organics that follow a scaled Draine (2003)-like size distribution. This distribution consists of several sub-distributions. The sub-distribution representing the largest particles is between 0.1 and 1 times the maximum size, with a MRN-like distribution between them. We refer to the latter as the major size for a mix, a_major . This major sub-distribution exists for both amorphous olivine and organics. The Draine (2003)-distribution contains an another sub-distribution of smaller, amorphous carbon species, for which we again used organics, which has a Gaussian radius size distribution centred at 1/30 of the major size. We refer to this component as the minor particle size, a_minor, which contains 1/30 of the large dust particles mass. Finally, water ice exists exclusively as coatings on the organics and olivine grains and distributes its total mass and volume equally on their surfaces.

We employed three types of such mixes in this work, each characterized fully by the preceding description and the statement of its a_major . The large dust mix, or ‘L’ in Table 1, has amajor = 1 mm and has been primarily designed to fit the average κ_R = 0.01 cm²g⁻¹ in our radiation hydrodynamic simulation. The small dust mix, or ‘S’, has a_major = 10 µm, and the ‘ISM’ mix has a maximum size of 1 µm.

The total dust mass that this scaling correpsonds to is 2 × 10⁻⁵m_⊕, deduced from the gas mass in our nominal CPD of 2 × 10⁻³m_⊕ with f_DG = 0.01. In our later simulations, runs 5–9, which emulate more mature CPDs (see also Table 1), both masses are increased artificially by a factor of 100, by upscaling all volumetric gas densities. The increased dust masses correspond to enough solid material to build an analogue of the Galilean moons.

We generated the dust density, ρ_i,D, of each dust subdistribution, i, from the gas density, ρ_G (described below), using the dust-to-gas ratio, f_DG , from Table 1, and a species fraction, f_i , so that it is $${\rho }_{i,\text{D}}={\rho }_{\text{G}}\times f_{\text{DG}}\times f_i\times 1/100,$$(2)

where the index, i, encapsulates the species; that is, i ∈ {Ol, major; C, major; C, minor; H₂O}. The numbers f_i are taken from Semenov et al. (2003) which are consistent with observations of the ISM, hence for dust of approximately micrometer in size, but lacking laboratory data on the relative growth speeds of dust compounds, we take those fractions for all dust mixes. Another reason to employ those mass fractions is because Draine (2003) did not measure water ice abundances, which we need for our modelling. This also implies that we do not prescribe any differential settling, which is consistent with the high a ≈ 10⁻² that was used to generate the hydro data.

The mass fractions are $f_{{\text{H}}_2\text{O}}=0.555,f_{\text{Si}}=0.264$, and f_C,major = 0.353, and we estimated from Draine (2003) that f_C,major = 30 × f_C,minor. Furthermore, the sublimation temperatures are $T_{{\text{H}}_2\text{O},\text{\hspace{0.17em}sub}}=155\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}K}$, T_{si, sub} =1500 K, and T_{C, sub} = 410 K.

The optical properties of the water ice coatings in our model are ignorant of their respective substrate. This does not reflect reality (Ossenkopf & Henning 1994), as an electromagnetic wave penetrating into an icy dust coating will interact with its substrate, if the coating is thin enough, but we nevertheless maintain this simplification here. In our implementation, the ice spectrum simply replaces that of olivine and organics in regions where $T<T_{{\text{sub,\hspace{0.17em}H}}_2\text{O}}$, which is the sublimation temperature of water ice. Hence, cold regions in which the water has not sublimated can be dominated by water ice emissions, depending on the exact model.

Computing the average Rosseland opacity now with those mixes of the specified sizes, the result would always be a factor of three times larger than the κ_R from the hydro data. One could fix this by adopting even larger particles, but this would erase all spectral features. Due to this reason we adjust the dust-to-gas ratio by this mismatch factor of 3.

The overall generation of spectra with radmc3d is performed without the use of scattering. A detailed investigation of scattering effects for a different set-up has been presented in Szulágyi et al. (2019). Their analysis concludes that scattering would eradicate features between 1–10 µm efficiently and overlay it with a scattering slope, but is naturally very weak for longer wavelengths. Because in this work, we are particularly interested in the continuum emission properties in the near-infrared to sub-millimetre range, we neglect scattering effects.

We treated gas opacities only in a very simplified manner, to imitate their effect as a background opacity source in very hot regions where all dust has evaporated. To this end, we set the gas opacity function to a constant grey value, κ_G(λ) ≡ 10⁻⁵cm²/g. However, for all practical purposes we find that the inclusion or exclusion of the gas in this manner has no effect on the spectra because the spectral signal is generally dominated by cool dust. Regions where hotter dust and gas emission might play a role, due to the sublimation of all dust components, are however hidden from the observer, as all τ = 1 surfaces in our models lie exterior the olivine sublimation front (cf. Fig. 2, dark blue line).

All this information was fed to the open source code RADMC3D (Dullemond et al. 2012) to generate the objects’ spectra. We employed this code without a stellar irradiation source, in which the source function for the photons is purely given by the temperature field from our radiation hydrodynamic simulations. The code then solved for the radiation transport through the density and opacity fields, fed into the code as separate input. The optical depth in each cell was chosen such that the optical depth of the dust mixture in each cell corresponds to the optical depth in our radiation hydrodynamic simulations. We shall now proceed to describe the radiation hydrodynamic data and how it was processed.

Fig. 1 Computed opacity functions for water ice, olivine and organics for individual components of single particle sizes (Left), scaled according to the available dust mass according to Eq. (2). The millimetre-sized particles make up the ‘L’ particle mix, whereas the 10 µm particles make up the ‘S’ particle mix (see also in the main text and Table 1). In general, it becomes evident that 10 µm particles are already fairly grey, except for the water ice. Furthermore, we show how this greyness, and particularly the disappearance of the prominent 10 µm feature for olivine, results from step-wise dust growth (Right). The opacity functions for the solid lines on the right are computed for pure olivine at single sizes only. The step-wise increase in size highlights the disappearance of the 10 µm-SiO feature in our spectra for dust sizes ≥10 µm, as then the size parameter 2πa/λ enters the interference regime of Mie-scattering.

Fig. 2 Cylindrically averaged radial-vertical map of κR computed only for ice and olivine components, centred on the planet for sim 1. This highlights the approximate 3D structure of ice and olivine sublimation ellipsoids in the envelope, as is seen via the corresponding opacity drops. Overlying black contours denote the gas density in g cm−3 to highlight the disc-structure of the hydrodynamic solutions. The Hill sphere is marked with a red circle. The actual τ = 1 surfaces at λ = (5, 30 , 300) µm computed with the full solution for sim 1 and viewing from z = +∞ are overlaid in dark blue, light blue, and cyan. Dashed colours denote the mirrored continuation of the surfaces below the disc midplane. The τ = 1 surfaces help to interpret spectral information.

2.2 Hydrodynamic data

In a previous work (Schulik et al. 2020), we have carried out global 3D radiation hydrodynamic simulations of the gaseous envelopes around Jupiter-mass planets. For an overview, we shall briefly summarize some of the most important points from that work. Our protoplanets are represented by a gravitationally smoothed potential. In units of r_H = 5.2 AU (m_planet/(3m_star))^1/3, the depth of the smoothed potential is r_s = 0.1 r_H, outside of which the potential is the planets real gravitational potential. The simulation domain extends well beyond 5r_H into the CSD. Gas is accreted from the CSD and onto the planet and the CPD via the action of gravity, angular momentum conservation and cooling. The dynamics in the disc is solved via the Fargo-algorithm (Masset 2000) in the code FargO CA (Bitsch et al. 2014; Lega et al. 2014). We employ viscous heating throughout the simulation domain with a constant physical viscosity, proportional to a viscous α of α ≈ 3 × 10⁻³ – 10⁻². The cooling was computed via radiation transport in the flux-limited diffusion approximation (Levermore & Pomraning 1981) with the flux-limiter by Kley (1989), for which κ_R, together with the gas density, ρ, and a fixed dust-to-gas mass ratio of ϵ = 0.01, plays the role of determining the optical depth, τ, of the dust and gas mix. Our previous analysis (Schulik et al. 2019) shows that a minimum resolution requirement must be observed in order for all those effects to play together and let accretion continue in a physically self-consistent manner, without numerical overheating of the envelope.

Those simulations show that a low κ_R is crucial for the formation of structures in protoplanetary gas envelopes and for the gas to reach Keplerian rotation. This can lead to a flattening of the gas into CPDs and rotation velocities of the gas surpassing υ_θ/υ_Kepler > 90% for values of κ_R = 0.01 cm²/g, while for high opacities of κ_R = 1.0 cm²/g the gaseous envelope remains fairly spherical and unable to form a rotating disc consistent with Szulágyi et al. (2016). We use the density-temperature data from the runs with those two constant Rosseland opacities for our synthetic observations. Nominally, the high κ_R simulation will be matched with the ‘S’ dust mix (sim 1) and the low κ_R will be matched with the ‘L’ dust misx (sim 2).

The structure of κ_R in our post-processed radmc3d runs is shown in Fig. 2 and reflects the temperature structure of the accreting const-κ_R gas giant envelope: vertically, above the planet, the accreting gas shocks and increases the temperature locally to ~500 K, leading to a deviation in the opacity and density contours. This effect is more extreme in the strongly cooling simulation with the low average κ_R = 0.01 cm²/g, which can be seen in Fig. 3. There, a hot funnel of gas empties out the icy dust above the accretion shock, making it transparent to escaping radiation. The relative patchiness of the opacity structures seen in Fig. 3 results from temperature fluctuations due to many complex adiabatic compression features and colliding flows in that simulation. In the midplane for both simulations, the density and opacity contours are aligned, as no disequilibrium radiation from the shock is present in that region. We note that in Schulik et al. (2020) we also ran a simulation with non-constant opacities given by Bell & Lin (1994), which would reproduce this radialvertical structure of κ_R more closely. However numerical speed and stability issues favoured our usage of the const-κ_R hydro data, also noting that the gas dynamics in the outer envelope does not change significantly.

To quantify the influence of the CPD on the overall spectra – that is, to assign ‘disc-dominated’ or ‘planet-dominated’ wavelengths – we first generated a spectrum of the pair CPD + planet and then substracted the CPD region from the data in a second step. Both those steps are performed in a face-on view (from the top of the disc) because otherwise one would peer sideways into the planet after subtracting the CPD, which would add excess contributions to the spectrum. This approach is in contrast to the synthetic image in Fig. 5, which is taken at a 50° inclination to mimic the view on PDS 70c.

The net effect of lowering κ_R on the optical depth via dust growth can be achieved by either increasing the dust sizes or by decreasing the overall number of dust grains. Hence we also process the κ_R hydrodynamic data with the ‘S’ and ‘ISM’ dust grains, but with 100 times (1000 times) lowered dust density, giving rise to sim 3 (sim 4). We call sims 1–4 the ‘early’ type sims due to their low overall masses. As is outlined in Schulik et al. (2020) our simulations feature only the initial few orbits of the formation stage of a CPD. Therefore our CPD dust masses are very low (of the order of a few 10⁻⁵ m_⊕). To obtain spectra of 100 times more massive CPDs, like those estimated in Isella et al. (2019), we simply scale all gas densities from sim 2 up by a factor of 100, and again use the different possibilities for obtaining a low κ_R, giving rise to sims 5–7, which we call ‘late’-type sims.

Additionally, under the hypothesis that late-type CPDs might be significantly cooled down compared to our simulations, we investigate another set of ‘late & cool’ CPDs, which feature a factor of 3 reduction in temperature for all cells. This procedure gives rise to sims 8 & 9. The late and cool sims are only processed with ‘L’ and ‘ISM’ dust mixes.

To simulate the effects of sublimation, we assume that ρ_i,D = 0 if the local temperature in a grid cell is above the sublimation temperature for that dust species T_i,sub. Ice coating sublimation is approximated by setting each dust species to their nominal values ρ_Ol, pc,major, ρ_C,minor as was described above, and ${\rho }_{{\text{H}}_2\text{O},\text{D}}=0$ at dust temperatures of $T>T_{{\text{H}}_2\text{O},\text{sub}}$. The temperature for all dust species is taken as equal to the gas temperature.

Fig. 3 Same as Fig. 2 but for the low-opacity hydro-simulations and sim 6, which has the same optical depth as sim 1 after a ×100 density increase. Again, the background colour is the Rosseland opacity computed only with olivine and water ice to highlight compositional differences. The τ = 1 emission surfaces were again computed with the full opacity model including all species, as was used in the paper to compute spectra.

Fig. 4 Intensity map of a synthetic observation of the planet and CPD at 5 µm. The image is as is to be expected from the computed 5 µm-emission surface in Fig. 2. The CPD is not visible, as at those wavelengths water ice and organics as opacity sources have sublimated away. This image is also representative of any simulation set-up where we find planet-dominated spectral regions, which can also happen up to the sub-millimetre. It is then usually the black-body-tail from the planetary contribution overshining the CPD signal. The latter is then weak, as the flux is Fλ ∝ κ(λ)Σ, with Σ being the surface density, and κ(λ) ∝ λ−2 in the far-infrared. Increasing the contrast to arbitrary, unrealistic levels would still make the CPD contribution visible (see text for more details).

3 Synthetic observations

3.1 General considerations – Emission surfaces leading to spectra

We begin our analysis with an investigation of the emission surfaces for simulation 1. The emission surfaces were computed by integrating the optical depth from infinity in a top-down view onto the CPD per individual wavelength. The surfaces were then found where the optical depth becomes larger than unity. They are a useful tool to understand the contribution of the CPD to the total spectrum, as well as the overall shape of the spectrum. This is due to the fact that an existing emission surface will pose a very bright area, radiating as black-body, whereas regions without an emission surface – optically thin regions – are not very bright Mihalas & Mihalas (1984).

In Fig. 2, we show a background colour map of κ_R, computed with only the olivine and water ice components for simulation 1, in order to clearly show the compositional transition at the water and olivine sublimation lines in the envelope of that particular giant planet. It becomes evident that the water ice line is in general more shaped like an ice ellipsoid. We overlay the gas density contours in black, which illustrate that the hydrodynamic data for κ_R = 1 cm²/g has difficulty forming a CPD. The hydrodynamic data used for sims 2–9 forms a flat CPD (which we show in Fig. 3 for completeness), but we found sim 1 to be an ideal illustrative example to study the emission properties of the giant planet envelopes.

Next, we show the emission surfaces at 5,30 and 300 µm now computed with the full opacity model, that is with all dust species included, in Fig. 2. It becomes evident that the emission surfaces (particularly around 30 µm) can follow a disc-like structure, owing to a significant opacity contribution from water ice and carbonates. The emission surfaces are representative of the final synthetic image at a given wavelength.

The blue line, representing the emission surface at 5 µm in Fig. 2, will be planet-dominated, as for temperatures corresponding to this wavelength $\left(T\approx 600\text{\hspace{0.17em}K}>T_{\text{sub},\text{C}},T_{\text{sub},{\text{H}}_2\text{O}}\right)$ both organics and water have sublimated. This results in a rather unspectacular appearance in the 5 µm image, Fig. 4, and leaves the silicate disc outside of r > 0.1 r_H transparent, as the silicate has an opacity minimum at this wavelength, as was previously seen in Fig. 1 (Left).

The light blue line, showing the emission surface at 30 µm, results in a more complex combination of all opacity carriers involved. As this emission surface clearly cuts through the water ice sublimation ellipsoid, a synthetic image at this wavelength merits some further study. In Fig. 5, we show such a synthetic image.

This image is generated at an angle of 50° similar to PDS 70c, in two different colour maps, to discuss the contrast of features seen in them. A typical ‘heat’ colour map, shown in Fig. 5 (Left), carries intuitive information pertaining to the relative brightness of features, but does not clearly show the CPD without hiding surrounding features. Therefore we also show in Fig. 5 (Right) an adjusted colour map, which can more clearly guide the expectations for the CPD contribution to the 30 µm spectrum. The CPD shows at its inner edge the water sublimation front, decreasing its brightness (marked as a white arrow). However, radially outwards from the planet, behind the water ice line, the CPD has a brightness of about 30% of the maximum value and a much larger surface area. Hence, for this wavelength we would expect the spectrum to have significant CPD contributions. This is not always given by the water ice, as olivine from inside the water ice line can also contribute importantly to the CPD contribution. Additionally to the planet and the CPD, other dense flow features radiate at this wavelength, such as the circumstellar spiral arms. Those spiral arms additionally show their complex 3D structure due to the inclination of the image.

The cyan line at 300 µm interestingly leaves the CPD midplane mostly transparent. This is because sim 1 uses the ‘S’ particles, which exhibit the opacity slope in this wavelength region, as is seen in Fig. 1 (Left). Only a weakly radiating CPD component from an emission surface below the midplane contributes to the spectral intensity at this wavelength. This is an important point of our analysis, as it explains why the CPD emissions can be dominated by the planet in the sub-millimetre for small particles, but the CPD contribution can be very substantial when using large particles at those long wavelengths.

As the 30 µm wavelength image turns out to exhibit significant CPD features, we now first continue analysing its compositional build-up, before turning to a joint analysis of composition and spectra.

Fig. 5 Synthetic intensity maps for the planet and CPD for sim 1 at 30 µm seen at 50° inclination, to mimic PDS 70c. Those maps assume a spatial resolution far superior to ALMA, but serve the purpose to illustrate the intensity contributions to the spectral flux at this wavelength. On the left we show a intuitive heat map, on the right with a more contrasted colour map. On the right, one can distinguish several regions according to their emission intensity I : Emission from the innermost planet, originating from olivine, I ≈ 1, the CPD with the water ice line I ≈ 0.1–0.6, also marked with the white arrows, the outer edge of the CPD I ≈ 0.03, various disc structures including spiral arms at I ≈ 0.01. The arrow marks the water ice sublimation line, which through its associated drop in opacity reveals some inner parts of the envelope.

3.2 The full model – Contributions of individual dust components to the emission surfaces and the spectrum

We next investigated the relative importance of the individual dust species in forming the emission surfaces and spectra. To do this, in Fig. 6 (Left) we show exemplary emission surfaces at 30 µm, and how they successively build up as we add the dust components, for the ‘S’ mixture in sim 1, which has maximum dust sizes of 10 µm. This is complemented by the spectrum of the entire object integrated over varying fractions of the Hill-radius, in Fig. 6 (Middle), which inform about the relative importance of the surface area with regard to the spectrum.

Olivines, sublimating at the highest temperatures, form the dominant opacity source near the planet. The olivine spectrum exhibits a peak shortwards of 10 µm, because the olivine opacity there decreases (see Fig. 1) and hence opens a transmission window into higher temperature regions of the protoplanetary envelope. It should also be noted that the olivine contribution remains important in this spectrum when including more and more surface area. This can be seen in Fig. 6 (Middle): when going from the dotted red line to the dashed red line to the solid red one, it becomes evident that sim 1 exhibits an olivine-dominated spectrum.

After adding water, the emission surface profile is only slightly changed inside the water-ice ellipsoid, due to contributions from water above the planet. Outside the midplane water-ice ellipsoid, there are strong contributions evident from the water ice, which then extends the emission surface nearly out into the CSD. In the spectrum, as water can exist only at cold temperatures, there is a 20% increase in emissivity at wavelengths longwards of 30 µm, but no strong water features are visible. This is due to the enormous surface area of the central olivine blob, which is outshining all the other contributions (for the same analysis of our low κ_R-hydro data, but with the same particle sizes and optical depth, see Fig. 6 (right)).

With the carbon added, there is a fairly uniform increase in optical depth, due to the greyness of the carbon spectra. Interestingly, the carbon sublimation at 410 K leaves only a very minuscule imprint on the emission surface, which may be surprising, and the spectrum overall. The main effect of the carbon on the spectrum is to move it to slightly colder temperatures and longer wavelengths, and provide extinction for the underlying olivine bump in the spectrum.

The general build-up of any spectrum will follow the same principles that are discussed here. Hence, for a low CPD contribution in any part of the spectrum, one can be sure that the emission surface will be similar to the red curve in Fig. 6 (left) and vice versa. We now turn to compare the ‘S’ and the ‘L’ dust mixtures with their corresponding hydrodynamic models.

Fig. 6 Detailed impact of including an increasing number of opacity carriers and radial integration regions onto the optically thick surface and the spectrum. Left and middle: effects of successively adding opacity carriers into sim 1. We show the emission surfaces in real space at 30 µm on the left, again integrated from z = +∞, with the radial distance on the x axis and the vertical distance on the z axis. The resulting full spectrum is shown on the right. The olivine (Ol) sets the dominant part of the spectrum, with the contribution from the water ice being minuscule mainly due to being outshone by the olivine close to the planet. Once the organics components are added (Cmajor and Cminor), the central emission surface is lifted up to colder temperatures (longer wavelengths), leading to a weaker signal in the NIR through extinction. Interestingly, the carbon greys out its own sublimation signature in the emission surface, whereas the sublimation line is clearly visible for water in the emission surface, but only weakly in the spectrum. When measuring the CPD contribution to the spectrum, we subtract everything outside of 0.1 rH from the simulation domain. This process can be seen on the right side, when following the dashed blue line as it changes into the dotted blue and then the solid blue line. It becomes evident that this spectrum, when all dust ingredients are included, gains a lot of flux for λ > 10 µm, while its gain is minor for λ < 10 µm. This leads to the conclusion that for λ > 10 µm this model is CPD-dominated. The same can be repeated for the olivine-only spectra. Those inform us that already the hot olivine blob seen on the left can be counted as CPD contribution. Right: same as on the left side, but for sim 6, deriving from the low κR-data, but having the same optical depth as sim 1 through a ×100 density increase. Water emission from the CPD plays a larger role in this spectrum, particularly between 30 µm and 100 µm.

3.3 Early type spectra – Sim 1 and sim 2 and CPD contributions

We now present an analysis of fractional flux emitted from the planet compared with that from the CPD for sims 1 and 2. Those two simulations, as was discussed earlier, were generated with average κ_R differing by a factor of 100, hence their average optical depths will vary accordingly, and so will their relative CPD contributions to the spectrum.

As can be expected from the preceding discussion, the CPD emission will contribute significantly to the total emission at wavelengths at which the CPD region is optically thick. Hence, when the emission surfaces cut through dense regions, including the hot midplane, we find CPD-dominated spectra. In Fig. 7 (left), we see how the CPD contribution dominates the blue curve due to the high κ_R contribution from small particles. A grey envelope around the plotted spectrum indicates the CPD flux contributing 10% or more of the total spectrum, and a thick coloured envelope around the spectrum indicates that the CPD is a dominant contributor (>50%) to the spectrum. It becomes evident that the CPD is generally dominating the spectrum long- wards of 10 µm until ~300µm, as should be expected by our preceding discussion. Longwards of 300 µm, the CPD becomes transparent due to the rapidly falling Mie-slopes in the opacity functions of small particles, and the corresponding disappearance of emission surfaces. The spectrum shortwards of 10 µm is planet-dominated, due to the minimum in olivine opacity, and the fact that the opacity carriers in the CPD evaporate away at higher temperatures.

Simulation 2, having a low average κ_R and the corresponding ‘L’ dust population with a maximum size of 1 mm, is shown as the solid red curve. As was seen previously, in Fig. 1, particles of this size are generally very grey, and different dust species only differ in their constant opacity level, with water ice having the strongest contribution, then the carbonates, followed by olivine. However, due to the overall lower level of opacity, the CPD is transparent at all wavelengths, as can be seen in Fig. 7 (Left, red curve).

There is still a weak (<10%) contribution to the overall width of this spectrum from some radiation in the CPD, as can be seen in Fig. 7 (left), through comparing the red curve with the width of black-body (3). Hence, it can be dangerous to assign the total spectral flux to a CPD, only because the CPD is, or is believed to be, transparent.

3.4 PDS 70c and black-body fits

In Fig. 7, we also show the data available at the moment from broadband photometry and sub-millimetre imaging for the planet and the putative CPD around PDS 70c (for details on the conversion process of data, see Appendix A). The wide spacing of the PDS 70c data seems to already be indicative of an at least two-black-body solution, when comparing black-body 1, which connects the L′-band and sub-millimetre point of T = 300 K, and black-body 2, connecting the Ks and L′ band points with T = 850 K. It seems, however, implausible that a CPD would be fitted by only one black-body.

Additionally, the planetary black-body fit can inform us about the true planetary radius (not the emission-surface radius in our hydro models), by converting the integral of the black-body flux into a surface area, which would correspond to R_P ≈ (4.5 ± 0.5)R_Jup. This radius is too large to be explained by current formation models (Mordasini et al. (2017)) at the age of PDS 70 for any planetary mass, and at odds with the findings in Wang et al. (2020). Further, Christiaens et al. (2024) caution that the nearby bands (Y, J, H) to the ones used in this work (Ks, L′) suffer from contamination due to scattering stellar light. A hotter planetary temperature, an even stronger dusty contribution to the spectrum than inferred here or stellar scattering are all possibilities to solve the discrepancy in planetary radius.

Hence, we shall next investigate how particle growth influences our model CPD spectra, with a focus on the spectral slope in the far-infrared. Once this is understood, we discuss whether it is possible to improve the planetary Ks and L′-band black-body fit by substracting a reasonable CPD model.

Fig. 7 Comparison of the spectral shape in the two main hydro simulations and studying the effects of dust growth. Left: spectral energy distributions for sims 1 and 2 with all dust species included and normalized to the distance of PDS 70, 113 pc. Furthermore, we plot the fluxes for PDS 70c in the Ks (2.18 μm), L′ (3.8 μm) and 855 μm bands reported by Haffert et al. (2019); Isella et al. (2019). For conversion details and flux errors, see Appendix A. The observational data as well as our model spectra are compared in width to black-bodies (BB), plotted as thin, dashed curves. There are BB (1) (T = 300 K), BB (2) (T = 850 K) and BB (3) (T = 1000 K). PDS 70c fails to be modelled by BBs (1) or (2) and due to this width, it reveals to possess a CPD component. Our model data compared with BB (3) shows it to have multiple BB contributions. Grey shading indicates a weak CPD contribution, i.e. FCPD > 0.1FPlanet and thick lines indicate a dominant CPD contribution, i.e. FCPD > FPlanet- Right: spectra for the low κR-simulations, emphasizing the effect of dust growth on the spectra. To keep the optical depth constant for the smaller dust particles with higher opacities, we adjusted fDG downwards accordingly. For all wavelengths, the CPD is transparent and the emission surfaces are convex features that are situated around the planet, similar to the olivine-only case in Fig. 6 (left). The comparison of spectral width to the black-body (dotted black line) shows that for the millimetre-particle case, the spectrum is a convolution of black-bodies over the 3D structure, while the origin of the width is less clear for the smaller particles. Furthermore, the growth from ISM-like towards millimetre-sized particles 1) greys the spectrum and 2) produces a FIR-tail that is significantly shallower than the expectation from single-temperature, single-small-particles, which would be a ~λ−5-slope.

3.5 Early-type spectra – Varying dust sizes, while keeping τ constant for sims 2–4

We focus now on the low κ_R simulation in order to eliminate complications due to changing density-temperature structures, when changing κ_R. We stress the possibility of generating the same optical depth, τ, at low κ_R with particles of smaller sizes, despite their higher opacities. For this, we decreased the dust-togas ratio accordingly, in order to leave the approximate optical depth overall the same. Specifically, this means that when using the 'S' mixture with a maximum size of 10 μm, which has a ~100 higher opacity, the dust-to-gas ratio was adjusted to 10⁻⁴, and accordingly even lower for the ISM dust.

In Fig. 7 (right), the changes in the spectrum when growing the dust particles becomes evident. Due to the low density and low average opacity the CPD remains transparent for all dust sizes, and the emission becomes planet-dominated. The 'ISM' particles of 1 μm size show essentially a pure olivine spectrum in emission, but its contrast is greyed out as for sim 1 through carbon. In the olivine-opacity minimum between of 3–10 μm (as is seen in Fig. 1), the entire envelope is transparent, as the opacity there is low enough to reach the olivine sublimation line. The slope of the spectrum in the FIR/submillimetre region is in general the sum of all radiating black-body contributions B_λ weighted with their opacities, λF_λ ∝ Σ,κ_λλB_λ(T) (e.g. Backman et al. 1992). A single, optically thick black body would have λB_λ ~ λ⁻³ and with the opacity slope for small particles being κ_λ ~ λ⁻² in that region, it should be λF_λ ∝ λ⁻⁵. As the spectrum for the small particles here measures however λF_λ ∝ λ^−43±0.1, we can attribute this flatter slope to the summation effect of black bodies to the spectrum.

For the dust that is 10 μm in size, the carbon plays an important effect in greying out the olivine minimum shortwards of 10 μm. This renders the entire dust mix near the planet optically thick and the olivine features appear only very weakly as a bump in the emission. The spectrum has changed its slope, in the FIR, with a break at ~300 μm, corresponding to the change in opacity slope at ~100 μm, as is seen in Fig. 7 (Right).

Finally, the millimetre-sized particles show a featureless convolution of black-bodies, with slope slightly steeper than λ⁻³, of about λ^−3.2±0.1. This is because at the sub-millimetre level the millimetre-particle opacities show a weak slope themselves.

3.6 Late-type spectra and a fit for PDS 70 c

We now turn to investigate the spectra of late-type CPDs, for which we have used the low κ_R hydrodynamic data, and multiplied the density by 100, in order to match the masses in our simulations with dust mass estimates for PDS70c from Isella et al. (2019). The first three simulations use the temperature structure given by the hydro data (sims 5–7), while we also artificially decrease the temperature in the entire domain by an arbitrary factor of 3, in order to imitate a massive, cooled down disc (sims 8 and 9).

From Fig. 8 (left), one can read off that dust growth leads to an increasing CPD contribution to the total spectrum, mainly due to the formation of extended emission surfaces, as is seen in sim 1, which can be attributed again to the greyness of the mm particles. This is an important finding, as it shows that our simulation data for the low κ_R , which actually forms CPD structures with important Keplerian rotation, can also be detectable in the late stages of their evolution.

The FIR slope of the spectrum evolves weaker than for the early spectra, from around λ^−3.5±0.1 for the ISM particles to λ^−3.1±0.1 for the millimetre particles. The less clear distinction between small and large particles stems from a stronger contribution of black-body mixtures to the spectrum.

The ‘cold’ simulations naturally have less total emission, but the disc contributions remain dominant longwards of 30 µm for the millimetre-sized particles (sim 9). From the preceding discussion, it is not surprising that the ISM particles (sim 8) are essentially planet-dominated, but this is strengthened due to the density structure for the following reason. We seed opacity carriers in cold regions that are proportional to the gas density; this gas density is, however, severely depleted due to the ongoing accretion process. Therefore, although the number of ice particles is enhanced in the cold case overall, the optical depth does not increase significantly above the CPD due to ice carriers.

Finally, it seems that we can obtain a fit to the broadness of the PDS 70c data with a two-component solution, comprising our CPD data for the millimetre-sized particles and for a planetary black-body. From Fig. 8 (Left) it becomes clear that the spectra for the CPD + planet for the 10 µm (sim 6) and 1 mm sized particles (sim 5) are natural candidates to fit the width of the PDS 70c spectrum. We emphasize here that this is an immediate result of our mass increase that we performed in order to be consistent with the dust mass estimate for this source.

An inspection of the 10 µm dust-size model spectrum in Fig. 8 (Left, solid blue line) reveals that this spectrum is planet- dominated at 855 µm and is therefore inconsistent with another planetary black-body, which would be needed to fit the Ks and L′ datapoints. Only the millimetre-sized particles (sim 5) possess the required sub-millimetre component to fit the spacing between the NIR and sub-millimetre data, without the need for another planetary black-body.

If we take the millimetre-sized late model, subtract our model planet and instead insert a more realistic planetary blackbody into it that would fit the Ks and Lp datapoints we would however still remain with an unrealistically large planet. A smaller, more realistic planet would require a stronger CPD contribution in the L′ band, which we found challenging to model. A natural way to increase the CPD contribution in the L band would be an improved version of our hydrodynamic data, because, as we outlined in Schulik et al. (2020), our modelling overestimates the planetary contribution.

Fig. 8 Spectral shapes of simulations with nominal and increased dust masses. Left: late models with a ×100 times increased dust mass show significantly increased CPD fractions in their spectra (solid lines), due to the formation of emission surfaces in the CPD. We keep the previously introduced black-bodies for very rough expectations on the planet and CPD emissions. FIR slopes are again black-body-like for millimetre-sized particles, while the overall steepness is not as dramatic as for the early-type small particles, making a fit more difficult. However, the width of the PDS 70c data can be approximately fitted when considering larger grains. The artificially cooled down envelopes (dashed lines) are too cool to fit the observations and particularly undershoot the ALMA measurement. Right: degeneracy check for the fits on the left. We attempted verify whether increases in dust mass of smaller particles can also fit the seen width, which we can confirm here. Hence more data is needed to distinguish between mass and particle size. Curves are plotted without the indicators for the CPD fraction, as the viable fits all are optically thick in the CPD region and hence the sub-millimetre signal is actually emitted by the CPD.

3.7 Uniqueness of the fit

In order to determine the robustness of the millimetre-sized particle fit, we attempted a number of other solutions to fit the width of the PDS 70c spectrum, as is shown in Fig. 8 (right). The submillimetre datapoint can be approximately fit when increasing the optical depth of the other particles sizes at 855 µm as well, showing that there is a linear degeneracy between particle sizes and mass.

The attempt to fit the data with the high opacity CPD from sim 1 fails for all attempts of mass/particles combinations. Thereason for this is that one always overshoots the near-infrared datapoints, indicating that this envelope is simply too hot to fit the data.

4 Summary and discussion

In this work, we have used data from two radiation hydrodynamic simulations that treat heating and cooling self-consistently (for technical details, see Schulik et al. (2019) and for the CPD simulations Schulik et al. (2020)), in order to generate their plausible spectra via post-processing with the opacities of various reasonable dust populations. The two hydro-simulations use two different constant average Rosseland opacities κ_R = (1.0, 0.01) cm² g⁻¹, but have otherwise identical parameters, in order to study the structure formation of gaseous envelopes around Jupiter-mass planets (see Table 1). We note that while early mass determination (Keppler et al. 2018) indicated planet masses for PDS 70b as m_b ≈ 5-8m_Jup, astrometry-based mass determination (Wang et al. 2020) results in lower planetary masses, towards m_b ≈ 2–4m_Jup and m_c ≈ 1-3m_Jup for PDS 70c.

Due to the nature of our flux-limited diffusion approach, which was used to generate the hydrodynamic data, it was not possible to obtain spectral information directly from those simulations. Therefore, we chose to match the average optical depths found in the hydro-simulations with rescaled variants of the dust properties known from the studies of the ISM (Draine 2003; Compiègne et al. 2011) and meteoritic as well as cometary material (Scott 2007; Brownlee 2014; Güttler et al. 2019). Motivated by those, we scaled the ISM distribution up so that the maximum particle sizes are between 1 µm and 1 mm, with the minimum sizes given by the corresponding Draine (2003)-distribution.

From the dust properties, we computed the wavelength resolved opacity functions, κ(λ), which are necessary to generate spectra and synthetic imaging. The geometric optics region of κ(λ), where those functions are near-constant, was then used to match the sum of all mass-weighted Rosseland mean opacities, κ_R, in order to be as self-consistent as possible.

We used the κ_R = 1.0 cm² g⁻¹ simulation from Schulik et al. (2020), which is similar in structure to the outcome of previous work (Ayliffe & Bate 2009b; Szulágyi & Mordasini 2017; Lambrechts et al. 2019), to introduce key concepts and explain when CPD contributions to the spectrum are important. However, as our results in our previous work indicated that mainly κ_R = 0.01 cm² g⁻¹ leads to the efficient formation of a CPD, spectral data based on this simulation is used for the remainder of the work that we present here.

Concerning the detectability of continuum features, we note how the growth from 1 µm to 10 µm in particle size already flattens the particle spectra significantly (see Fig 1). Water nominally should still exhibit features, but those are obscured by the grey carbonates. Carbonates are furthermore important for the extinction in our models, as was already noted in Ossenkopf & Henning (1994). The main remaining feature that models with 10 µm sized particles exhibit is the olivine emission feature. For millimetre-sized particles, the NIR and MIR regions are featureless but clearly broader than a single black-body.

Our modelling generally shows that the spectral continuum background in the wavelength range of 30–500 µm is a good hunting ground for CPD contributions. Further details would need to be elucidated on a case-by-case basis. The NIR shows generally planet-dominated emissions, tentatively validating the approach by Haffert et al. (2019), but we caution and note that, due to our limited gravitational potential in the hydro-simulations, we are likely overestimating the planetary contribution. Therefore, the L band could in reality be significantly more dominated by the CPD signal, which was also the conclusion arrived at in Szulágyi et al. (2019).

Our key finding is that the greyness of millimetre-sized particles translates into a black-body-like FIR slope of λF_λ ~ λ⁻³. While one could expect to find shallower slopes from 2D, optically thin disc models (Backman et al. 1992) or 3D effects, those do not show up in our simulations. This is to be contrasted with the FIR slopes of 1-10 µm particles, which should give λF_λ ~ λ⁻⁵−λ⁻⁴. Based on the preceding discussion of low opacity for small particles in the FIR due to the scattering slope, the λF_λ ~ λ⁻³ slope should be detectable even in the presence of smaller particles according to an MRN distribution, and should thus be a clear signal of dust growth.

It is well known, however, that realistic icy coatings and crystalline structures can reduce the opacity slope to κ_R ~ λ⁻¹ (Ossenkopf & Henning 1994; Jager et al. 1998). Additionally, fluffy growth and warming effects from 20 to 60 K change the FIR slopes of olivines through the activation of phonon modes (Henning & Mutschke 1997; Mennella et al. 1998; Demyk et al. 2017). It is unlikely, however, that all those effects can conspire to produce an entirely grey opacity function in the FIR. Studies like those performed for pre-stellar cores (e.g. Chacón-Tanarro et al. 2019), in which a change in the submillimetre extinction slope can be detected in one object, should be performed with ALMA and would be able to give the final word on whether grey – and hence millimetre-sized – particles exist in objects like PDS 70c.

Dust continuum modelling is a particularly strong tool for getting the overall spectral shapes correct, as opposed to, for example, the modelling of individual gas lines. The infrared data from Haffert et al. (2019) and the submillimetre data from Isella et al. (2019) allowed us to attempt exactly this. The late-type models that were designed to fit the derived dust masses in Isella et al. (2019) move the CPD spectrum into a reasonable range between the data points. As we used a factor-of-three decrease in the effective dust-to-gas ratios to match the overall optical depth of the hydro runs with the dust populations, this means that our preferred dust mass for the PDS 70c disc is one third less than the previously reported value. With the modelling performed, it seems that only using millimetre-sized particles results in a spectrum that is broad enough to fit this data and consistent with a planetary body and at a similar mass. The residual planetary black-body fit resulting from this is excessively large (R_P ≈ (4.5 ± 0.5)R_Jup), as this is an unrealistic planet size Mordasini et al. (2017) at the PDS 70 age of 3–5 Myr (Keppler et al. 2018)). We note, however, that increasing the dust mass arbitrarily, we can also fit PDS 70c with smaller dust sizes, but this requires the CPD to be optically thick.

This can be taken as an indicator of ongoing dust growth in PDS 70c for the following reason. The final width of our CPD spectra will be less wide at the NIR end than the spectra presented here, due to our overestimation of the planetary contribution (the simulated contribution, not the black-body fit of the data). Therefore, more realistic CPD simulations will require grown particles to fit both ends of the spectrum simultaneously. While our overall fit seems to indicate that the dust mass estimate for PDS 70c is correct, we caution that the late-type simulations are not fully self-consistent simulations. Realistic CPDs of this mass might be hotter, but unfortunately it is unfeasible at the moment to construct such CPDs in simulations through the successive accretion of gas, due to prohibitively long simulation timescales. Therefore, as the radiation hydrodynamic models need some improvement, we need more data points in order to ascertain the uniqueness of the fit.

The interpretation of what real age the late runs correspond to for an object like PDS 70c is unclear. While in our simulations the CPDs have only just formed and started to accrete, a PDS 70c-like mass would be in a late stage of evolution, and in the interpretation of Isella et al. (2019), a PDS 70c-like mass would be late because they expected more massive objects and interpret the majority of the mass as having already accreted into the planet. This ambiguity is obviously not settled and awaits future work.

Beyond this, Shibaike & Mordasini (2024) recently inferred in 1D in situ dust growth models that dust growth to 1 mm in size is plausible in the dense CPD environment, and Karlin et al. (2023) showed via multi-fluid dusty 3D simulations that direct delivery of dust of up to 100 µm in size to the CPD is realistic. While the latter simulations, similarly to our earlier work in Schulik et al. (2020), were limited in their radial extent surrounding the planet due to the 3D approach, the 2.5D radiation hydrodynamic simulations of Marleau et al. (2023) resolved the planet. They show that the circumstellar midplane gas flow, which should carry the largest particles, can be distributed over large radii in the circumplanetary midplane. While all those results were obtained in different set-ups and under different assumptions, and a fully self-consistent 3D radiation hydrodynamic multi-fluid dust-growth solution of the problem is still lacking, a plausible synthesis seems to be that millimetresized dust can either grow or be delivered to the CPD. The remaining outstanding question is whether the small dust can be removed quickly enough, either through growth or drift, to allow for the spectral distribution to be consistent with the mean optical depths assumed throughout this work. Nevertheless, our study shows that the spectral width from 3µm-1mm for an object like the protoplanetary disc surrounding PDS 70c can be used to infer active dust growth, which is a necessary step in the formation of satellites.

Acknowledgements

MS wants to thank Sebastiaan Haffert and Per Bjerkeli for discussing flux conversion at various wavelengths and Sebastian Lorek for extensive discussion on dust populations in the solar system. MS was supported by a project grant from the Swedish Research Council (grant number 2014-5775). AJ wants to thank the support by the Knut and Alice Wallenberg Foundation (grant number 2017.0287), the European Research Council (ERC Consolidator Grant 724687-PLANETESYS) and VR grant 2018-04867. ML thanks the Knut and Alice Wallenberg Foundation under grant 2017.0287. All the simulations presented in this work were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lunarc in Lund University, Sweden, and the entire team is grateful for being supported with their expertise.

Appendix A Conversion of observational data

The comparison of our simulated spectral data with the real observations of PDS 70c requires some conversion operations, which are described in the following. Radmc3d puts out spectral data F_v in erg s⁻¹ cm⁻² Hz⁻¹ at 1 pc distance. The conversion of those fluxes to the universal quantity νF_v = λF_λ is trivial via F_λ = F_v × λ² /c, and then rescaled with the distance of PDS70 (113 pc).

The infrared data is given as magnitude differences, Δm_B, for band B with regard to the host star (Keppler et al. 2018; Müller et al. 2018; Haffert et al. 2019) where the latter is the magnitude m_PDS70,B , and the stellar magnitudes are listed in Keppler et al. (2018). Conversion to fluxes is then performed via the fundamental relation $$F_{\text{B}}=F_{0,\text{\hspace{0.17em}B}}\times \overline{\lambda }\times {10}^{-\left(\delta m_{\text{B}}+m_{\text{PDS}70,\text{\hspace{0.17em}B}}\right)/2.5}$$(A.1)

where $\overline{\lambda }$ is the band-midpoint, and F_0,B are the standardized flux-zeropoints for Vega, and therefore we use F_0,Ks = 4.3 × 10⁻⁷erg/cm²/µm, F_0,L′ = 5.26 × 10⁻⁸erg/cm²/µm. While those are not the most precise values (e.g. for NACO/SPHERE a more involved calibration process is described in Vigan et al. (2016)), they reproduce the values in Müller et al. (2018) for PDS 70b, and the aim of our fit is the broadness of the spectrum. Then, with m_Ks = 7.75, m_L′ = 8.84, δm_Ks = 8.8, δm_L′ = 6.6, and the previous specifications, we get λF_λ(2.18µm) = 1.1 × 10⁻¹³erg s⁻¹ cm⁻², λF_λ(3.8µm) = 1.1 × 10⁻¹³erg s⁻¹ cm⁻².

The sub-millimetre data from Isella et al. (2019) is taken as 0.1mJy/beam, which is mentioned as the peak value, and we estimate 1 beamwidth as setting most of the flux. Hence $v{F}_v(855\mu m)=\frac{c}{855\mu m}{10}^{-27}{\text{erg\hspace{0.17em}s}}^{-1}{\text{\hspace{0.17em}cm}}^{-2}{\text{\hspace{0.17em}Hz}}^{-1}=3.5\times {10}^{-16}{\text{erg\hspace{0.17em}s}}^{-1}{\text{\hspace{0.17em}cm}}^{-2}.$. The conversion factor from erg s⁻¹ cm⁻² to W m⁻² is 10⁻³ .

Error bars for the Ks, Lp bands are taken from the literature, the error bar for the submillimetre flux has been estimated the 2- sigma width of the data, assuming the CPD signal is a Gaussian.

References

All Tables

All Figures

Fig. 1 Computed opacity functions for water ice, olivine and organics for individual components of single particle sizes (Left), scaled according to the available dust mass according to Eq. (2). The millimetre-sized particles make up the ‘L’ particle mix, whereas the 10 µm particles make up the ‘S’ particle mix (see also in the main text and Table 1). In general, it becomes evident that 10 µm particles are already fairly grey, except for the water ice. Furthermore, we show how this greyness, and particularly the disappearance of the prominent 10 µm feature for olivine, results from step-wise dust growth (Right). The opacity functions for the solid lines on the right are computed for pure olivine at single sizes only. The step-wise increase in size highlights the disappearance of the 10 µm-SiO feature in our spectra for dust sizes ≥10 µm, as then the size parameter 2πa/λ enters the interference regime of Mie-scattering.

In the text

Fig. 2 Cylindrically averaged radial-vertical map of κR computed only for ice and olivine components, centred on the planet for sim 1. This highlights the approximate 3D structure of ice and olivine sublimation ellipsoids in the envelope, as is seen via the corresponding opacity drops. Overlying black contours denote the gas density in g cm−3 to highlight the disc-structure of the hydrodynamic solutions. The Hill sphere is marked with a red circle. The actual τ = 1 surfaces at λ = (5, 30 , 300) µm computed with the full solution for sim 1 and viewing from z = +∞ are overlaid in dark blue, light blue, and cyan. Dashed colours denote the mirrored continuation of the surfaces below the disc midplane. The τ = 1 surfaces help to interpret spectral information.

In the text

Fig. 3 Same as Fig. 2 but for the low-opacity hydro-simulations and sim 6, which has the same optical depth as sim 1 after a ×100 density increase. Again, the background colour is the Rosseland opacity computed only with olivine and water ice to highlight compositional differences. The τ = 1 emission surfaces were again computed with the full opacity model including all species, as was used in the paper to compute spectra.

In the text

Fig. 4 Intensity map of a synthetic observation of the planet and CPD at 5 µm. The image is as is to be expected from the computed 5 µm-emission surface in Fig. 2. The CPD is not visible, as at those wavelengths water ice and organics as opacity sources have sublimated away. This image is also representative of any simulation set-up where we find planet-dominated spectral regions, which can also happen up to the sub-millimetre. It is then usually the black-body-tail from the planetary contribution overshining the CPD signal. The latter is then weak, as the flux is Fλ ∝ κ(λ)Σ, with Σ being the surface density, and κ(λ) ∝ λ−2 in the far-infrared. Increasing the contrast to arbitrary, unrealistic levels would still make the CPD contribution visible (see text for more details).

In the text

Fig. 5 Synthetic intensity maps for the planet and CPD for sim 1 at 30 µm seen at 50° inclination, to mimic PDS 70c. Those maps assume a spatial resolution far superior to ALMA, but serve the purpose to illustrate the intensity contributions to the spectral flux at this wavelength. On the left we show a intuitive heat map, on the right with a more contrasted colour map. On the right, one can distinguish several regions according to their emission intensity I : Emission from the innermost planet, originating from olivine, I ≈ 1, the CPD with the water ice line I ≈ 0.1–0.6, also marked with the white arrows, the outer edge of the CPD I ≈ 0.03, various disc structures including spiral arms at I ≈ 0.01. The arrow marks the water ice sublimation line, which through its associated drop in opacity reveals some inner parts of the envelope.

In the text

Fig. 6 Detailed impact of including an increasing number of opacity carriers and radial integration regions onto the optically thick surface and the spectrum. Left and middle: effects of successively adding opacity carriers into sim 1. We show the emission surfaces in real space at 30 µm on the left, again integrated from z = +∞, with the radial distance on the x axis and the vertical distance on the z axis. The resulting full spectrum is shown on the right. The olivine (Ol) sets the dominant part of the spectrum, with the contribution from the water ice being minuscule mainly due to being outshone by the olivine close to the planet. Once the organics components are added (Cmajor and Cminor), the central emission surface is lifted up to colder temperatures (longer wavelengths), leading to a weaker signal in the NIR through extinction. Interestingly, the carbon greys out its own sublimation signature in the emission surface, whereas the sublimation line is clearly visible for water in the emission surface, but only weakly in the spectrum. When measuring the CPD contribution to the spectrum, we subtract everything outside of 0.1 rH from the simulation domain. This process can be seen on the right side, when following the dashed blue line as it changes into the dotted blue and then the solid blue line. It becomes evident that this spectrum, when all dust ingredients are included, gains a lot of flux for λ > 10 µm, while its gain is minor for λ < 10 µm. This leads to the conclusion that for λ > 10 µm this model is CPD-dominated. The same can be repeated for the olivine-only spectra. Those inform us that already the hot olivine blob seen on the left can be counted as CPD contribution. Right: same as on the left side, but for sim 6, deriving from the low κR-data, but having the same optical depth as sim 1 through a ×100 density increase. Water emission from the CPD plays a larger role in this spectrum, particularly between 30 µm and 100 µm.

In the text

Fig. 7 Comparison of the spectral shape in the two main hydro simulations and studying the effects of dust growth. Left: spectral energy distributions for sims 1 and 2 with all dust species included and normalized to the distance of PDS 70, 113 pc. Furthermore, we plot the fluxes for PDS 70c in the Ks (2.18 μm), L′ (3.8 μm) and 855 μm bands reported by Haffert et al. (2019); Isella et al. (2019). For conversion details and flux errors, see Appendix A. The observational data as well as our model spectra are compared in width to black-bodies (BB), plotted as thin, dashed curves. There are BB (1) (T = 300 K), BB (2) (T = 850 K) and BB (3) (T = 1000 K). PDS 70c fails to be modelled by BBs (1) or (2) and due to this width, it reveals to possess a CPD component. Our model data compared with BB (3) shows it to have multiple BB contributions. Grey shading indicates a weak CPD contribution, i.e. FCPD > 0.1FPlanet and thick lines indicate a dominant CPD contribution, i.e. FCPD > FPlanet- Right: spectra for the low κR-simulations, emphasizing the effect of dust growth on the spectra. To keep the optical depth constant for the smaller dust particles with higher opacities, we adjusted fDG downwards accordingly. For all wavelengths, the CPD is transparent and the emission surfaces are convex features that are situated around the planet, similar to the olivine-only case in Fig. 6 (left). The comparison of spectral width to the black-body (dotted black line) shows that for the millimetre-particle case, the spectrum is a convolution of black-bodies over the 3D structure, while the origin of the width is less clear for the smaller particles. Furthermore, the growth from ISM-like towards millimetre-sized particles 1) greys the spectrum and 2) produces a FIR-tail that is significantly shallower than the expectation from single-temperature, single-small-particles, which would be a ~λ−5-slope.

In the text

Fig. 8 Spectral shapes of simulations with nominal and increased dust masses. Left: late models with a ×100 times increased dust mass show significantly increased CPD fractions in their spectra (solid lines), due to the formation of emission surfaces in the CPD. We keep the previously introduced black-bodies for very rough expectations on the planet and CPD emissions. FIR slopes are again black-body-like for millimetre-sized particles, while the overall steepness is not as dramatic as for the early-type small particles, making a fit more difficult. However, the width of the PDS 70c data can be approximately fitted when considering larger grains. The artificially cooled down envelopes (dashed lines) are too cool to fit the observations and particularly undershoot the ALMA measurement. Right: degeneracy check for the fits on the left. We attempted verify whether increases in dust mass of smaller particles can also fit the seen width, which we can confirm here. Hence more data is needed to distinguish between mass and particle size. Curves are plotted without the indicators for the CPD fraction, as the viable fits all are optically thick in the CPD region and hence the sub-millimetre signal is actually emitted by the CPD.

In the text