# Radiative Scale Height and Shadows in Protoplanetary Disks

Matías Montesinos, Nicolás Cuello, Johan Olofsson, Jorge Cuadra, Amelia Bayo, Gesa H. -M. Bertrang, Clément Perrot

► To cite this version:

Matías Montesinos, Nicolás Cuello, Johan Olofsson, Jorge Cuadra, Amelia Bayo, et al. Radiative Scale Height and Shadows in Protoplanetary Disks. The Astrophysical Journal, 2021, 910, 10.3847/1538-4357/abe3fc. insu-03705322

HAL Id: insu-03705322

https://insu.hal.science/insu-03705322v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Radiative Scale Height and Shadows in Protoplanetary Disks

Matías Montesinos¹,²,³ ID, Nicolás Cuello³,⁴,⁵, Johan Olofsson¹,³ ID, Jorge Cuadra³,⁶, Amelia Bayo¹,³ ID,
Gesa H.-M. Bertrang⁷ ID, and Clément Perrot¹,³,⁸

¹Instituto de Física y Astronomía, Universidad de Valparaíso, Chile; matías.montesinos@uv.cl

²Chinese Academy of Sciences South America Center for Astronomy, National Astronomical Observatories, CAS, Beijing 100012, People's Republic of China

³Núcleo Milenio de Formación Planetaria (NPF), Chile

⁴Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁵Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Santiago, Chile

⁶Departamento de Ciencias, Facultad de Artes Liberales, Universidad Adolfo Ibáñez, Av. Padre Hurtado 750, Viña del Mar, Chile

⁷Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

⁸LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France

Received 2020 September 15; revised 2021 January 29; accepted 2021 February 4; published 2021 March 24

## Abstract

Planets form in young circumstellar disks called protoplanetary disks. However, it is still difficult to catch planet formation in situ. Nevertheless, from recent ALMA/SPHERE data, encouraging evidence of the direct and indirect presence of embedded planets has been identified in disks around young stars: co-moving point sources, gravitational perturbations, rings, cavities, and emission dips or shadows cast on disks. The interpretation of these observations needs a robust physical framework to deduce the complex disk geometry. In particular, protoplanetary disk models usually assume the gas pressure scale height given by the ratio of the sound speed over the azimuthal velocity $H/r = c_s/v_k$. By doing so, radiative pressure fields are often ignored, which could lead to a misinterpretation of the real vertical structure of such disks. We follow the evolution of a gaseous disk with an embedded Jupiter-mass planet through hydrodynamical simulations, computing the disk scale height including radiative pressure, which was derived from a generalization of the stellar atmosphere theory. We focus on the vertical impact of the radiative pressure in the vicinity of circumplanetary disks, where temperatures can reach $\ge 1000$ K for an accreting planet and radiative forces can overcome gravitational forces from the planet. The radiation pressure effects create a vertical, optically thick column of gas and dust at the protoplanet location, casting a shadow in scattered light. This mechanism could explain the peculiar illumination patterns observed in some disks around young stars such as HD 169142 where a moving shadow has been detected or the extremely high aspect ratio $H/r \sim 0.2$ observed in systems like AB Aur and CT Cha.

Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Hydrodynamical simulations (767); Planetary-disk interactions (2204); Radiative transfer simulations (1967)

## 1. Introduction

During the past few years, high angular resolution imaging of disks around young stars has revealed different patterns in the gas and dust structures. These structures are usually observed in gas-rich disks with an inner cavity around low- and intermediate-mass stars, classified as transition disks. Spatially resolved observations reveal astonishing images of circumstellar disks. For instance, in the disk around HL Tau, ALMA unveiled for the first time a series of gaps at millimeter wavelengths (ALMA Partnership et al. 2015). Several other remarkable examples from high angular resolution images (at near-infrared and sub-millimeter wavelengths) show multiple gaps, spirals, and rings. It is thought that embedded planets cause most of these features (e.g., Long et al. 2018a; Keppler et al. 2018). A complete review of disk substructures' observations at high angular resolution can be found in the DSHARP project paper series (Andrews et al. 2018; Huang et al. 2018a, 2018b).

Besides these astonishing gas/dust structures, intriguing illumination features—both in scattered light and thermal (sub) millimeter emission—have been observed. For instance, Avenhaus et al. (2014) reported two intensity nulls seen in the circumbinary disk of HD 142527. These were later identified as shadows cast by a tilted inner circumprimary disk (Casassus et al. 2015a; Marino et al. 2015). Moreover, multiple shadow features had been discovered through high-contrast polarimetric differential imaging

with VLT/SPHERE in the outer disk of HD 135344B (Stolker et al. 2016) and the disk around TW Hya (Debes et al. 2017).

Interestingly, Quanz et al. (2013) reported VLT/NACO observations of the disk around the isolated Herbig Ae/Be star HD 169142 showing a narrow emission dip located at $\sim 80^{\circ}$. More recently, Bertrang et al. (2018, 2020) presented VLT/SPHERE/ZIMPOL observations of the same system at higher spatial resolution. The latter reveals an inner ring at 24 au from the star and a narrow emission dip in surface brightness located now at $\sim 50^{\circ}$. Remarkably, the moving shadow could be caused by a rotating optically thick bump possibly located at $\sim 12$ au from the star, as suggested by Bertrang et al. (2018).

Most of these features seem to be related to planet formation processes, in which planet-disk interactions sculpt the disk surface and midplane. However, direct detection of planets embedded in disks remains elusive. To date, only one planetary-mass companion has been imaged coexisting with a protoplanetary disk (Keppler et al. 2018; Müller et al. 2018), corresponding to the discovery of a forming planet within the gap in PDS 70, while some planet candidates have been proposed through kinematic signatures (Pinte et al. 2019, 2020).

