Thermal Structure and Millimeter Emission from a Protoplanetary Disk with Embedded Protoplanets from Radiative Transfer Modeling

The discovery of protoplanets and circumplanetary disks provides a unique opportunity to characterize planet formation through observations. Massive protoplanets shape the physical and chemical structure of their host circumstellar disk by accretion, localized emission, and disk depletion. In this work, we study the thermal changes induced within the disk by protoplanet accretion and synthetic predictions through hydrodynamical simulations with postprocessed radiative transfer with an emphasis on radio millimeter emission. We explored distinct growth conditions and varied both planetary accretion rates and the local dust-to-gas mass ratios for a protoplanet at 1200 K. The radiative transfer models show that beyond the effect of disk gaps, in most cases, the circumplanetary disk (CPD) and the planet’s emission locally increase the disk temperature. Moreover, depending on the local dust-to-gas depletion and accretion rate, the presence of the CPD may have detectable signatures in millimeter emission. It also has the power to generate azimuthal asymmetries that are important for continuum subtraction. Thus, if other means of detection of protoplanets are proven, the lack of corresponding evidence at other wavelengths can set limits on their growth timescales through a combined analysis of the local dust-to-gas ratio and the accretion rate.


Introduction
The observational detection of protoplanets through different means has become a reality in the present decade.A handful of protoplanets have been confirmed or nominated through direct imaging methods, including PDS 70 b and c (Keppler et al. 2018;Haffert et al. 2019), AB Aur b (Currie et al. 2022), and HD 169142 b (Gratton et al. 2019;Hammond et al. 2023).Additionally, a larger number of protoplanet candidates have been proposed through velocity "kinks" and footprints in disk gas kinematics (Pinte et al. 2019(Pinte et al. , 2020) ) and more recently via detections of circumplanetary disks (CPDs) in gas or dust emission (Benisty et al. 2021;Bae et al. 2022).While existing facilities (e.g., Atacama Large Millimeter/submillimeter Array/ALMA, the James Webb Space Telescope/JWST, and ground-based optical/infrared telescopes) will continue to push this field, the upcoming large ground-based optical observatories, such as the European Extremely Large Telescope, the Thirty Meter Telescope, and the Giant Magellan Telescope, will further revolutionize our knowledge of planet formation.This increased observational capacity can potentially help to achieve characterization of young planet-forming systems at multiple wavelengths.This means we will be able to observe and analyze (proto)planets and their natal disks across a wide range of wavelengths, ranging from optical light to infrared.
During their formation, giant planets change the physical and chemical conditions in their surroundings by carving gaps in the gas and dust, which leads to radial differences in the temperature and radiation fields (Facchini et al. 2018;van der Marel et al. 2018;Alarcón et al. 2022;Broome et al. 2023).They cause chemical exchange between layers through distinct meridional flows of gas and dust (Szulágyi et al. 2014(Szulágyi et al. , 2022;;Dong et al. 2019) and locally elevate the temperature in the surrounding disk, potentially inducing changes in the chemical composition (Cleeves et al. 2015).Much of the extant work in the references above is based on ALMA images of gas and dust in disk combined with theoretical models.However, additional observational constraints are becoming available.As an example, the detection of UV continuum toward PDS 70 b supplies information about its accretion rate and possibly about its CPD (Zhou et al. 2022).From a chemical point of view, determining the effective accretion rate is important for determining the UV and Hα luminosity, which sets up the thermal and photochemical equilibrium of the planet-feeding gas, where atomic tracers may be dominant (Alarcón et al. 2022).
In this investigation, we explore the thermal effects of a growing planet and the expected emission at millimeter wavelengths by postprocessing models with different accretion rates and local levels of dust depletion both inside and outside the millimeter/pebble disk.We then use these predictions to study the effective protoplanet footprint within the current capabilities of ALMA.
This paper is organized as follows: we present the methodology applied in the work to study the effect of the radiative transfer of a planetary source in a circumstellar disk in Section 2. In Section 3 we present the main results of our simulations, exploring different parameters relevant in planetforming settings.We discuss the main findings of our work in Section 4 and we summarize the present paper in Section 5.

Methods
We explore the thermal differences associated with the emission from a protoplanet within a protoplanetary disk.We first run 3D hydrodynamical simulations with FARGO3D (Benítez-Llambay & Masset 2016), and then postprocess these simulations with RADMC-3D (Dullemond et al. 2012) by including the blackbody emission from a protoplanet, accretion shocks, and the spectrum of a circumplanetary disk.We then produced synthetic continuum emission images of the hosting disk at λ = 1.3 mm observable by ALMA at typical angular resolutions of dust continuum emission, ∼3 au (Andrews 2020).

Hydrodynamical Simulation Setup
We ran two 3D hydrodynamical simulations with FARGO3D (Benítez-Llambay & Masset 2016), which is built by expanding the FARGO algorithm (Masset 2000), which is particularly efficient for advection flows in Keplerian disks.Using a dedicated hydrodynamics code provides a more realistic physical setup for the radiative transfer models.One of them replicated the best-fit 2D model of the dust continuum emission from HD 163296 from the DSHARP survey (Andrews et al. 2018;Zhang et al. 2018).We also embedded a planet beyond the edge of the observable millimeter disk for the other run.We decided to replicate the HD 163296 disk as our starting point because it presents a well-known and thoroughly characterized protoplanetary disk, particularly at high resolution in millimeter wavelengths as part of both the DSHARP and MAPS surveys (Andrews et al. 2018;Öberg et al. 2021).The HD 163296 disk presents a rich substructure in dust emission with crescents, dark dust gaps, and bright emission rings (Isella et al. 2018).Given this structure, it has been proposed in the literature that the HD 163296 disk hosts several protoplanets, which are predicted via simulations that match the structures seen in dust continuum emission (Isella et al. 2018;Zhang et al. 2018;Garrido-Deutelmoser et al. 2023).Additional evidence for protoplanets is found in the gas kinematics acting as planetary footprints (Pinte et al. 2019(Pinte et al. , 2020;;Teague et al. 2019;Izquierdo et al. 2022Izquierdo et al. , 2023) ) and from the emission from atomic tracers (Alarcón et al. 2022).In particular, for our hydrodynamical run, we located three planets of 0.7, 2.18, and 0.14 M Jup at 10, 48, and 86 au respectively, although we focus on the effect of the planet at 48 au, which is the most massive inside the millimeter disk.
In our second modeling setup, we analyze the effects of a circumstellar disk with a massive planet located beyond the edge of the millimeter dust disk, i.e., where dust emission falls below a given sensitivity limit.In this model, we embedded a 2 M Jup planet at 150 au beyond the observable millimeter disk.It is expected that the effects, in this case, will be more noticeable and easier to decouple from the stellar-disk interactions as the disk becomes colder in the outer regions.A mass tapering for planetary growth in both models was included over the first 100 orbits.
Our models assume a constant α-viscosity prescription (Shakura & Sunyaev 1973) across the disk with α = 10 −3 .The 3D grids for the two hydrodynamical simulations have (256, 512, 64) cells in spherical coordinates (r, f, θ).The colatitude f spans from 1.17 to π/2, assuming mirror symmetry with respect to the midplane, while the azimuthal coordinate covers the whole range from 0 to 2π.The radial coordinate ranged from 7 to 200 au for the HD 163296 disk model and from 20 to 400 au for the single-planet model with a logarithmic spacing.The radial velocity has antisymmetric boundary conditions at the radial edges and symmetric at both ends of the colatitude range, while the colatitude velocity is symmetric in radius and antisymmetric at the colatitude boundaries.The hydrodynamical simulation ran for 500 orbits of the most massive planet at 48 au in the HD 163296 model and for the same number of orbits for the model with one planet at 150 au.An aspect ratio of 0.08 with a flaring index ψ = 0.08 was considered in the initial density distribution of the disk.The planets were input in fixed circular orbits using the planet-to-star mass ratio, q, assuming a stellar mass of 2.1 M e .We show the dust density fields of our 3D hydrodynamical simulation of HD 163296 in Figure 1.In the simulation, we see the density depletion in the gaps as expected at the location of the massive planets along with the planet-driven spirals.

Radiative Transfer Setup
We studied the effects of the emission from a protoplanet and its CPD by postprocessing the output from the hydrodynamical simulations of FARGO3D with RADMC-3D (Dullemond et al. 2012;Benítez-Llambay & Masset 2016).Within the heating components in the modeling, we took into account the viscous heating associated with the stellar accretion rate, assuming a stellar accretion rate of M M 0.69 10 yr ´-- (Donehew & Brittain 2011).The viscous heating was modeled by adding concentric rings following the standard temperature distribution of an accretion disk (Calvet et al. 1994): We used 10 8 photons to calculate the thermal equilibrium.The stellar spectrum is the one corresponding to the HD 163296 star given that we set up the hydrodynamical simulations based on that disk.The input spectra for the star and the protoplanet are shown in Figure 2.
We modeled the dust opacity in the disk using two dust populations for the radiative transfer modeling.The first population corresponds to large grains holding 99% of the disk dust mass, and the other population represents only small grains with 1% of the disk dust mass.The two populations have the same composition, which corresponds to a water ice mass fraction of 20%, 40% in organics, 33% astrosilicates, and ∼7.8% in troilites.The opacity values were taken from Warren & Brandt (2008) for water ice, Henning & Stognienko (1996) for organics and troilite, and Draine (2003) for the astrosilicates.
We set a constant dust-to-gas mass ratio across the disk of 0.01 (Sandstrom et al. 2013) with the exception of the Hill sphere around the planet at 48 au, which is a free parameter in our modeling that will be detailed below.Zhu et al. (2012) demonstrate that planet-induced gaps filter the drifting of large grains; additionally, Zhu et al. (2018) show that the drifting timescale of millimeter and submillimeter dust grains in massive CPDs can be as short as t drift ∼ 1000 yr.Thus, we explored different scenarios, varying the local depletion of large grains, which is reflected in the local dust-to-gas mass ratio, assuming values of 10 −2 , 10 −3 , and 10 −4 .As part of the radiative transfer simulation with RADMC-3D, anisotropic scattering with the Henyey-Greenstein phase function was considered (Henyey & Greenstein 1941).The dust size distribution for the two populations follows the Mathis et al. (1977) power law proportional to the dust grain size, i.e., with a the grain size, and p the power-law index, which we set to 3.5.The first dust population ranges from a min = 0.05 μm to a max = 2.5 mm, while the second spans sizes between a min = 0.05 μm and a max = 5 μm.

Protoplanetary Spectra
We added an additional point source representing the emission from the protoplanet and its associated CPD.The spectrum of the planet is a combination of three components: 1. Thermal emission coming from the protoplanet itself, which is modeled as a blackbody emitting with a temperature of T p = 1200 K. 2. The emission coming from the accretion of the CPD using the viscous heating from Equation (1) inside the Hill sphere of the planet.For this case, viscous heating will be dependent on the mass of the planet, M P , the planetary accretion rate, M P  , treated as a free parameter, the radius of the planet, R P , and the separation between the star and the planet, R. 3. The shock emission is estimated following the prescription of Calvet & Gullbring (1998) and Zhu (2015).The temperature of the shocked gas, T sh , is where m H , μ, and k are the hydrogen mass, the mean molecular weight, and the Boltzmann constant respectively; v s is the freefall velocity at the planetary surface: with G the gravitational constant, R P the planetary radius, and R i the inner radius of the accretion disk.
We added an extra emitting point source in the radiative transfer modeling to explore the effects of thermal emission from a massive protoplanet, the CPD, and the emission of shocks.The additional point source corresponds to the most massive planet at the same location as the hydrodynamical simulation, i.e., at 48 au for the HD 163296 model and at 150 au for the single-planet one.The spectral energy distribution of the planets is dependent on multiple factors such as the accretion rate, M  , and the assumed effective shock area over the surface of the planet.Hence, we explored three different values for the accretion rate: M  =10 −7 , 10 −6 , and 10 −5 M e yr −1 .
To account for the drift of large dust (Zhu 2015) we also change the local dust-to-gas ratio inside the typical size of a CPD, which is around one-third of the planet's sphere of influence, i.e., the Hill sphere (Ayliffe & Bate 2009;Szulágyi et al. 2014;Perez et al. 2015).The Hill radius, r H , in the circular orbits of the embedded planets in our simulations is characterized by where a is the orbital radius of the planet, M P the planetary mass, and M * the stellar mass, which in our case matches the mass of HD 163296, i.e., M * = 2.1 M e (Öberg et al. 2021).Given our assumption of a planet of mass 2 M Jup at 48 au and at 150 au, the Hill sphere has a size of 3.27 au and 10.2 au for each model respectively.We used three different values for the dust-to-gas ratio inside the CPD: an interstellar medium value of d2g = 10 −2 , a moderate depletion value d2g = 10 −3 , and a severely dust-depleted case with d2g = 10 −4 .

Thermal Equilibrium and Synthetic Images
We postprocessed the hydrodynamical simulations by creating synthetic images using RADMC-3D (Dullemond et al. 2012) to study the observability of the distinct accretion and CPD setup at radio wavelengths.We created 1.3 mm dust continuum images for ALMA assuming a disk position angle of 43°.4 and an inclination of 46°.6.These angles follow the respective fitted values from observations of high-resolution gas kinematics in the HD 163296 disk (Teague et al. 2021).We used the SIMIO package1 to reproduce the observing conditions of the HD 163296 disk, matching the antenna configuration, the coverage of the uv space, and the time integration.
It is expected that, due to the lower optical depth, the temperature inside dust-depleted gaps is higher than that in a smooth disk (Facchini et al. 2018;van der Marel et al. 2018;Alarcón et al. 2020), which has been observationally confirmed to be true for HD 163296 (Calahan et al. 2021).Thus, to isolate and avoid mixing the effect of the gap itself, we compare the thermal fields with the hydrodynamical simulation with the planet embedded, turning on and off the CPD and planetary emission.We apply the same criterion for the synthetic ALMA continuum images at 1.3 mm where those effects have already been considered.
The resulting synthetic ALMA images were then subtracted and compared with a postprocessed hydrodynamical simulation with planet-driven dust-depleted gaps but without the presence of a secondary emitting source.

Temperature Structure
We recover the expected thermal changes of previous studies by postprocessing hydrodynamical simulations with radiative transfer, i.e., the temperature increases around the planet's location by tens of kelvin (Cleeves et al. 2015;Montesinos et al. 2015;Oberg et al. 2022;Portilla-Revelo et al. 2022).Our work expands on those studies by splitting the additional energy source between the blackbody emission from the planet, accretion shocks, and the CPD itself (see Figure 2).
In Figure 3 we illustrate the thermal fields in a meridional cut (r, z) at the planet's location and the (x, y) thermal field in the midplane from the models with the planet's temperature T P = 1200 K and a high planetary accretion rate, M 10 5  = -M Jup yr −1 .The point-source emission from the planet effectively changes the physical temperature in the surroundings on physical scales up to ∼5 au at 48 au, i.e., at least one disk scale height.The extent of this thermal effect is comparable to the size of the Hill sphere.
We show a comparison of the models with the extra emitting sources at 48 au for the different setups in Figures 4 and 5. From the thermal distributions, we observe that in all cases there is a noticeable effect in azimuth and height for accretion rates M 10 7  > -M e yr −1 , indistinctive of the local large dust depletion.The expected trend in temperature with the planetary accretion rate is confirmed.As the planetary accretion rate increases, the CPD and shock emission do so as well, producing an extra local heating effect.For M 10 7  = -M e yr −1 , the effects are almost negligible with very small changes in the vicinity of the planet.It is noticeable that very high accretion rates M 10 6  > -M e yr −1 produce strong thermal signatures, i.e., the temperature increases by at least 20 K.Such thermal increases are significant enough to change the chemical balance in the gaseous and solid composition.When looking at the trend of different local dust-to-gas ratios, the effect is slightly more significant for the models with severe dust depletion values.As dust depletion increases, the material surrounding the planet becomes more transparent, increasing the mean free path of photons emitted by the planet or the CPD and making the thermal gradient flatter, although the overall extent of the thermal changes does not change drastically.However, the thermal changes within the Hill sphere are more uniform for the lowest dust-to-gas ratios.It is worth mentioning that the same increased transparency will cause an enhanced local UV flux with a considerable extent.This topic and a more careful UV radiative transfer will be the subject of a subsequent paper.The thermal increases in the disk may have a significant effect that can create a localized sublimation effect, further creating azimuthal asymmetries in the disk.However, such thermal changes are not critical when considered together with the temperature gains due to the gap itself (van der Marel et al. 2018; Alarcón et al. 2020).Figures 4 and 5 also show that the extent of the thermal differences is dependent on the parameters of the CPD, and there is an overall thermal increase in the disk due to the protoplanetary emission; however, this net thermal increase should not be significant enough because it is of the order of 1-2 K in regions far from the planet.When the temperature field is compared between azimuthally opposed extremes of the disk, the difference is only noticed as a bright hot spot surrounding the planet.Therefore, the overall 1-2 K increase across the disks within 20-70 au from the star is coming from the blackbody emission from the protoplanet.This effect is not really detectable with current state-of-the-art observatories, but it points to an overall heating of the disk just produced by the localized blackbody emission from the planet.
In Figures 6 and 7, we show the same thermal differences but compare the single-planet model at 150 au instead of the HD 163296 disk.We observe that many of the trends observed in the simulation with the massive planet inside the pebble disk persist.One of those trends is that by boosting the planetary accretion rate, the local temperature rises as the overall energy input from the shocks and CPD increases.However, there are also important distinctions related to the extent of these changes and their absolute values.
One of the main differences is the trend with the local dustto-gas ratio.With a similar planetary accretion rate, local dust depletion has the potential to produce drastic changes in the disk temperature around the planet that is found beyond the edge of the dust disk.Such changes can be extended to more than 10 au with temperature increases over 25 K. Thus, these changes are within achievable spatial scales resolvable with ALMA.Additionally, in the presence of a massive gaseous CPD, or detectable line emission, molecular tracers could probe changes in the local line ratios or in the abundance of volatile species.
When the planet is beyond the pebble disk, the overall heating of 1-2 K over the surrounding material disappears.The probable reason behind this phenomenon is the lack of processes scattering the blackbody emission due to a much lower dust density, i.e., the mean free path of dust scattering is large enough for the thermally emitted photons to escape the disk.

Millimeter Synthetic Emission
In our synthetic predictions of the 1.3 mm dust continuum emission in Figure 8, under the right combination of parameters, the CPD will not show a significant difference in yr −1 the thermal increase is noticeable, while for accretion rates  M M 10 7 Jup  yr −1 it is probable that the effect is milder and not significant, even taking into account the blackbody emission from the protoplanet at 1200 K.
emission with respect to the surrounding material, i.e., the local dust depletion counteracts the viscous heating coming from the accretion.Only for the strong accretion scenario without dust depletion in the CPD is there a strong local enhancement of millimeter emission at the planet's location.Nevertheless, as we show in the azimuthal cuts of Figure 9, in the cases of low accretion and local dust depletion, the emission resembles an azimuthal absorption feature in the millimeter dust continuum at the location of the planet.However, given the current observational capabilities, those changes can be distinguishable from the noise in the observations for deep integrations.With the exception of a high accretion rate, M  > 10 −5 M Jup yr −1 , the emission enhancements are not above the 3% level.In contrast, when we assume a dust depletion in the CPD there is a significant decrement in the local emission to a level of 10%.
The same exercise for the single planet at 150 au under ideal conditions, i.e., without observational noise, would produce changes of ∼10%; however, the disk is optically thin at millimeter wavelengths due to being embedded beyond the edge of the pebble disk, and has very low millimeter optical depth.even though the planet and accretion irradiation can produce significant relative changes, their absolute emission falls below the sensitivity threshold of a current state-of-the-art radio facility such as ALMA.In this case, no footprint or continuum emission feature is noticeable or detectable, although it may show an important effect in line emission.
We compare the peak expected deviation in brightness temperature from our synthetic millimeter models to a planet at 50 au with the rms values reported in the literature for deep high-resolution dust continuum observations with ALMA in Figure 10 ( Pérez et al. 2020;Andrews et al. 2021;Benisty et al. 2021).The figure shows that the expected deviations are marginally located at the 3σ level of detection in some cases, even with very deep integrations.However, it is possible that some of the marginal proposed CPDs in the residual images of Andrews et al. (2021) may be tracing different regimes in accretion and dust content within CPDs.
In the case of PDS 70 c, even though models in the literature predict the planet to be a super-Jupiter, the observed brightness temperature shown in Figure 10 is still consistent with the values we predict.From millimeter emission alone, it is hard to fully disentangle the accretion rate in the CPD of PDS 70 c, but it is consistent with a dust-undepleted CPD.Nevertheless, multiwavelength observations will help to constrain the accretion rate and the CPD content itself by breaking the degeneracies between those two variables.
In the past, Wolf & D'Angelo (2005) predicted a local bright emission considering a massive planet accreting at 5 au with current ALMA capabilities.However, in the bright sources observed with ALMA, such massive CPDs have not been observed.From current deep high-resolution observations, we can discard the presence of multiple high-accretion-rate scenarios without dust depletion because their local changes  > -M Jup yr −1 the thermal increase is noticeable, while for accretion rates  M M 10 7 Jup  yr −1 it is probable that the effect is milder and not significant even taking into account the blackbody emission from the protoplanet at 1200 K.

6
The Astrophysical Journal, 967:144 (11pp), 2024 June 1 are within detectable sensitivity thresholds (Andrews et al. 2018).Even though dust substructures can have multiple causes, under the assumption that a protoplanet causes them, the low detection rate of CPDs with current facilities and surveys implies two possible scenarios: 1.The first explanation is that high-accretion-rate states are temporary and relatively short, which can be explained by planets being formed on very short timescales.Lubow & Martin (2012) predict possible timescales if planetary growth undergoes accretion outbursts similar to FU Ori objects.They proposed that massive planets can undergo these states on timescales of 10 3 -10 5 yr.As the detection of more protoplanets and characterizations of their masses and accretion rates become more ubiquitous, constraints on the timescales of their formation will be more accurate, which will lead to a better understanding of planet formation.If we take a closer look at the protoplanet's orbital radius in the HD 163296 model by making an azimuthal cut at a separation of 2 M Jup , we show that the planet creates a global emission effect at all azimuths, which follows from the global thermal increase.Nevertheless, there are still local modulations at the location of the planet that depend on the accretion rate and local dust-to-gas ratio.As mentioned above, when there is strong dust depletion within the Hill sphere of the planet, the local emission is lower, matching that from the models without the planet present, which may look like an effective "absorption" in the azimuthal profile.If we focus on the accretion rate, the viscous heating and shock emission only play a major role if the accretion rate is very high and the local dust-to-gas ratio is kept similar to the fiducial value of the disk.The lack of detections of CPDs in dust continuum images suggests that either planets do not have constant high accretion rates, i.e., they go through accretion outbursts analogous to FU Ori objects, or there is a constant replenishment of dust grains to the CPD, keeping the dust-to-gas ratio roughly constant inside it.Doing the same exercise with the planet embedded beyond the pebble disk, where the thermal changes are more significant, the relative changes in emission become more important as well.However, the low local dust surface density causes the emission to fall below the detectable sensitivity thresholds of ALMA.Thus, even though there may be a significant temperature change beyond the pebble disk, it will not be noticeable by ALMA observations of the dust continuum.Nevertheless, such changes will manifest through line ratios and thermal changes in gas emission, which are potentially detectable under the right conditions.

Observational Signatures of Protoplanets and/or CPDs
Former models explored the effect of depleted planet-carved gaps in protoplanetary disks and their observability (Rosotti et al. 2016).Previous work in the literature (Wolf & D'Angelo 2005) showed that CPDs close to the star, if present, should be easier to detect at millimeter wavelengths with the best current observatories, which has not been the case except for PDS 70 b (Isella et al. 2019;Benisty et al. 2021).Further, Szulágyi et al. (2018) tested the emission coming from planets at ALMA wavelengths, showing that massive planets have detectable footprints at short wavelengths with ALMA.Such an approach was further expanded by Binkert et al. (2021), assuming a dust component and using a range of Stokes numbers.However, the physical conditions of material surrounding protoplanets can span orders of magnitudes, broadening the range of Stokes numbers that a single particle size goes through while being accreted into the planet.Zhu et al. (2018) mentioned that CPDs may be detectable as long as they get replenished with large grains from the host disk.Andrews et al. (2021) localized negative and positive residuals that can be consistent with the presence of CPDs; however, it is still hard to differentiate between real features and some artifacts not related to CPDs.Deeper resolved observations with higher sensitivity at different epochs should be able to confirm/discard the prevalence of the structures and their possible origin as CPDs.
Other methods of finding protoplanet candidates or CPDs through disk kinematics, line emission, or small-grain content will be key in constraining the dust filtering and growth of planets in their early stages.It is expected that massive planets will produce velocity kinks in the channel maps of line emission (Pérez et al. 2015(Pérez et al. , 2018;;Pinte et al. 2019).The velocity kinks allow us to constrain the mass of a present protoplanet or at the very least to place an upper limit.The asymmetries in line emission can help to constrain the dust depletion at the protoplanetary gap due to local chemical changes at these locations (Cleeves et al. 2015).
Looking for CPD signatures beyond the pebble disk in probes of millimeter dust continuum emission is hard given the low optical depth at millimeter wavelengths.Regardless of the lack of millimeter continuum emission, the models presented above show that there are significant thermal changes produced by the protoplanet emission, the CPD, and possible accretion onto the planet.Therefore, such changes may still be traceable by other means; temperature increases have the potential to lead to thermal desorption of otherwise absent volatiles, making them easier to observe due to the weak continuum emission and isolated emission.Moreover, the local density of dust can be low enough that the millimeter emission is negligible while making the CPD brighter than its surroundings at infrared wavelengths.A more thorough analysis of this behavior will be the topic of a subsequent paper.
Beyond the edge of pebble (millimeter) disks, It is possible that the changes linked to the planet's emission may present themselves through the sublimation of key volatiles (Cleeves et al. 2015).More sensitivity and deeper studies of line emission in their optically thin regions may point to the presence of otherwise hidden CPDs.Shocked gas during the accretion of massive planets can reach very high temperatures (>3000 K), producing optically thick atomic blankets that hide emission from the protoplanet (Szulágyi & Mordasini 2017).Shocks can also play a major role in local photochemistry through photodissociation and photodesorption as well, which could be probed by the emission from other tracers such as SO or CI (Alarcón et al. 2022;Law et al. 2023).A tentative detection of a CPD beyond the pebble disk has been suggested in the AS 209 disk, where there is a local enhancement of the predicted gas density, but also a slight increase in the temperature of the emission measured by CO isotopologues (Bae et al. 2022).The measurement of this thermal increase is a key ingredient in determining the CPD mass, and therefore the properties of the protoplanet.
The observable effects in the millimeter dust continuum emission from the Jovian planet at 150 au are negligible.However, in most cases, when there is a lack of significant replenishment and/or a strong accretion rate, there will be a significant increase in gas temperature, which should make a strong imprint on the local gas emission.So, even though the dust millimeter continuum emission will not show strong observable footprints, the footprints of a planet should be traceable by analysis of molecular line emission from different transitions for the same tracers as well.
Overall, when high accretion rates and strong dust depletion are combined, important effects are expected in both physics 3 mm dust continuum images using SIMIO between the different models and a standard model without protoplanetary emission for ∼1 hr integration time.Apart from one case, there is no strong additional signature of the CPD.Furthermore, when the CPD is severely depleted of large grains, the millimeter emission from dust is lower than that from the surroundings.From the models with low accretion rate we also infer that the 1200 K blackbody protoplanetary spectrum does not produce a significant change in synthetic ALMA images.and chemistry.Those changes include multiple conditions such as high-energy radiation fields, local temperature, thermal desorption rates, line ratios, line emission, and reaction rates of the gas chemistry.

The Low Number of CPD Detections
The low number of CPDs detected to date through millimeter continuum observations can point to a lack of millimeter grains if massive planets are actually carving the dust gaps or to low accretion rates.A possible reason for the lack of CPD detections is that massive planets accrete the solid reservoir in CPDs on very short timescales, together with some level of dust replenishment or very efficient dust growth.Nevertheless, the net effect of the blackbody emission on their surroundings, in particular in dust-depleted gaps, may be leading to an overestimation of the local dust surface density.This means that the actual local depletion may be stronger, which can have strong implications for our understanding of the values of disk viscosity, disk shape, and predicted planetary masses from observations.Another possible explanation for the lack of CPD detections is that most of the observed dust substructures are not actually being carved by massive planets, but are being caused by other mechanisms and/or physical instabilities (see Bae et al. 2023 for a detailed overview of substructure-forming mechanisms).Many of these instabilities will lead to structures resembling the ones produced by massive planets, i.e., gaps, rings, and spirals.It is only through deep observation and/or multiple analytical approaches involving continuum and line emission coupled with modeling that the presence of a planet, or the absence of instabilities, can be tackled.

Hydrodynamics Simulation: Caveats
In our initial hydrodynamics simulations we assumed isothermal structures.It is important to note that the actual distribution of infalling material into the planet will depend on whether an adiabatic or isothermal equation of state is chosen, producing a more spherical or a flatter matter distribution respectively (Szulágyi et al. 2016;Fung et al. 2019;Paardekooper et al. 2023).The same changes in the assumed thermodynamics caused different densities around the planet.The isothermal solution produced lower densities in the CPD than the adiabatic regime.Thus, depending on the thermodynamic behavior of the gas around the protoplanet, the signature may change due to optical depth effects.
We also opted not to use FARGO3D multifluid module to study the dust distribution, in particular, because we focus on the regions close to the protoplanet.Multifluid modules generally rely on the Stokes number, St, although that depends on the stopping time and Keplerian velocity, and it does not trace a uniform grain size.Since the stopping time depends on the local surface density and the Keplerian velocity, the grain size can range over more than one order of magnitude for the same Stokes number.Thus, it is not a reliable tracer of the dust content in the CPD, unless a dedicated study is conducted, which goes beyond the scope of this paper.Krapp et al. (2022Krapp et al. ( , 2024b) ) (Andrews et al. 2021).The ellipse in the bottom left illustrates the beam size of the observations.Bottom: azimuthal cut inside the gap at 48 au for the different models.It is shown that the undepleted case with a high accretion rate produces a noticeable bright structure in the disk.Moreover, there is some degeneracy between the local depletion in the CPD and the accretion rate that hides the presence of the CPD at millimeter wavelengths, making it not noticeable.presented a new solver for studies of the particle size distribution in multifluid modules to study dust dynamics.

Summary
We present the thermal changes produced by including the emission from an accreting planet and a simple prescription for its associated shocks on its host protoplanetary disks for cases when the planet is embedded inside the optically thick pebble disk, and one when the planet is located in the more diffuse, less dense outer disk.
From postprocessing 3D hydrodynamics with radiative transfer, we show that the radiation from a massive planet increases further the local and global temperature beyond the effect of the dust substructures it produces.For low values of accretion rate, the changes in local millimeter emission are not noticeable.Moreover, we found that a high accretion rate combined with a high depletion of large grains can hide the emission from a hotter CPD.Hence, the brightness temperature of a CPD may appear higher/lower than that of its surroundings depending on the local level of dust depletion.
The temperature differences become more important beyond the millimeter disk as the stellar viscous heating becomes negligible and the disk gets cooler.However, the dust surface density is low enough that the emission is optically thin.Hence, the millimeter emission and its changes fall below the sensitivity threshold of current state-of-the-art facilities such as ALMA.Nevertheless, the thermal effects can potentially be detectable in gas emission or temperature tracers from protoplanetary disks.

Figure 1 .
Figure 1.Top: density field in the midplane of the disk after 500 orbits of the HD 163296 model using the FARGO3D code (Benítez-Llambay & Masset 2016).Bottom: meridional cut of the density field along the azimuth of the planet at 48 au showing the gap and density structure of the disk in the radial and vertical directions.

Figure 2 .
Figure 2. Input spectra of the model for different planetary accretion rates, which include the circumplanetary disk and the accretion shocks.The stellar spectrum of the HD 163296 star is also shown for reference and comparison with the protoplanetary spectra.

Figure 3 .
Figure 3. Dust temperature in the disk for a planet having an accretion rate of M M 10 5 Jup  = - yr −1 , showing the local increase in the temperature and its extent in the three directions where the protoplanet is located.Left: dust temperature in a meridional cut along the azimuthal location of the protoplanet.Right: dust temperature in the midplane of the disk.

Figure 4 .
Figure 4. Changes in dust temperature in the disk caused by the additional input spectra for the different conditions of the protoplanet and protoplanetary accretion rates along meridional cuts at the location of the planet.The differences show that for accretion rates M M 10 7 Jup  > - yr −1 the thermal increase is noticeable, while for accretion rates  M M 10 7

Figure 5 .
Figure 5. Changes in dust temperature in the disk caused by the additional input spectra for the different conditions of the protoplanet and protoplanetary accretion rates in the midplane of the disk.The differences show that for accretion rates M 10 7 > -M Jup yr −1 the thermal increase is noticeable, while for accretion rates  M M 10 7 2. The second scenario is that planets accrete large grains in their CPDs on relatively short timescales.Significant depletion of large grains within translates into lower emission from optically thin CPDs, which would seem like azimuthally darker spots inside dust gaps due to the lack of effective millimeter emitters.Despite the fact there is one detection of a millimeter CPD from PDS 70 c (Benisty et al. 2021), it required deep integrations to observe a massive planet in a severely depleted disk with confirmed emission from accretion tracers, so they are likely caught in their ongoing assembly.Regardless of PDS 70, in many other observations, deep campaigns for the detection of CPDs in dust continuum emission have been unsuccessful.Thus, millimeter emission from deep dust gaps, associated with massive planets, is already low, making the detection of azimuthal dips within them a hard task(Andrews et al. 2021).

Figure 6 .
Figure6.Meridional cut at the azimuth where the planet is located.The shift in trends from the embedded planet compared to when is located at 48 au is caused mainly by changes in the dust optical depth.

Figure 7 .
Figure 7. Changes in midplane temperature with the simulation of a single planet embedded at 150 au.The dust temperature changes when the planet is embedded beyond the edge of the pebble disk at 150 au.

Figure 8 .
Figure8.Comparison of synthetic 1.3 mm dust continuum images using SIMIO between the different models and a standard model without protoplanetary emission for ∼1 hr integration time.Apart from one case, there is no strong additional signature of the CPD.Furthermore, when the CPD is severely depleted of large grains, the millimeter emission from dust is lower than that from the surroundings.From the models with low accretion rate we also infer that the 1200 K blackbody protoplanetary spectrum does not produce a significant change in synthetic ALMA images.
emphasize the need for global radiation hydrodynamics to fully understand the dust distribution and dynamics in the vicinity of the planet.Krapp et al. (2024a) also

Figure 9 .
Figure 9. Top: 1.3 mm continuum emission image from the ALMA DSHARP Large program(Andrews et al. 2021).The ellipse in the bottom left illustrates the beam size of the observations.Bottom: azimuthal cut inside the gap at 48 au for the different models.It is shown that the undepleted case with a high accretion rate produces a noticeable bright structure in the disk.Moreover, there is some degeneracy between the local depletion in the CPD and the accretion rate that hides the presence of the CPD at millimeter wavelengths, making it not noticeable.

Figure 10 .
Figure 10.Expected deviations in brightness temperature in our suite of models with respect to the baseline emission compared with deep highresolution continuum data from ALMA programs.Black squares show that our models predict localized deviations with a brightness temperature lower than 2 K. Blue triangles show the 3σ levels from the Andrews et al. (2021) limits, while the red triangle is the one from Pérez et al. (2020) and the red circle is the observed emission from PDS 70 c.