# Thermal Structure and Millimeter Emission of Protoplanetary Disk with
embedded protoplanets from radiative transfer modeling

Felipe Alarcón  Department of Astronomy, University of Michigan, 323 West
Hall, 1085 S University, Ave., Ann Arbor, MI 48109, USA  Edwin A. Bergin
Department of Astronomy, University of Michigan, 323 West Hall, 1085 S
University, Ave., Ann Arbor, MI 48109, USA

###### Abstract

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
post-processed 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 CPD and the planet’s emission locally increase the
disk temperature. Moreover, depending on the local dust-to-gas depletion and
accretion rate, the CPD presence may have detectable signatures in millimeter
emission. It also has the power to generate azimuthal asymmetries 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.

planet–disk interactions — protoplanetary disks — — accretion disks –
submillimeter: planetary systems

††software: astropy (Astropy Collaboration et al., 2013, 2018, 2022), FARGO3D
(Benítez-Llambay & Masset, 2016), RADMC-3D (Dullemond et al., 2012), numpy
(Harris et al., 2020), matplotlib (Hunter, 2007), cmasher (van der Velden,
2020)

##  1 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, 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 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 ELT, 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 multiple
wavelength characterization of young planet-forming systems. 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 gas and
dust meridional flows (Szulágyi et al., 2014; Dong et al., 2019; Szulágyi et
al., 2022) 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 towards PDS 70 b supplies information about its accretion rate
and possibly about its Circumplanetary Disk or CPD (Zhou et al., 2022). From a
chemical point of view, determining the effective accretion rate is important
to determine the UV and H$\alpha$ 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 post-processing models with
different accretion rates and local levels of dust depletions 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 radiative transfer effect 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 planet-forming
settings. We discuss the main findings of our work in Section 4 and we
summarize the present paper in Section 5.

##  2 Methods

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

###  2.1 Hydrodynamical Simulation Setup

We ran two 3D hydrodynamical simulations with FARGO3D (Benítez-Llambay &
Masset, 2016), which is built expanding the FARGO algorithm (Masset, 2000),
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
HD 163296 dust continuum emission 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 as 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
that are predicted via simulations that match the structures seen in dust
continuum emission (Zhang et al., 2018; Isella 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, 2020;
Teague et al., 2019; Izquierdo et al., 2022, 2023); and from the emission of
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 MJup 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.

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 direction.

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 MJup 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 $\alpha$-viscosity prescription (Shakura &
Sunyaev, 1973) across the disk with $\alpha=10^{-3}$. The 3D grid for the two
hydrodynamical simulations have (256,512,64) cells in spherical coordinates
($r$,$\phi$,$\theta$). The colatitude $\phi$ spans from 1.17 to
$\frac{\pi}{2}$ assuming mirror symmetry with respect to the midplane, while
the azimuthal coordinate covers the whole range from 0 to 2$\pi$. 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 in the radial edges and
symmetric in both ends of the colatitude range, while the colatitude velocity
is symmetric in radius and antisymmetric in 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 $\psi=$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⊙. We show the dust density fields of our HD
163296 3D hydrodynamical simulation in Figure 1. In the simulation, we see the
density depletion in the gaps as expected in the location of the massive
planets along with the planet-driven spirals.

###  2.2 Radiative Transfer Setup

We studied the effects of the emission of a protoplanet and respective CPD by
post-processing the output from the hydrodynamical simulations of FARGO3D with
RADMC-3D (Benítez-Llambay & Masset, 2016; Dullemond et al., 2012). 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 $\dot{M}=0.69\times 10^{-7}M_{\odot}\ $yr-1 (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),

$q_{\rm visc}=\frac{3GM_{P}\dot{M}_{P}}{8\pi
R^{3}}\Big{[}1-\big{(}\frac{R_{P}}{R}\big{)}^{1/2}\Big{]}.$ (1)

We used 108 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 protoplanetar are shown in Figure 2.

Figure 2: Input spectra of the model for different planetary accretion rates,
which includes 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.

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 a 1% of the disk dust mass. Both populations have the same
composition, which corresponds to a water ice mass fraction of 20%, 40% in
organics, 33% astrosilicates, and $\sim$ 8.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 sub-millimeter dust grains in massive CPDs can be as short as
$t_{\rm drift}\sim$ 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 RADMC3D, anisotropic scattering with the Henvey-
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.,

$n_{0}(a)\propto a^{-p}da,$ (2)

with $a$ the grain size, and $p$ the power law index, which we set to 3.5. The
first dust population ranges from $a_{\rm min}$ = 0.05 $\mu$m to $a_{\rm max}$
= 2.5 mm, while the second dust population spans sizes between $a_{\rm min}$ =
0.05 $\mu$m to $a_{\rm max}$ = 5 $\mu$m.

####  2.2.1 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. 1.

Thermal emission coming from the protoplanet itself which is modeled as a
blackbody emitting with a temperature of Tp = 1200 K.

  2. 2.

The emission coming from the accretion of the CPD using the viscous heating
from Equation 1 inside the Hill’s sphere of the planet. For this case, viscous
heating will be dependent on the mass of the planet, $M_{P}$, the planetary
accretion rate, $\dot{M}_{P}$, treated as a free parameter, $R_{P}$, the
radius of the planet, and the separation between the star and the planet, $R$.

  3. 3.