An important quantity that characterizes the circumstellar disk geometry is its pressure scale height, which is needed to explain several observations (especially the non-resolved ones). For instance, transitional disks exhibit an excess of NIR emission related to unusually high aspect ratios ($H/r \sim 0.2$) at
the inner rim, where a puffed-up wall is created by dust evaporation (Dullemond et al. 2001; Natta et al. 2001; Dullemond & Monnier 2010; Olofsson et al. 2013). There is also some inferred difference in the aspect ratio between the disks of groups I and II stars (Meeus et al. 2001). Variations on the vertical scale height in flat or only moderately flared disks are also supposed to be responsible for self-shadowing effects (Garufi et al. 2014; Stolker et al. 2016)

The vertical disk extension is directly related to the fluid pressure field and gravity. In general, circumstellar gaseous disks are expected to be gas pressure-dominated. This condition is used to assume equilibrium between gravity and gas pressure, leading to the well-known relation for the aspect ratio $H/R = c_s/v_φ$, where $c_s$ and $v_φ$ are the sound speed and the azimuthal velocity of the gas, respectively.

However, it is worth asking if this condition holds everywhere in the disk—especially at high temperatures. In particular, upon which circumstances (if any) does the scale height formula break down? In such a case, what kind of equation should be used instead?

One important modification should be implemented in the vicinity of a planet where its gravity dominates and must be taken into account (Müller et al. 2012). Also, very high temperatures are expected to develop in some regions of protoplanetary disks. One of these hot regions arises at the inner boundaries of transitional disks where the dust sublimates, reaching temperatures in the range of 1000–1500 K, which produces a puffed-up inner rim (Natta et al. 2001). Also, highly luminous ($10^{-6}-10^{-3} L_⊙$) forming planets are expected to form during early stages (~1 Myr) of planet formation (Mordasini et al. 2012). From numerical simulations of disks with embedded luminous planets, Montesinos et al. (2015) showed that the circumplanetary disk (CPD) region reaches temperatures of ~1000 K or more, with thermal emission peaking in L' in concordance with observations (Spiegel & Burrows 2012; Zhu et al. 2015).

At high temperatures (≥1000 K), the radiative pressure may become relevant, where processes similar to stellar atmospheres arise (Mihalas & Weibel Mihalas 1984). A simplified calculation of the scale height, including radiative pressure, can be found in Montesinos & de Freitas Pacheco (2011). In that work, the authors also consider contributions from pressure fields related to viscous turbulence and the disk self-gravity, in the context of accretion disks around super-massive black holes.

In this paper, we adopt a generalization of the gray model for stellar atmospheres applied to disk geometry introduced by Hubeny (1990) to include radiation pressure when computing the vertical structure of the disk. We also include the gravitational influence of an embedded planet in that calculation. Around hot accreting planets, radiative forces can dominate over the planet's gravity, creating an optically thick bump above the circumplanetary region. We link such bumps with self-shadowing effects through radiative transfer calculations. This could explain, for instance, the moving shadow observed in HD 169142 by associating it with the presence of an inner planet (Bertrang et al. 2018, 2020).

In Section 2, we discuss the physical background of the atmosphere model for computing the disk scale height. In Section 3, we run numerical hydro-simulations following the evolution of a protoplanetary disk with the adopted formulation. In Section 4, we use our hydro-models to interpret observations

by performing radiative transfer calculations on them. Finally, in Section 5 we discuss our results, which are summarized in Section 6.

## 2. Scale height

The vertical structure of a disk can be computed under the assumption of hydrostatic equilibrium between pressure forces and gravity, i.e., the vertical gradient pressure is balanced by the vertical tidal acceleration toward the midplane. Then, we have:

$$ \frac{dP}{dz} = -\int \rho \nabla \Psi dz, \qquad (1) $$

where $z$ corresponds to the distance from the midplane, $P$ and $\rho$ are the total pressure and density, respectively, and $\Psi$ is the gravitational potential given by contributions from the star and the embedded planet, namely, $\Psi = \Psi_* + \Psi_{\text{planet}}$.

### 2.1. Scale height Near the Planet

Taking into account the gravitational potential from the star and the planet, the vertical hydrostatic Equation (1) can be written in cylindrical coordinates as:

$$ \frac{1}{\rho} \frac{\partial P}{\partial z} = - \frac{GM_*}{(r^2 + z^2)^{(3/2)}} - \frac{GM_p}{(s^2 + z^2)^{(3/2)}}, \qquad (2) $$

where $M_*$, and $M_p$ are the mass of the star and planet, respectively, and $s$ is the distance from the projected position of the planet ($z=0$ plane) to a disk element.

Assuming pure gas pressure given by $P_g = \rho c_s^2$, the scale height $H$ in the absence of a planet is defined as:

$$ H = \frac{c_s}{\Omega}, \qquad (3) $$

where $c_s$ is the vertically isothermal sound speed and $\Omega$ is the angular velocity of the disk.

The embedded planet changes the disk structure, leading to a reduced thickness near it. We follow the procedure by Müller et al. (2012) to compute the scale height $H_{\text{planet}}$, and the density $\rho_{\text{planet}}$ of the disk in the vicinity of the planet. For a vertically isothermal disk, Equation (2) can be integrated:

$$ \rho_{\text{planet}} = \rho_0 \exp \left\{ - \left( \frac{1}{2} \frac{z^2}{H^2} + \frac{|z|}{H_{\text{planet}}} \right) \right\}, \qquad (4) $$

where the scale height near the planet is given by,

$$ H_{\text{planet}} = \frac{4s^2H^2}{qr^3}, \qquad (5) $$

where $q = M_p/M_*$, and $H$ defines the scale height in the absence of the gravitational influence of the planet (Equation (3)).

The planet influence is characterized by the transition distance (Müller et al. 2012):

$$ s_t = \frac{1}{2} \left( \frac{r^3 q}{H} \right)^{1/2}. \qquad (6) $$

When $s < s_t$, the planet dominates and therefore $H = H_{\text{planet}}$ (Equation (5)); otherwise $H$ is given by Equation (3).
## 2.2. Radiative Scale height

In some cases, radiation pressure may be relevant, such as in circumplanetary disks around accreting planets where temperatures could reach high temperatures above 1000 K. We use an analytical model of the scale height derived by Hubeny (1990) to compute the vertical structure, obtained from a generalization of the classical stellar atmospheric theory applied to a plane-parallel geometry compatible with a disk structure. We combine this procedure with the calculations of the scale height in the vicinity of the planet (Equation (5)).

The total pressure is assumed to be the contribution of the gas and radiation pressure, i.e., $P = P_g + P_r$, respectively. The gas pressure can be written as:

$$P_g = \rho c_s^2. \qquad (7)$$

It is important to note that the isothermal sound speed $c_s$ is related to the gas pressure $P_g$ only and not to the total pressure $P$. The radiation pressure $P_r$ can be written as:

$$P_r = \frac{4\pi}{c} K_{\nu}, \qquad (8)$$

where $c$ is the speed of light, and $K_{\nu}$ is the second momentum radiation flux (proportional to the radiation pressure tensor) and defined as:

$$K_{\nu} \equiv \frac{1}{4\pi} \int I_{\nu} \mu^2 d\Omega, \qquad (9)$$

where $I_{\nu}$ is the radiation intensity field in the solid angle $\Omega$ and its direction is specified by $\mu \equiv \cos(\theta)$ ($\theta$ being the elevation angle respect to the surface normal). For a semi-isotropic radiation field we have $\mu \approx 1.9$.

Assuming an energy balance in the vertical direction in which the energy dissipated per unit volume $Q^+$ equals the net radiation loss per unit volume $Q^-$, from the total pressure $P$ of the fluid and the definition of the moment equation $K_{\nu}$ from the stellar atmospheres theory (Mihalas 1978), the hydrostatic equilibrium condition in Equation (1) can be rewritten as a density equation (Hubeny 1990):

$$\frac{d\rho}{dz} = -\frac{2}{H_g^2}\rho z + \frac{2H_r}{H_g^2}\rho[1 - \theta(z)], \qquad (10)$$

where $H_g$ is defined as the pure gas pressure scale height, given by:

$$H_g \equiv \begin{cases} H_{\text{planet}} & (\text{Equation (5)}), \quad \text{if } s < s_t \\ H_g & (\text{Equation (3)}), \quad \text{if } s \ge s_t, \end{cases} \qquad (11)$$

and $H_r$ represents a pure radiative pressure scale height, given by:

$$H_r \equiv (\sigma/c) T_{\text{eff}}^4 \kappa / \Omega_K^2, \qquad (12)$$

where $T_{\text{eff}}$ is the disk effective temperature, $\Omega_K$ the Keplerian velocity, $\kappa$ the flux mean opacity, and $\sigma$ the Stefan-Boltzmann constant. The effective temperature is related with the midplane temperature through:

$$T_{\text{eff}}^4 = \frac{T^4}{\tau_{\text{eff}}}, \qquad (13)$$

where we use for the effective optical depth $\tau_{\text{eff}}$ (valid for the optically thick case, Hubeny 1990):

$$\tau_{\text{eff}} = \frac{\sqrt{3}}{4} + \frac{3\tau}{8} + \frac{1}{4\tau}, \qquad (14)$$

and the optical depth is obtained from the mean opacity $\kappa$:

$$\tau = \frac{1}{2} \kappa \rho. \qquad (15)$$

To solve Equation (10) one can use a monotonically increasing function between 0 and 1 for $\theta(z)$:

$$\theta(z) = \begin{cases} 1 - (z/H), & \text{if } z < H \\ 0, & \text{if } z \ge H, \end{cases} \qquad (16)$$

where $H$ may be called the density scale height, given by:

$$H = \Sigma / (\sqrt{2\pi} \rho_0), \qquad (17)$$

where $\Sigma = \int_0^\infty \rho(z) dz$ corresponds to the integrated surface density and $\rho_0 = \rho(z=0)$.

$$\rho(z) = \begin{cases} \rho_0 \exp\left\{-\left[\frac{1}{2} - \frac{z^2}{H_g^2}\left(1 - \frac{H_r}{H}\right) + \frac{|z|}{H_{\text{planet}}}\right]\right\}, & \text{if } z < H \\ \\ \rho_0 \exp\left\{-\left[\frac{1}{2}\left(\frac{z - H_r}{H_g}\right)^2 + \frac{|z|}{H_{\text{planet}}}\right]\right\} \exp\left\{-\frac{(H - H_r)H_r}{H_g}\right\}, & \text{if } z \ge H. \end{cases} \qquad (18)$$

Using Equations (17) and (16), the differential density Equation (10) can be solved to obtain $\rho(z)$: The term $|z|/H_{\text{planet}}$ in the above equation is only considered in regions where the planet gravity dominates, i.e., when $s < s_t$.

We have now all the ingredients to derive the density scale height $H$ of the disk. We follow the treatment from Hubeny (1990) by introducing the dimensionless parameters:

$$h \equiv H/H_g, \quad h_r \equiv H_r/H_g, \qquad (19)$$