The shock emission is estimated following the prescription of Calvet &
Gullbring (1998) and Zhu (2015). The temperature of the shocked gas, $T_{\rm
sh}$, is:

$T_{\rm sh}=\frac{3}{16}\frac{\mu m_{H}}{k}v_{s}^{2}=3.44\times 10^{6}\
\mathrm{K}\Big{(}\frac{M_{P}}{M_{\odot}}\Big{)}\Big{(}\frac{R_{P}}{R_{\odot}}\Big{)}^{-1},$
(3)

where $m_{H}$, $\mu$, $k$ are the hydrogen mass, the mean molecular weight,
and the Boltzmann constant respectively; $v_{s}$ is the free-fall velocity at
the planetary surface:

$v_{s}=\Big{(}\frac{2GM}{R_{P}}\Big{)}^{1/2}\Big{(}1-\frac{R_{P}}{R_{i}}\Big{)}^{1/2},$
(4)

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 a massive protoplanet thermal emission, the CPD, and
shocks emission. 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 SED of
the planets is dependent on multiple factors such as the accretion rate,
$\dot{M}$, and the assumed effective shock area over the surface of the
planet. Hence, we explored three different values for the accretion rate:
$\dot{M}$ =10-7, 10-6, and 10-5 M⊙ yr-1.

Figure 3: Dust temperature in the disk for a with a planet having an accretion
rate of $\dot{M}=10^{-5}\ M_{\mathrm{Jup}}$ yr-1 showing the local increase of
the temperature and its extension 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.

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:

$r_{H}=a\sqrt[3]{\frac{M_{P}}{3M_{*}}},$ (5)

were $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_{\odot}$ (Öberg et al., 2021). Given our assumption of
a 2 MJup mass planet at 48 and 150 au, the Hill sphere has a size of 3.27 au
and 10.2 au for each model respectively. We used different values for the
dust-to-gas ratio inside the CPD, adopting three different values; an ISM
value of d2g=10-2, a moderate depletion value d2g=10-3, and a severely dust-
depleted case with d2g=10-4.

####  2.2.2 Thermal Equilibrium and Synthetic Images

We post-processed 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 in radio wavelengths. We created 1.3
millimeter dust continuum images ALMA assuming a disk PA of 43.4 degrees and
an inclination of 46.6 degrees. These angles follow the respective fitted
values from high-resolution gas kinematics observations of the HD 163296 disk
(Teague et al., 2021). We used the SIMIO package
111https://github.com/nicokurtovic/SIMIO 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 hotter than the one of 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 hydrodynamic 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
post-processed hydrodynamic simulation with planet-driven dust-depleted gaps
but without the presence of a secondary emitting source.

##  3 Results

###  3.1 Temperature Structure