which allows us to construct an algebraic equation to compute $H$ in the form of the dimensionless parameter $h$:

$$h = \frac{\sqrt{\pi}}{2} \left( \frac{h}{h - h_r} \right)^{1/2} \left[ 1 - \text{erf}\left( \{h(h - h_r)\}^{1/2} \right) \right] + \text{erf}(h - h_r) \exp[-(h - h_r)h_r], \qquad (20)$$

where $\text{erf}(x) = \int_x^\infty \exp(-t^2) dt$ is the error function.

<footnote n="9">⁹ For a detailed discussion on stellar atmospheres and the radiative transfer equations consult Mihalas (1978) and Mihalas & Weibel Mihalas (1984).</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Figure 1. Temperature profile of the disk after $10^4$ yr. The left panel shows a model without planet feedback $L_p = 0$ (temperature peak ~200 K), while the right panel includes a $L_p = 1 \times 10^{-3} L_\odot$ luminous planet (temperature peak ~1060 K).</figcaption></figure>

Once a numerical solution for $h$ is found, the final scale height is simply given by:

$$H \equiv H_{g+r} = hH_g, \qquad (21)$$

where $H_g$ is defined in Equations (11). We call $H_{g+r}$ the density scale height obtained from Equations (21) and (20), which corresponds to the scale height, containing contributions from both gas and radiation pressure. It is worth mentioning that $H_r$ in Equation (12) represents the semi-thickness of the disk only when radiative forces are dominant. When the radiation pressure is negligible, i.e., $H_r \ll H_g$, the scale height is given by:

$$H \approx H_g, \qquad (22)$$

recovering the standard calculation defined in Equation (11).

In this limit, the vertical density profile in Equation (18) has the standard form:

$$\rho(z) \approx \begin{cases} \rho_0 \exp\left\{-\frac{1}{2}(z/H_g)^2\right\}, & \text{if } s \ge s_t \\ \rho_0 \exp\left\{-\left[\frac{1}{2}(z/H_g)^2 + |z|/H_{\text{planet}}\right]\right\}, & \text{if } s < s_t, \end{cases} \qquad (23)$$

where $\rho_0 = \Sigma/(\sqrt{2\pi}H_g)$.

Summarizing, $H_{g+r}$ (Equation (21)) represents the circumstellar scale height when gas and radiation pressures are taken into account, and $H_g$ (Equation (11)) when only gas pressure is included.

### 3. Numerical Hydro-simulations

We ran a set of hydro-simulations to follow the evolution of a gaseous disk with an embedded Jupiter-mass planet, in which the planet radiates away a fraction of its internal energy at a constant rate. We modify the public FARGO-ADSG code (Baruteau & Masset 2008) following the same procedure as in Montesinos et al. (2015, 2016). Assuming hydrostatic equilibrium, we computed the vertical scale height following the expression for $H_{g+r}$ (Equation (21)) to recreate a 3D structure with an initial aspect ratio $H/r$ set to 0.05 for each model. The gravitational force of the planet uses a softening length $\epsilon = 0.6$ over which the potential is smoothed to guarantee a good

agreement between 2D and 3D simulations (Masset 2002; Baruteau & Masset 2008).

We solve a non-stationary energy equation considering a heating term due to shear viscosity, where we adopt the $\alpha$-viscosity prescription (Shakura & Sunyaev 1973) with $\alpha = 10^{-4}$. For the radiative cooling mechanism we assume blackbody emission $Q^{-} = 2\sigma T_{\text{eff}}^4$. We also take into account a radiative heating source (feedback) from the planet as in Montesinos et al. (2015).

The disk extends from 4 to 25 au, with an initial surface density given by:

$$\Sigma(r) = \Sigma_0 \frac{r_p}{r}, \qquad (24)$$

where $\Sigma_0 = 30 \text{ g cm}^{-2}$ corresponds to the density at the planet location $r_p$. The total disk mass is therefore $M_{\text{disk}} \approx 5 \times 10^{-3} M_\odot$. The embedded luminous planet ($L_p = 10^{-3} L_\odot$) is located at 10 au.

The grid resolution for all the simulations was set to $n_r = 1024$ logarithmically spaced radial sectors and $n_\theta = 1024$ azimuthal equally spaced sectors over $2\pi$. We initially put a planet located at $r_p = 10$ au, $\phi_p = 270^\circ$ on a circular orbit ($\phi = 0^\circ$ is the north). The planet is not allowed to migrate. We follow the evolution of the system in the reference frame of the planet. We present here two cases: one in which the planet feedback is turned off (i.e., $L_p = 0$), and another one in which the planet radiates at a constant luminosity $L_p = 1 \times 10^{-3} L_\odot$, which we have shown is achievable even when considering accretion feedback (Gárate et al. 2021).

Figure 1 shows the midplane temperature field after $10^4$ yr of disk evolution (or ~316 planetary orbits) at which the system has reached a quasi-steady state. The left panel shows a model without feedback ($L_p = 0$). In this case, the gas temperature at the planet location is about 200 K. When the feedback is included, the temperature peaks at 1060 K at the circumplanetary regions (right panel).

We are interested in the temperature increment at the planet location. Hence, in Figure 2 we present a slice corresponding to the last simulation orbit of the temperature profile as a function of the radial distance $r$ passing through the planet at fixed $\phi_p$. As in Figure 1, we observe that when the feedback is turned on, the temperature peaks at 1060 K, while it peaks at 200 K when the planet is not emitting.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Slice of the temperature field along the radial direction—passing through the planet located at $r_p = 10$ au—for models with and without feedback, in red and blue, respectively.</figcaption></figure>

In Figure 3 we show how the disk aspect ratio $H/r$ is affected by temperature changes due to the planet feedback. We plot $H/r$ for a slice passing over the planet, using the pure gas scale height $H_g$ model (Equation (11)) and the radiative scale height $H_{g+r}$ (Equation (21)). The left panel of Figure 3 corresponds to a model with no feedback. In this case, both calculations ($H_g$ and $H_{g+r}$) matched, indicating that radiative pressure is negligible when the feedback is turned off. In the right panel of Figure 3, we plot the model with a luminous planet. In this case, the radiative pressure locally overcomes the planet gravity, enhancing the scale height above the planet. The CPD region increased its height—with respect to $H_g$—by $(H_r - H_g)/H_g \sim (0.2 - 0.01)/0.01 = 19$. Radiative forces are dominant above the CPD when the planet feedback is activated.

# 4. Observational Effects of the Radiative Scale height

## 4.1. Vertical Optical Depth

To create a 3D distribution, we vertically extend the disk using Equation (21), where the volume density $\rho(r, \phi, z)$ is computed from Equation (18). A scale height enhancement such as the one observed in Figure 3 (right panel) could act as a bump blocking or scattering a fraction of the stellar radiation. Therefore, it is expected to cast a shadow into the outer disk in the radial direction away from the planet. However, this phenomenon only occurs if, and only if, the “bump” is optically thick in the vertical and radial directions.<fnref n="10" />

We compute the optical depth in Equation (15) from the volume density given by Equation (18), using for $\kappa$ the Rosseland mean opacity appropriate for protoplanetary disks (Semenov et al. 2003). Figure 4 shows a slice of the vertical optical depth $\tau(z, r) = \int_z^\infty \frac{1}{2}\kappa(r)\rho(r, z)dz$, passing through the planet for the model with planetary feedback. We note that, when using the radiative formulation $H_{g+r}$ (Equation (21)), an optically thick “bump” of about $\sim 2$ au of height forms above the CPD (the planet is at $r_p = 10$ au), compared with a no-bump situation if the pure gas formulation $H_g$ (Equation (11)) is used, meaning that planet gravity dominates over gas pressure alone.

As mentioned before, if the “bump” is also optically thick in the radial direction, then a shadow is expected to be cast. We explore this possibility in the next section through radiative transfer calculations.

## 4.2. Monte Carlo Radiative Transfer: Dust Vertical Structure

We perform radiative transfer calculations using the radiative transfer code RADMC3D (Dullemond et al. 2012). We feed RADMC3D with the gas distribution obtained from our hydro-simulations (see section Section 3), corresponding to an evolutionary step after $\approx 10^4$ yr. The dust of the disk is assumed to be composed of astrosilicates of intrinsic density $\rho_{\text{intr}} = 2$ g cm$^{-3}$ following a power-law size distribution $dn(a) \propto a^{-3.5}da$, with $a$ ranging from 0.1 µm to 1 cm. The dust density is normalized by imposing a gas-to-dust ratio of 100. We assume local thermo-dynamic equilibrium, i.e., $T_{\text{gas}} = T_{\text{dust}}$. The absorption efficiencies were computed using Mie theory (Bohren & Huffman 1983).

We decompose the dust size range into 12 logarithmically spaced bins representing 12 dust i-species. We obtain a dust density distribution $\Sigma_i$ for each i-bin in such a way that the sum of individual species' surface density gives the total dust density, i.e., $\Sigma_{\text{dust}}(r, \phi) = \sum_i \Sigma_i(r, \phi)$, where the sum is performed from $a_{\text{min}}$ to $a_{\text{max}}$. The Stokes number is computed from $\text{St} = \rho_{\text{intr}} a \Omega_K / \rho c_s$, where $\rho_{\text{intr}}$ is the intrinsic density of particles, $a$ the particle radius, $\Omega_K$ the Keplerian velocity, $\rho$ the volume density, and $c_s$ the isothermal sound speed.

The dust scale height is assumed to be $H_d = H_g \sqrt{\alpha/(\alpha + \text{St}_i)}$,<fnref n="11" /> where $\alpha$ is the turbulent viscosity, $\text{St}_i$ the average Stokes number of the i-th-species, and $H_g$ the density scale height computed from $H_g$ (Equation (11)) or $H_{g+r}$ (Equation (21)). This method allows us to take into account different vertical distributions for different dust species i (Youdin & Lithwick 2007).

In Figure 5, we plot the dust scale height $H_d(r, z)$ for the case with planetary feedback. The top panel corresponds to computations with $H_g$ (Equation (11)), while the bottom one used the generalized $H_{g+r}$ expression (Equation (21)). When using $H_g$ alone, the bump over the CPD disappears.

For the 3D radiative transfer calculation, we assume that the central star is of solar-type with an effective temperature of 6000 K. We incline the disk by 13°, and we set the system at 140 pc from the Earth.

In Figure 6 we show the relative intensity ($I/I_{\text{max}}$) for $\lambda = 1$ µm. We compare the effect of the feedback (switched on/off) when using the different scale height formulations ($H_g$ and $H_{g+r}$) to compute the 3D distribution as an input in RADMC3D. We remark that:

1. A shadow is cast in scattered light from the CPD projected to the outer disk only when the feedback is activated and the formulation $H_{g+r}$ is used.

2. No shadows are observed when $L_p = 0$, whether $H_g$ or $H_{g+r}$ is used to compute the synthetic image.

The features described above can be explained as follows. At the CPD location, the gravity from the planet dominates over gas pressure forces. When the feedback is activated (producing a local temperature increment from 200 to 1060 K), radiative forces play a relevant role. They can overcome the gravity from the planet and the material is puffed-up in the vertical direction