We recover the expected thermal changes of previous studies by post-processing
hydrodynamical simulations with radiative transfer, i.e., the temperature
increases around the planet’s location by tens of Kelvin degrees (Montesinos
et al., 2015; Cleeves et al., 2015; Portilla-Revelo et al., 2022; Oberg 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, $\dot{M}=10^{-5}$ MJup yr-1. The planet point source emission
effectively changes the physical temperature in the surroundings in physical
scales up to $\sim$ 5 au at 48 au, i.e., at least one disk scale height. The
extension of this thermal effect is comparable to the size of the Hill’s
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 $\dot{M}>10^{-7}$ M⊙ yr-1, indistinct 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 does as well, producing an extra local heating effect.
For $\dot{M}=10^{-7}$ M⊙ 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 $\dot{M}>10^{-6}$M⊙ yr-1 produce strong thermal signatures,
i.e., the temperature increases by at least 20 K. Such thermal increases are
significant to change the chemical balance in the gaseous and solid
composition.

Figure 4: Dust temperature changes 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 $\dot{M}>10^{-7}M_{\mathrm{Jup}}$
yr-1 the thermal increase is noticeable while for accretion rates $\dot{M}\leq
10^{-7}\ M_{\mathrm{Jup}}$ yr-1 is probable that the effect is milder and not
significant even taking into account the blackbody emission of the protoplanet
at 1200 K.

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 surrounding material around the planet becomes
more transparent increasing the mean free path of emitted photons from the
planet or the CPD making the thermal gradient flatter, although the overall
extension 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 extension. 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 extension
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
as 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 of 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.

Figure 5: Dust temperature changes 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 $\dot{M}>10^{-7}$ MJup yr-1 the thermal increase is noticeable
while for accretion rates $\dot{M}\leq 10^{-7}\ M_{\mathrm{Jup}}$ yr-1 is
probable that the effect is milder and not significant even taking into
account the blackbody emission of the protoplanet at 1200 K.

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 extension of these changes and their absolute values.

Figure 6: Meridional cut at the azimuth where the planet is located. The
trends shift from the embedded planet compared to when is located at 48 au is
caused mainly by changes in the dust optical depth. Figure 7: Midplane
temperature changes with the simulation of a single planet embedded at 150 au.
Dust temperature changes when the planet is embedded beyond the edge of the
pebble disk at 150 au.

One of the main differences is the trend with the local dust-to-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
degrees over the surrounding material disappears. The probable reason behind
this phenomenon is the lack of scattering processes of 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.

###  3.2 Millimeter synthetic emission

In our synthetic predictions of the 1.3 millimeter dust continuum emission in
Figure 8, under the right combination of parameters, the CPD will not show a
significant difference in 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, there is 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 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, $\dot{M}>$10-5 MJup
yr-1, the emission enhancement are not above the 3% level. On the contrary,
when we assume a dust depletion in the CPD there is a significant decrement in
the local emission to the levels of 10%.

The same exercise for the single-planet at 150 au under ideal conditions,
i.e., without observational noises, would produce changes of $\sim$10%;
however, the disk is optically thin in millimeter wavelengths due to being
embedded beyond the edge of the pebble disk, having very low millimeter
optical depth. Hence, even though the planet and accretion irradiation can
produce significant relative changes, their absolute emission falls below the
sensitivity threshold of the 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 reported rms values
in the literature for deep high-resolution dust continuum observations with
ALMA in Figure 10 (Pérez et al., 2020; Benisty et al., 2021; Andrews et al.,
2021). The Figure shows that the expected deviations are marginally located at
the 3$\sigma$ level 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 70c, 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 reported values by our predictions.
From millimeter emission alone, it is hard to fully disentangle the accretion
rate in the CPD of PDS 70c, but it is consistent with a dust-undepleted CPD.
Nevertheless, multi-wavelength 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 are within detectable sensitivity thresholds
(Andrews et al., 2018). Even though dust substructures can have multiple
causes, under the assumption that a protoplanet is the cause of dust
substructures, the low detection rate of CPDs with current facilities and
surveys implies two possible scenarios:

  1. 1.

The first explanation is that either high accretion rate states are temporary
and relatively short, which can be explained by planets being formed in 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
103-105 yr. As the detection of more protoplanets and characterizations of
their masses and accretion rate become more ubiquitous, constraining the
timescales of their formation will be more accurate, which will lead to a
better understanding of planet formation.

  2. 2.

The second scenario is that planets accrete large grains in their CPDs on
relatively short scales. 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 70c (Benisty et al., 2021), it required deep integrations to
observe a massive planet in a severely depleted disk with confirmed accretion
tracers emission, 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 8: Comparison of synthetic 1.3-millimeter dust continuum images using
SIMIO between the different models and a standard model without protoplanetary
emission for $\sim$1 hour integration time. Except one case, there is no
strong additional signature of the CPD. Furthermore, when the CPD is severely
depleted of large grains, the dust millimeter emission is lower than the
surroundings. From the low accretion rate models we also infer that the 1200 K
blackbody protoplanetary spectrum does not produce a significant change in
synthetic ALMA images

If we take a closer look at the protoplanet’s orbital radius in the HD 163296
model by making an azimuthal cut at the separation of the 2 MJup, 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 depending on the accretion rate and
local dust-to-gas ratio. As mentioned above, when there is a strong dust
depletion within the Hill sphere of the planet, the local emission is lower,
matching the one 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 analog 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 the CPD.

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 dust continuum ALMA
observations. Nevertheless, such changes will manifest through line ratios and
thermal changes in gas emission, which have the potential to be detectable
under the right conditions.

Figure 9: Top: 1.3 millimeter 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 accretion rate that hides the presence of the CPD at millimeter
wavelengths, making it not noticeable Figure 10: Expected brightness
temperature deviations in our suite of models with respect to the baseline
emission compared with deep high-resolution 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$\sigma$
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 70c.

##  4 Discussion

###  4.1 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 they were present, should be easier to detect in 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 in the
short wavelengths with ALMA. Such an approach was further expanded by Binkert
et al. (2021) assuming a dust component through using a range of Stokes
numbers. However, the physical conditions of material surrounding protoplanets
can span over 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 of 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 (Perez 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 to
put at the very least an upper limit. The line emission asymmetries 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 millimeter dust continuum
emission probes 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 in 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, circumplanetary disks. Shocked gas during the accretion of
massive planets can reach very hot temperatures ($>$3000 K) producing
optically thick atomic blankets hiding the protoplanet emission (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 of 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 of 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 gas temperature increase which should have 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 of 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 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.

###  4.2 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 low accretion rates. A possible reason
behind the lack of CPD detections is that massive planets accrete the solid
reservoir in CPDs in 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 actually 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.

###  4.3 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 flat 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 thermodynamical 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 the Stokes number 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
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. (2022) and
Krapp et al. (2024b) 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 presented a new solver for particle-size
distribution studies in multifluid modules to study dust dynamics.

##  5 Summary

We present the thermal changes produced by including the emission of 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 post-processing 3D hydrodynamics with radiative transfer, we show that
the radiation of 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 of a hotter CPD. Hence, CPDs
brightness temperature may look hotter/colder than 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.

F.A, E.A.B. acknowledge support from NSF AAG Grant #1907653.

##  Appendix A Appendix information

## References

  * Alarcón et al. (2022) Alarcón, F., Bergin, E. A., & Teague, R. 2022, ApJ, 941, L24, doi: 10.3847/2041-8213/aca6e6
  * Alarcón et al. (2020) Alarcón, F., Teague, R., Zhang, K., Bergin, E. A., & Barraza-Alfaro, M. 2020, ApJ, 905, 68, doi: 10.3847/1538-4357/abc1d6
  * Andrews (2020) Andrews, S. M. 2020, ARA&A, 58, 483, doi: 10.1146/annurev-astro-031220-010302
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Andrews et al. (2021) Andrews, S. M., Elder, W., Zhang, S., et al. 2021, ApJ, 916, 51, doi: 10.3847/1538-4357/ac00b9
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Ayliffe & Bate (2009) Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 397, 657, doi: 10.1111/j.1365-2966.2009.15002.x
  * Bae et al. (2023) Bae, J., Isella, A., Zhu, Z., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Astronomical Society of the Pacific Conference Series, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 423
  * Bae et al. (2022) Bae, J., Teague, R., Andrews, S. M., et al. 2022, ApJ, 934, L20, doi: 10.3847/2041-8213/ac7fa3
  * Benisty et al. (2021) Benisty, M., Bae, J., Facchini, S., et al. 2021, ApJ, 916, L2, doi: 10.3847/2041-8213/ac0f83
  * Benítez-Llambay & Masset (2016) Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11, doi: 10.3847/0067-0049/223/1/11
  * Binkert et al. (2021) Binkert, F., Szulágyi, J., & Birnstiel, T. 2021, MNRAS, 506, 5969, doi: 10.1093/mnras/stab2075
  * Broome et al. (2023) Broome, M., Kama, M., Booth, R., & Shorttle, O. 2023, MNRAS, 522, 3378, doi: 10.1093/mnras/stad1159
  * Calahan et al. (2021) Calahan, J. K., Bergin, E. A., Zhang, K., et al. 2021, ApJS, 257, 17, doi: 10.3847/1538-4365/ac143f
  * Calvet & Gullbring (1998) Calvet, N., & Gullbring, E. 1998, ApJ, 509, 802, doi: 10.1086/306527
  * Calvet et al. (1994) Calvet, N., Hartmann, L., Kenyon, S. J., & Whitney, B. A. 1994, ApJ, 434, 330, doi: 10.1086/174731
  * Cleeves et al. (2015) Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015, ApJ, 807, 2, doi: 10.1088/0004-637X/807/1/2
  * Currie et al. (2022) Currie, T., Lawson, K., Schneider, G., et al. 2022, Nature Astronomy, 6, 751, doi: 10.1038/s41550-022-01634-x
  * Donehew & Brittain (2011) Donehew, B., & Brittain, S. 2011, AJ, 141, 46, doi: 10.1088/0004-6256/141/2/46
  * Dong et al. (2019) Dong, R., Liu, S.-Y., & Fung, J. 2019, ApJ, 870, 72, doi: 10.3847/1538-4357/aaf38e
  * Draine (2003) Draine, B. T. 2003, ApJ, 598, 1017, doi: 10.1086/379118
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, record ascl:1202.015.  http://ascl.net/1202.015
  * Facchini et al. (2018) Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104, doi: 10.1051/0004-6361/201731390
  * Fung et al. (2019) Fung, J., Zhu, Z., & Chiang, E. 2019, ApJ, 887, 152, doi: 10.3847/1538-4357/ab53da
  * Garrido-Deutelmoser et al. (2023) Garrido-Deutelmoser, J., Petrovich, C., Charalambous, C., Guzmán, V. V., & Zhang, K. 2023, ApJ, 945, L37, doi: 10.3847/2041-8213/acbea8
  * Gratton et al. (2019) Gratton, R., Ligi, R., Sissa, E., et al. 2019, A&A, 623, A140, doi: 10.1051/0004-6361/201834760
  * Haffert et al. (2019) Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nature Astronomy, 3, 749, doi: 10.1038/s41550-019-0780-5
  * Hammond et al. (2023) Hammond, I., Christiaens, V., Price, D. J., et al. 2023, MNRAS, 522, L51, doi: 10.1093/mnrasl/slad027
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Henning & Stognienko (1996) Henning, T., & Stognienko, R. 1996, A&A, 311, 291
  * Henyey & Greenstein (1941) Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70, doi: 10.1086/144246
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Isella et al. (2019) Isella, A., Benisty, M., Teague, R., et al. 2019, ApJ, 879, L25, doi: 10.3847/2041-8213/ab2a12
  * Isella et al. (2018) Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49, doi: 10.3847/2041-8213/aaf747
  * Izquierdo et al. (2022) Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2, doi: 10.3847/1538-4357/ac474d
  * Izquierdo et al. (2023) Izquierdo, A. F., Testi, L., Facchini, S., et al. 2023, arXiv e-prints, arXiv:2304.03607, doi: 10.48550/arXiv.2304.03607
  * Keppler et al. (2018) Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44, doi: 10.1051/0004-6361/201832957
  * Krapp et al. (2024a) Krapp, L., Garrido-Deutelmoser, J., Benítez-Llambay, P., & Kratter, K. M. 2024a, ApJS, 271, 7, doi: 10.3847/1538-4365/ad14f9
  * Krapp et al. (2022) Krapp, L., Kratter, K. M., & Youdin, A. N. 2022, ApJ, 928, 156, doi: 10.3847/1538-4357/ac5899
  * Krapp et al. (2024b) Krapp, L., Kratter, K. M., Youdin, A. N., et al. 2024b, arXiv e-prints, arXiv:2402.14638, doi: 10.48550/arXiv.2402.14638
  * Law et al. (2023) Law, C. J., Booth, A. S., & Öberg, K. I. 2023, ApJ, 952, L19, doi: 10.3847/2041-8213/acdfd0
  * Lubow & Martin (2012) Lubow, S. H., & Martin, R. G. 2012, ApJ, 749, L37, doi: 10.1088/2041-8205/749/2/L37
  * Masset (2000) Masset, F. 2000, A&AS, 141, 165, doi: 10.1051/aas:2000116
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
  * Montesinos et al. (2015) Montesinos, M., Cuadra, J., Perez, S., Baruteau, C., & Casassus, S. 2015, ApJ, 806, 253, doi: 10.1088/0004-637X/806/2/253
  * Öberg et al. (2021) Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1, doi: 10.3847/1538-4365/ac1432
  * Oberg et al. (2022) Oberg, N., Kamp, I., Cazaux, S., Woitke, P., & Thi, W. F. 2022, A&A, 667, A95, doi: 10.1051/0004-6361/202244092
  * Paardekooper et al. (2023) Paardekooper, S., Dong, R., Duffell, P., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Astronomical Society of the Pacific Conference Series, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 685
  * Pérez et al. (2018) Pérez, S., Casassus, S., & Benítez-Llambay, P. 2018, MNRAS, 480, L12, doi: 10.1093/mnrasl/sly109
  * Perez et al. (2015) Perez, S., Dunhill, A., Casassus, S., et al. 2015, ApJ, 811, L5, doi: 10.1088/2041-8205/811/1/L5
  * Pérez et al. (2020) Pérez, S., Casassus, S., Hales, A., et al. 2020, ApJ, 889, L24, doi: 10.3847/2041-8213/ab6b2b
  * Pinte et al. (2019) Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nature Astronomy, 3, 1109, doi: 10.1038/s41550-019-0852-6
  * Pinte et al. (2020) Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJ, 890, L9, doi: 10.3847/2041-8213/ab6dda
  * Portilla-Revelo et al. (2022) Portilla-Revelo, B., Kamp, I., Rab, C., et al. 2022, A&A, 658, A89, doi: 10.1051/0004-6361/202141764
  * Rosotti et al. (2016) Rosotti, G. P., Juhasz, A., Booth, R. A., & Clarke, C. J. 2016, MNRAS, 459, 2790, doi: 10.1093/mnras/stw691
  * Sandstrom et al. (2013) Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5, doi: 10.1088/0004-637X/777/1/5
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Szulágyi et al. (2022) Szulágyi, J., Binkert, F., & Surville, C. 2022, ApJ, 924, 1, doi: 10.3847/1538-4357/ac32d1
  * Szulágyi et al. (2016) Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853, doi: 10.1093/mnras/stw1160
  * Szulágyi et al. (2014) Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65, doi: 10.1088/0004-637X/782/2/65
  * Szulágyi & Mordasini (2017) Szulágyi, J., & Mordasini, C. 2017, MNRAS, 465, L64, doi: 10.1093/mnrasl/slw212
  * Szulágyi et al. (2018) Szulágyi, J., Plas, G. v. d., Meyer, M. R., et al. 2018, MNRAS, 473, 3573, doi: 10.1093/mnras/stx2602
  * Teague et al. (2019) Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378, doi: 10.1038/s41586-019-1642-0
  * Teague et al. (2021) Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18, doi: 10.3847/1538-4365/ac1438
  * van der Marel et al. (2018) van der Marel, N., Williams, J. P., & Bruderer, S. 2018, ApJ, 867, L14, doi: 10.3847/2041-8213/aae88e
  * van der Velden (2020) van der Velden, E. 2020, The Journal of Open Source Software, 5, 2004, doi: 10.21105/joss.02004
  * Warren & Brandt (2008) Warren, S. G., & Brandt, R. E. 2008, Journal of Geophysical Research: Atmospheres, 113, doi: https://doi.org/10.1029/2007JD009744
  * Wolf & D’Angelo (2005) Wolf, S., & D’Angelo, G. 2005, ApJ, 619, 1114, doi: 10.1086/426662
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47, doi: 10.3847/2041-8213/aaf744
  * Zhou et al. (2022) Zhou, Y., Sanghi, A., Bowler, B. P., et al. 2022, ApJ, 934, L13, doi: 10.3847/2041-8213/ac7fef
  * Zhu (2015) Zhu, Z. 2015, ApJ, 799, 16, doi: 10.1088/0004-637X/799/1/16
  * Zhu et al. (2018) Zhu, Z., Andrews, S. M., & Isella, A. 2018, MNRAS, 479, 1850, doi: 10.1093/mnras/sty1503
  * Zhu et al. (2012) Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ, 755, 6, doi: 10.1088/0004-637X/755/1/6


# Histogram- and Diffusion-Based Medical Out-of-Distribution Detection

Evi M.C. Huijben∗ Department of Biomedical Engineering, Eindhoven University
of Technology, Eindhoven, The Netherlands  Sina Amirrajab E.M.C. Huijben and
S. Amirrajab contributed equally Department of Biomedical Engineering,
Eindhoven University of Technology, Eindhoven, The Netherlands  Josien P.W.
Pluim  Department of Biomedical Engineering, Eindhoven University of
Technology, Eindhoven, The Netherlands

###### Abstract

Out-of-distribution (OOD) detection is crucial for the safety and reliability
of artificial intelligence algorithms, especially in the medical domain. In
the context of the Medical OOD (MOOD) detection challenge 2023, we propose a
pipeline that combines a histogram-based method and a diffusion-based method.
The histogram-based method is designed to accurately detect homogeneous
anomalies in the toy examples of the challenge, such as blobs with constant
intensity values. The diffusion-based method is based on one of the latest
methods for unsupervised anomaly detection, called DDPM-OOD. We explore this
method and propose extensive post-processing steps for pixel-level and sample-
level anomaly detection on brain MRI and abdominal CT data provided by the
challenge. Our results show that the proposed DDPM method is sensitive to blur
and bias field samples, but faces challenges with anatomical deformation,
black slice, and swapped patches. These findings suggest that further research
is needed to improve the performance of DDPM for OOD detection in medical
images.

##  1 Introduction

Out-of-distribution (OOD) detection, the task of identifying inputs that do
not belong to the data distribution on which the model was trained [15], is
crucial for ensuring the safety and reliability of machine learning systems in
the medical domain, where unexpected inputs can lead to misdiagnosis and
adverse outcomes. Generative models are powerful tools for learning the
underlying distribution of complex medical data, which in turn can be used for
OOD detection [16, 18, 2, 7, 13].

Recently, diffusion models, introduced as a powerful category of generative
models [4, 6], have shown impressive results in generating realistic and
diverse images from various domains in computer vision [3] and have found
their way into unsupervised anomaly detection in medical images [14, 5, 12, 8,
1].

To develop and benchmark methods for OOD detection, the Medical Out-of-
Distribution Analysis Challenge (MOOD)111http://medicalood.dkfz.de/web/ [17],
which includes sample-level and pixel-level anomaly detection for brain MRI
and abdominal CT datasets, was first organized in conjunction with MICCAI 2020
and will be organized again in conjunction with MICCAI 2023.

##  2 Dataset

The MOOD2023 dataset consists of 800 brain MRI scans ($256\times 256\times
256$) and 550 abdominal CT scans ($512\times 512\times 512$), which we
randomly divided into a training (90%) and validation (10%) set. Additionally,
a small toy set with one in-distribution sample and three samples with a
homogeneous sphere inserted to represent an anomaly was provided for both
regions.

To evaluate our model, we transformed 40 in-distribution validation set
samples for both regions by applying the transformations described in Table 1
four times.

Table 1: Transformations used to create the OOD validation sets (based on torchio.transforms, except for RandomBlackSlice). Transformation | Parameters | Value brain | Value abdomen
---|---|---|---
RandomElasticDeformation | max_displacement | 30 | 50
max_displacement | 40 | 80
RandomBlur | std | 2 | 3
std | 4 | 5
RandomBiasField | coefficients | 1 | 1
coefficients | 2 | 2
RandomSwap (1 iteration) | patch_size | 30 | 50
patch_size | 80 | 100
RandomBlackSlice | slice_thickness | 1 | 1
slice_thickness | 5 | 7

##  3 Method

Figure 1: Anomaly detection pipeline using histogram OOD detection method
(blue) and diffusion OOD detection method (green). The same pipeline is used
for brain MRI and abdominal CT data.

The proposed pipeline is shown in Figure 1222Code available at
https://github.com/evihuijben/mood_challenge. In case of the abdominal
dataset, the input images are first resized to $256\times 256\times 256$
before being processed through the pipeline. The histogram OOD detection
module is then used for each input image to detect homogeneous anomalies,
shown in blue. These anomalies are identified as specific intensity peaks in
the image histogram. If no anomalies are detected by the histogram module, the
image is passed on to the diffusion module. At this stage, the reconstruction
error is measured using a pixel-wise structural similarity index measure
(SSIM) map. The output is then post-processed to generate the final
prediction, which is shown in green.

###  3.1 Histogram OOD detection

We propose a histogram-based OOD detection method that can detect only
homogeneous anomalies, such as the toy examples in the challenge. We first
compute the histogram mean and standard deviation (4096 bins) for the brain
and abdomen based on the 3D training images. When calculating the histogram of
an OOD test sample, homogeneous anomalies result in a peak with a much higher
frequency than in the training set. As shown in Figure 2, the mean of the
training set with four times the standard deviation is used to suppress the
histogram of normal brain areas, leaving only the peak at intensity value 0.24
for toy example one. One spike appeared consistently for all brain cases, as
seen in green in Figure 2. We found that this spike originated from
homogeneous masking artifacts at the border of the applied brain mask.
However, since this bin has a large standard deviation, it is never detected.
In addition, we discard the zero intensity peak to avoid detection due to the
large homogeneous background. Finally, we increase the standard deviations to
64 and 128 for brain and pelvis, respectively, to avoid false positives in
intensities with little variation. The final detected peak is used to
threshold the image and create an anomaly mask. We further process the mask
using erosion and dilation operations with a structure size of six to remove
small erroneous objects.

Figure 2: Histogram OOD detection to detect homogeneous regions in toy
examples. All histograms are computed 3D images. We first subtract the
histogram of the toy example from the mean histogram of the training data, and
then detect the peak location (intensity value of 0.24) above four times the
standard deviation.

###  3.2 Diffusion OOD detection

####  3.2.1 Diffusion model

Diffusion models are latent variable models, a variant of variational
autoencoders, where the latent variable matches the dimension of the input
data [6]. They include a forward diffusion process that incrementally adds
noise to the data, and a backward diffusion process that uses deep learning-
based denoising to generate realistic images.

Our diffusion-based OOD detection method is based on the work of Graham et al.
[5]333https://github.com/marksgraham/ddpm-ood, employing denoising diffusion
probabilistic models (DDPMs) [6] for reconstruction-driven OOD detection. The
key concept involves training the DDPM on the training data and then use it to
reconstruct input images at test time. If the input image is within the
training distribution, the model can reconstruct a high quality image.
However, if the input image is OOD, the model yields a poor reconstruction,
which is detectable by a reconstruction score. We trained a time-conditional
2D U-Net with attention mechanisms based on [11, 10] to denoise images at each
step in the backward process. Two identical models were trained independently
for brain and abdomen, using the parameters described in Table 2.

Table 2: Diffusion model parameters, built on top of [5] and MONAI generative AI [9]. Noise | Simplex [14]
---|---
Diffusion steps | 1000
Scheduler | Scaled linear [11]
($\beta_{start}$,$\beta_{end}$) | (0.001 , 0.015)
Optimizer | Adam
Learning rate | $2.5e-5$
Loss | $MSE({noise}_{added},{noise}_{predicted})$
Batch size | 4
Epoch | 60
Intensity normalization | [0, 1]
Data augmentation | No
Input size (brain & abdomen) | 2D axial slices $265\times 256$

####  3.2.2 Inference and post-processing

Given the trained DDPM, we can perform multiple reconstructions at different
diffusion steps, as suggested in [5], or determine the optimal step, as
proposed in [14]. Smaller steps (e.g., t=10 or 50) introduce minimal noise,
resulting in anatomically similar reconstructions with faithful representation
of anomalous regions. Conversely, larger steps (e.g., t=500 or 800) heavily
distort the image content, causing anomalous regions to heal during
reconstruction. Moreover, the reconstruction time becomes unacceptably long
with larger steps. After experimental evaluation on our OOD validation set
(see section 2), we chose t=200.

To obtain an OOD prediction score from the reconstruction, we first compute a
body mask using SimpleITK algorithms. This involves Otsu thresholding on the
input image, followed by binary dilation (radius = 5) and selection of the
largest connected component. We then apply binary morphological closure
(radius = 5) and hole filling to ensure complete object coverage. This body
mask is then used to process the reconstruction image into a pixel-level OOD
prediction. To achieve this, we compute a per-pixel SSIM map that compares the
3D masked reconstruction to the 3D input image. A 3-pixel border padding is
added to the SSIM map to ensure accurate SSIM computation while minimizing
border effects. Subsequently, we apply a Gaussian filter with a standard
deviation of 15 to blur the SSIM map and threshold it at 0.5. The resulting
binarized 3D map serves as the pixel-level prediction score for the MOOD
challenge, with the sample-level OOD score set to 1 if at least one voxel is
predicted to be OOD. The pixel-level predictions for the abdominal dataset are
resized to the original resolution in the final post-processing step.

Figure 3: The qualitative results for pixel-level and sample-level predictions
using our proposed histogram OOD detection (blue) and diffusion OOD detection
(green). Figure 4: Sample-level sensitivity grouped by OOD type for brain
cases. Figure 5: Sample-level sensitivity grouped by OOD type for abdomen
cases.

##  4 Results

Figure 3 shows qualitative results from a selection of the toy and OOD
validation datasets for brain images. The first row (in blue) illustrates a
toy example correctly identified by our histogram-based OOD detection. The
last three rows (in green) depict one in-distribution and two OOD samples. As
observed, the in-distribution sample is correctly classified, while the
blurred example is accurately identified as OOD. However, our method fails to
detect the elastic deformation. Figure 4 presents the sample-level sensitivity
categorized by different OOD groups, demonstrating that toy, blur, and bias
field samples are easily detected, while deformation, swap, and black slice
samples pose greater challenges. Furthermore, our method achieves a sample-
level specificity of 1.00 for the in-distribution brain validation set.
Similar patterns are observed for the abdomen dataset, both visually and
quantitatively, as depicted in Figure 5, which shows the sensitivity of the
OOD groups. Additionally, the sample-level specificity yields a similar value
of 0.96.

##  5 Discussion and Conclusion

In this paper, we introduced a two-step approach for OOD detection in medical
images: 1) histogram-based and 2) diffusion model-based. The histogram-based
method excelled in identifying homogeneous anomalies with high accuracy, while
the diffusion model-based approach detected OOD cases missed by the histogram-
based method. However, we observed that the diffusion model lacked anatomical
awareness, resulting in poor detection of cases involving deformations (e.g.
elastic deformations or random patch swaps).

For future research, we recommend exploring anatomically aware models by using
2.5D or 3D models, and utilizing models that better capture spatial features,
such as latent diffusion models [11]. Additionally, since real-time OOD
detection is essential, addressing the longer inference time inherent in
diffusion models should be a priority for future work.

In conclusion, our proposed two-step OOD approach effectively identifies
specific OOD cases in brain MRI and abdominal CT scans, which has the
potential to enhance the application of machine learning in healthcare.

## References

  * [1] Finn Behrendt, Debayan Bhattacharya, Julia Krüger, Roland Opfer, and Alexander Schlaefer.  Patched diffusion models for unsupervised anomaly detection in brain mri.  arXiv preprint arXiv:2303.03758, 2023.
  * [2] Xiaoran Chen and Ender Konukoglu.  Unsupervised detection of lesions in brain mri using constrained adversarial auto-encoders.  arXiv preprint arXiv:1806.04972, 2018.
  * [3] Florinel-Alin Croitoru, Vlad Hondru, Radu Tudor Ionescu, and Mubarak Shah.  Diffusion models in vision: A survey.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2023\.
  * [4] Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  Advances in neural information processing systems, 34:8780–8794, 2021.
  * [5] Mark S Graham, Walter HL Pinaya, Petru-Daniel Tudosiu, Parashkev Nachev, Sebastien Ourselin, and Jorge Cardoso.  Denoising diffusion models for out-of-distribution detection.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2947–2956, 2023.
  * [6] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  Advances in neural information processing systems, 33:6840–6851, 2020.
  * [7] Sergio Naval Marimont and Giacomo Tarroni.  Anomaly detection through latent space restoration using vector quantized variational autoencoders.  In 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI), pages 1764–1767. IEEE, 2021.
  * [8] Walter HL Pinaya, Mark S Graham, Robert Gray, Pedro F Da Costa, Petru-Daniel Tudosiu, Paul Wright, Yee H Mah, Andrew D MacKinnon, James T Teo, Rolf Jager, et al.  Fast unsupervised brain anomaly detection and segmentation with diffusion models.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 705–714. Springer, 2022.
  * [9] Walter HL Pinaya, Mark S Graham, Eric Kerfoot, Petru-Daniel Tudosiu, Jessica Dafflon, Virginia Fernandez, Pedro Sanchez, Julia Wolleb, Pedro F da Costa, Ashay Patel, et al.  Generative ai for medical imaging: extending the monai framework.  arXiv preprint arXiv:2307.15208, 2023.
  * [10] Walter HL Pinaya, Petru-Daniel Tudosiu, Jessica Dafflon, Pedro F Da Costa, Virginia Fernandez, Parashkev Nachev, Sebastien Ourselin, and M Jorge Cardoso.  Brain imaging generation with latent diffusion models.  In MICCAI Workshop on Deep Generative Models, pages 117–126. Springer, 2022.
  * [11] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 10684–10695, 2022.
  * [12] Julia Wolleb, Florentin Bieder, Robin Sandkühler, and Philippe C Cattin.  Diffusion models for medical anomaly detection.  In International Conference on Medical image computing and computer-assisted intervention, pages 35–45. Springer, 2022.
  * [13] Julia Wolleb, Robin Sandkühler, and Philippe C Cattin.  Descargan: Disease-specific anomaly detection with weak supervision.  In Medical Image Computing and Computer Assisted Intervention–MICCAI 2020: 23rd International Conference, Lima, Peru, October 4–8, 2020, Proceedings, Part IV 23, pages 14–24. Springer, 2020.
  * [14] Julian Wyatt, Adam Leach, Sebastian M Schmon, and Chris G Willcocks.  Anoddpm: Anomaly detection with denoising diffusion probabilistic models using simplex noise.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 650–656, 2022.
  * [15] Jingkang Yang, Kaiyang Zhou, Yixuan Li, and Ziwei Liu.  Generalized out-of-distribution detection: A survey.  arXiv preprint arXiv:2110.11334, 2021.
  * [16] Chong Zhou and Randy C Paffenroth.  Anomaly detection with robust deep autoencoders.  In Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining, pages 665–674, 2017.
  * [17] David Zimmerer, Peter M Full, Fabian Isensee, Paul Jäger, Tim Adler, Jens Petersen, Gregor Köhler, Tobias Ross, Annika Reinke, Antanas Kascenas, et al.  Mood 2020: A public benchmark for out-of-distribution detection and localization on medical images.  IEEE Transactions on Medical Imaging, 41(10):2728–2738, 2022.
  * [18] David Zimmerer, Simon AA Kohl, Jens Petersen, Fabian Isensee, and Klaus H Maier-Hein.  Context-encoding variational autoencoder for unsupervised anomaly detection.  arXiv preprint arXiv:1812.05941, 2018.