<footnote n="10">¹⁰ In this case “radial” refers to a spherical coordinate, rather than the cylindrical $r$ coordinate used in the simulations.</footnote>

<footnote n="11">¹¹ From this we have that small grains with small Stokes numbers ($\text{St} \ll \alpha$) have scale heights $\sim$gas scale height.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Radial slices of the aspect ratio $h = H/r$—passing through the planet location—computed with the pure gas formula $h_g = H_g/r$ (blue) from Equation (11), and the radiative one $h_{g+r} = H_{g+r}/r$ (red) from Equation (21). The left panel corresponds to a model without feedback. In this case both curves ($h_g$ and $h_{g+r}$) match. The middle panel corresponds to a model with planet feedback. The radiative component in the expression $h_{g+r}$ overcomes gravitational forces from the planet, locally enhancing the aspect ratio. The right panel is a zoom of the vicinity of the planet (around ~10 au) for the feedback model. The observed fluctuation around ~9.95 au corresponds to a perturbation in the temperature field due to planetary feedback.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Optical depth $\tau(r, z)$ of the disk as a function of the vertical and radial direction for a model with $L_p = 10^{-3} L_\odot$. The top panel was computed using the pure gas formulation $H_g$ (Equation (11)), the bottom panel uses $H_{g+r}$ (Equation (21)). At the planet location ($r_p = 10$ au) an optically thick column of gas (a bump) is created due to the gas pressure of the circumplanetary disk. When radiation pressure is included, the *bump* is enhanced by a factor of 19.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Grain size ($10^{-4}$–10 cm) distribution in $r, z$ for a model with $L_p = 10^{-3} L_\odot$. The slice passes through the planet located at $r_p = 10$ au. The top figure shows a model computed with $H_g$; the bottom one includes radiative pressure $H_{g+r}$. The upper layers of the disk are populated with micron-sized particles as expected, while the biggest particles (~cm) settle to the midplane.</figcaption></figure>

(Figure 3, left panel). This material is optically thick, producing a narrow shadow cast from the CPD to the outer disk regions (Figure 6).

We repeat the same calculation for $\lambda = 10 \mu$m. In this case, the shadows disappear since large dust tends to settle, making the bump optically thin at those frequencies. The shadow is observed only in reflected light and not in the dust temperature. This suggests that the “optimal” wavelength for shadow detection should be dominated by scattered light (~1–3 $\mu$m) rather than thermal emission (>10 $\mu$m). However, to validate these results, proper dust modeling is required with large grain sizes (≥100 $\mu$m) where the dust may not be well coupled to the gas as assumed in this work.

It is worth noting in Figure 6 that in models with $L_p = 10^{-3} L_\odot$ (bottom panels) some azimuthal features appear, such as wakes within the gap and the disk. Their origin comes from the planetary feedback, which heats up the gap, producing noticeable turbulence (see Montesinos et al. 2015; Gárate et al. 2021 for a discussion about the feedback).

In Figure 7, we plot the same disk emissivity reported for $\lambda = 1 \mu$m (Figure 6), but following a slice arc between 260° to

290° for different radii ranging from 16 to 20 au. The top panels of Figure 7 correspond to models without feedback, where the emissivity was calculated using $H_g$ (top left) and $H_{g+r}$ (top right). The bottom panels include the $L_p = 10^{-3} L_\odot$ feedback model with $H_g$ (left) and $H_{g+r}$ (right). Looking at the feedback model with $H_{g+r}$ (bottom right panel) an emission dip is observed up to 20 au.

## 5. Discussion

We model the evolution of a protoplanetary disk with an embedded Jupiter-mass planet in which the planet presents an intrinsic luminosity of $10^{-3} L_\odot$. When this feedback is activated, the region around the planet reaches a midplane temperature of about ~1060 K; such high luminosities are supposed to be present during the early phases of planet formation (~ Myr old planet) and can last for about ~1 Myr (Mordasini et al. 2012).

The disk's vertical structure is described by taking into account the pressure terms related to the gas and the isotropic
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Images for $\lambda = 1 \mu$m obtained from RADMC3D. Panel (a): model without feedback ($L_p = 0$), using the pure gas model $H_g$. Panel (b): model with $L_p = 0$, but using the new formula $H_{g+r}$. Panel (c): model with feedback $L_p = 1 \times 10^{-3} L_\odot$, using the pure gas formula $H_g$. Panel (d): same feedback as (c), but using $H_{g+r}$. A shadow cast from the planet location is observed in this model.</figcaption></figure>

radiation inside the fluid. The gravity of the planet is also taken into account to counterbalance the pressure gradient from the CPD (Müller et al. 2012).

The radiation pressure treatment follows an analytic expression adopted from classical stellar atmospheric pressure readjusted to disk geometry (Hubeny 1990). This treatment relies on the assumption of a vertical isothermal sound speed, which is an adequate approximation when the disk interior is optically thick (Armitage 2015).

It is worth mentioning that an equivalent but simplified calculation of the radiative scale height can also be found in Montesinos & de Freitas Pacheco (2011) where the scale height in that work considers the contribution from radiation pressure, disk self-gravity, and turbulent pressure. Results derived from

both approximations (Hubeny 1990 and Montesinos & de Freitas Pacheco 2011) applied to the present disk model are equivalent.

We show that in local thermodynamic equilibrium, the radiative pressure that arises in the vicinity of an accreting protoplanet with midplane temperatures above 1000 K can overcome the planet's gravity, producing an optically thick bump in its surroundings. For a non-accreting planet ($L_p = 0$), the radiative pressure can be neglected, as both approaches, $H_g$ (Equation (11)) and $H_{g+r}$ (Equation (21)), produce the same scale height (see Figure 3), and no shadows are observed.

In protoplanetary disks, hot regions can be reached at least in two different locations: at the inner rim of a transitional disk and a planet-forming region. The inner rim of a transitional
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Slice cuts between 260° to 290° degrees for different radial cuts—ranging from 16.1 to 20.5 au—of the emissivity reported for $\lambda = 1 \mu$m in Figure 6. The top panels correspond to models without feedback, where the emissivity was calculated using $H_g$ (top left) and $H_{g+r}$ (top right), respectively. The bottom panels include the $10^{-3}L_\odot$ feedback models with $H_g$ (left) and $H_{g+r}$ (right). The feedback model reports a strong brightness dip when $H_{g+r}$ is used.</figcaption></figure>

disk is expected to be at temperatures of about ~1500 K (Dullemond & Monnier 2010; Vinković 2014) which could result in a non-negligible contribution of the vertical radiation pressure impacting its scale height. Here, we do not attempt to describe the inner rim of disks around pre-main sequence stars. However, we suggest the use of $H_{g+r}$ (Equation (21)) to explain the extremely high aspect ratios $H/r \sim 0.2$ needed to interpret several non-resolved observations of disks (Natta et al. 2001; Olofsson et al. 2013).

When the planet feedback is included, the midplane of the CPD reaches temperatures of about ~1060K, and the inclusion of radiation pressure in the scale height equation (Equation (21)) dramatically affects the vertical disk geometry. More specifically, by surpassing the planet's gravity, radiation pressure locally enhances its height by 19% compared to a pure gas pressure model (left versus right panels in Figure 3).

From radiative transfer calculations, we showed that the “bump” created above the CPD is optically thick, casting a shadow in scattered light, which extends to the outer regions of the disk in about ~20 au from the planet. The shadows are only observed at wavelengths of ~1–3 $\mu$m, mostly because small particles (micron-size) get lifted to the upper disk layers making the bump optically thick at those wavelengths (Figure 4). At longer wavelengths, the shadows disappear, as thermal emission tends to dominate over scattered light. For close-in planets, even though they are relatively bright at 1 micron (Figure 6), their direct

detection can be harmed by speckles and other observational artifacts. On the other hand, since shadows extend much further away in the disk, they are more favorable for detection.

To conciliate our findings with observations, recent ZIMPOL/ SPHERE images obtained by Bertrang et al. (2018, 2020) revealed small-scale structures could be due to planet–disk interactions. The observations also show a moving narrow surface brightness dip with an azimuthal width of ~20°. The moving shadow seems to be cast by a large amount of optically thick material that blocks a fraction of the stellar emission, which could be interpreted as an undetected CPD (Bertrang et al. 2018). Such dips are precisely the kind of predictions obtained with our prescription for the scale height.

A compelling case is the first direct image of a forming planet within the gap of PDS 70 (Keppler et al. 2018; Müller et al. 2018), where an accreting $10 M_J$ mass planet surrounded by a protoplanetary disk was suggested (Keppler et al. 2019). Also, Christiaens et al. (2019) recently presented observational evidence of the presence of a CPD around the protoplanet PDS 70 b. Under these circumstances, a shadow should be cast by the putative CPD. In fact, some shadowing has been observed in the outer disk of PDS 70 (Long et al. 2018b). Its origin, however, remains unclear. High-contrast observations of micron-sized dust from PDS 70b could, in principle, reveal new shadows arising from the planet-forming region. This system is an exceptional laboratory to
test our ideas, where a dedicated model would be needed to predict possible shadowing effects.

Other interesting illumination patterns—non-related to accreting planets—have also been reported. For instance, in HD 142527 two diametrically opposed shadows have been discovered in polarized scattered light by Avenhaus et al. (2014). These intensity nulls were later explained by the presence of a misaligned inner disk able to block a fraction of the stellar emission (Casassus et al. 2015b; Marino et al. 2015; Price et al. 2018). The main difference with the shadows obtained in our simulations is that the latter cause both a deeper intensity null and a broader azimuthal extension (e.g., ~30° for HD 142527 as opposed to ~10° in our models). Moreover, it has been shown that these strong shadows can trigger spirals in the gas and dust distribution (Montesinos et al. 2015; Cuello et al. 2019), respectively. Also, if the inclined inner disk is precessing, then the projected shadow rotates as well (Facchini et al. 2018). Remarkably, this could create planetary-like spirals in the gaseous disk at the co-rotating region with the shadow (Montesinos & Cuello 2018). Further observations of intriguing moving shadows have been reported by Stolker et al. (2017) and Debes et al. (2017). However, there are still doubts about their origin. They could be caused by precessing inner warped disks (Nealon et al. 2019) or by moving optically thick clouds/bumps at the inner disk rims. In this work, the shadows from the CPD move at Keplerian rotation and we do not expect any impact on the gas dynamics due to such illumination effects. We emphasize that the obtained dips are not deep enough to cool down the gas and trigger azimuthal perturbations in the density field, as in Montesinos et al. (2016).

## 6. Summary

In this work we showed that for any region in a circumstellar disk with temperatures larger than ~1000 K, the vertical radiative pressure effect should be taken into account to compute the disk scale height, which can be approximated by an analytical expression (Equation (21)) derived from stellar atmospheres theory adapted to a disk geometry by Hubeny (1990).

In the vicinity of an accreting protoplanet, radiative forces surpass the planet's gravitational force, producing an optically thick bump that may cast a narrow shadow into the outer disk. Such a shadow should move at Keplerian speed with the planet. For a non-accreting planet, radiative pressure can be neglected and no shadows are observed. Due to the narrow width of the predicted shadows (~10°), observations at high angular resolution are required to reach high signal-to-noise levels. These dips should be better observed in scattered light at wavelengths close to ~1–3 µm. The shadows are expected to disappear at longer wavelengths since the bump above the CPD becomes optically thin for lower frequencies and the disk thermal emission contaminates any residual shadowing pattern.

The authors thank the referee for constructive comments and recommendations, which help to improve this paper. M.M. acknowledges financial support from the Chinese Academy of Sciences (CAS) through a CAS-CONICYT Postdoctoral Fellowship administered by the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. N.C. acknowledges financial support provided by FONDECYT grant 3170680. G.H.M.B. acknowledges funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No.

757957). J.O. acknowledges financial support from Fondecyt (grant 1180395). A.B. acknowledges support from FONDECYT grant 1190748. C.P. acknowledges financial support from FONDECYT (grant 3190691). M.M., N.C., J.C., J.O., A.B., and C.P. acknowledge support from Iniciativa Científica Milenio via the Núcleo Milenio de Formación Planetaria. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 210021.

## ORCID iDs

Matías Montesinos https://orcid.org/0000-0001-9789-5098

Johan Olofsson https://orcid.org/0000-0003-4475-3605

Amelia Bayo https://orcid.org/0000-0001-7868-7031

Gesa H.-M. Bertrang https://orcid.org/0000-0001-5127-0172

## References

* ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJL, 808, L3

* Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

* Armitage, P. J. 2015, arXiv:1509.06382

* Avenhaus, H., Quanz, S. P., Schmid, H. M., et al. 2014, ApJ, 781, 87

* Baruteau, C., & Masset, F. 2008, ApJ, 678, 483

* Bertrang, G. H.-M., Avenhaus, H., Casassus, S., et al. 2018, MNRAS, 474, 5105

* Bertrang, G. H. M., Flock, M., Keppler, M., et al. 2020, arXiv:2007.11565

* Bohren, C. F., & Huffman, D. R. 1983, Absorption and Scattering of Light by Small Particles (New York: Wiley)

* Casassus, S., Marino, S., Pérez, S., et al. 2015a, ApJ, 811, 92

* Casassus, S., Wright, C. M., Marino, S., et al. 2015b, ApJ, 812, 126

* Christiaens, V., Cantalloube, F., Casassus, S., et al. 2019, ApJL, 877, L33

* Cuello, N., Montesinos, M., Stammler, S. M., Louvet, F., & Cuadra, J. 2019, A&A, 622, A43

* Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205

* Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957

* Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, ascl:1202.015

* Dullemond, C. P., & Monnier, J. D. 2010, ARA&A, 48, 205

* Facchini, S., Juhász, A., & Lodato, G. 2018, MNRAS, 473, 4459

* Gárate, M., Cuadra, J., Montesinos, M., & Arévalo, P. 2021, MNRAS, 501, 3113

* Garufi, A., Quanz, S. P., Schmid, H. M., et al. 2014, A&A, 568, A40

* Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018a, ApJL, 869, L42

* Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018b, ApJL, 869, L43

* Hubeny, I. 1990, ApJ, 351, 632

* Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

* Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118

* Long, F., Pinilla, P., Herczeg, G. J., et al. 2018a, ApJ, 869, 17

* Long, Z. C., Akiyama, E., Sitko, M., et al. 2018b, ApJ, 858, 112

* Marino, S., Perez, S., & Casassus, S. 2015, ApJL, 798, L44

* Masset, F. S. 2002, A&A, 387, 605

* Meeus, G., Waters, L. B. F. M., Bouwman, J., et al. 2001, A&A, 365, 476

* Mihalas, D. 1978, Stellar Atmospheres (San Francisco, CA: W. H. Freeman)

* Mihalas, D., & Weibel Mihalas, B. 1984, Foundations of Radiation Hydrodynamics (New York: Oxford Univ. Press)

* Montesinos, M., Cuadra, J., Perez, S., Baruteau, C., & Casassus, S. 2015, ApJ, 806, 253

* Montesinos, M., & Cuello, N. 2018, MNRAS, 475, L35

* Montesinos, M., & de Freitas Pacheco, J. A. 2011, A&A, 526, A146

* Montesinos, M., Perez, S., Casassus, S., et al. 2016, ApJL, 823, L8

* Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012, A&A, 547, A111

* Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2

* Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123

* Natta, A., Prusti, T., Neri, R., et al. 2001, A&A, 371, 186

* Nealon, R., Pinte, C., Alexander, R., Mentiplay, D., & Dipierro, G. 2019, MNRAS, 484, 4951

* Olofsson, J., Benisty, M., Le Bouquin, J.-B., et al. 2013, A&A, 552, A4

* Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJL, 890, L9

* Pinte, C., van der Plas, G., Ménard, F., et al. 2019, NatAs, 3, 1109
Price, D. J., Cuello, N., Pinte, C., et al. 2018, MNRAS, 477, 1270

Quanz, S. P., Avenhaus, H., Buenzli, E., et al. 2013, ApJL, 766, L2

Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A,
410, 611

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174

Stolker, T., Dominik, C., Avenhaus, H., et al. 2016, A&A, 595, A113

Stolker, T., Sitko, M., Lazareff, B., et al. 2017, ApJ, 849, 143

Vinković, D. 2014, A&A, 566, A117

Youdin, A. N., & Lithwick, Y. 2007, Icar, 192, 588

Zhu, Z., Stone, J. M., Dong, R., Rafikov, R., & Bai, X.-N. 2015, AAS/
Division for Extreme Solar Systems Abstracts, 3, 105.09