© ESO 2020

1 Introduction

Planets are born in protoplanetary disks. While they do not always make their presence clear like in the case of PDS 70b and c (Keppler et al. 2018; Haffert et al. 2019), the ALMA observations of the DSHARP survey (Andrews et al. 2018) have provided theorists with high-fidelity datasets for testing the planet–disk interaction theory. It is fascinating how well a single planet can reproduce the annular substructures in some of the observed systems, for instance, Zhang et al. (2018).

Most of the aforementioned systems show resolved features at the 100 au scale, suggesting that the disk temperature profile is set by a balance between heating by stellar irradiation and thermal cooling. This makes modeling these systems with a locally isothermal equation of state quite attractive because the cooling timescale in the optically thin irradiation-dominated outer disk is commonly sufficiently short (shorter than a hundredth of an orbit) to render radiative effects negligible. At the same time, the very low level of effective viscosity that is inferred for these disks (Zhang et al. 2018) indicates minuscule contributions by viscous heating to the thermal budget of the disk, such that the locally isothermal assumption seems further justified.

However, an embedded planet also interacts with the disk gravitationally, forming spiral arms (Ogilvie & Lubow 2002; Rafikov 2002) that permeate the disk and can steepen into shocks (Zhu et al. 2015). Several studies have shown that heating by these spiral shocks is another significant heat source at distances of several au for solar-type stars (Richert et al. 2015; Lyra et al. 2016; Rafikov 2016; Ziampras et al. 2020). While the planet’s contribution to the energy budget decreases with increasing distances from the star, there is a gray area at a few tens of au where shock heating could operate to some degree, even when the cooling timescale is a small fraction of the local orbital period.

In addition, Miranda & Rafikov (2019) showed that even in seemingly locally isothermal scenarios, assuming an adiabatic equation of state results in fundamentally different physics concerning the angular momentum flux in protoplanetary disks, and it therefore leads to noticeable changes in dust continuum profiles. They quoted distances of about 80 au as the lower limit beyond which a locally isothermal equation of state can indeed be justified, and showed that deviations of about 10⁻³ from γ = 1 can lead to noticeable differences (of about 10%) in simulated outcomes and therefore observables such as continuum emission intensity.

This subject is quickly gaining clarity as more effort is made to understand planet-induced gap opening in order to interpret observations and constrain the properties of protoplanetary disks. During the reviewing process for this paper, two publications by Miranda & Rafikov (2020) and Zhang & Zhu (2020) were made available, which discuss the effect of a finite cooling timescale on the location and number of gaps that can be produced by a single planet. Our study does not reiterate these new results, but rather enriches them by handling radiative effects differently. It serves as an additional point of view in arriving at a similar conclusion: radiative cooling is of key importance in determining the radial structure of a disk with an embedded planet.

We use numerical simulations to address the importance of proper treatment of radiative effects for two systems that mimic the properties of HD 163296 and AS 209 (Andrews et al. 2018; Zhang et al. 2018) and show that an adiabatic equation of state can producedifferent results with respect to a locally isothermal model when the gap-opening capabilities of a planet are modeled.

In Sect. 2 we present our physical and numerical setup in modeling these two systems. We present our results in terms of disk structures in Sect. 3, and discuss their implications in Sect. 4. Finally, we summarize our results and conclude this study in Sect. 5.

Fig. 1 ALMA continuum emission observations of the two systems we used as testbeds for our numerical models, as shown in Huang et al. (2018). Solid and dotted white arcs mark the location of bright rings and dark gaps, respectively, along with their distance from the host star in au. It is likely that a planet is responsible for one or more of the gaps observed in either system. Left: HD 163296, a 12.6 Myr old likely isolated system that features at least two clearly visible gaps at 48 and 86 au. Our models focus on the 20–60 au range, with a planet carving a gap at 48 au. Right: AS 209, a 1 Myr old system in the Ophiuchus region that is rich in annular structures, boasting five rings between 20 and 130 au. This source has been modeled by Zhang et al. (2018), who found that a fit with a single planet at 99 au could reconstruct all five of these rings.

2 Sources, physics, and numerics

In this section we describe the physical and numerical modeling of HD 163296 and AS 209. We first present the two systems that we use as our reference. We then describe the source terms in the vertically integrated hydrodynamics equations, our physical assumptions in terms of star, planet, and disk properties, and the initial and boundary conditions for our simulations.

2.1 Sources: HD 163296 and AS 209

Rather than carrying out a study on an arbitrary toy model of a protoplanetary disk, we decided to model two of the sources targeted by the DSHARP survey (Andrews et al. 2018), namely HD 163296 and AS 209 (see Fig. 1). Both of these sources feature annular structures in the form of bright rings and dark gaps in dust continuum emission images, with the former also showing a crescent at roughly 50 au and the latter standing out with its ring-rich emission profile. These systems have been modeled in the past and their various features have been studied (Huang et al. 2018; Dullemond et al. 2018; Fedele et al. 2018; Zhang et al. 2018; Zhang & Zhu 2020), and using them as testbeds for our study can provide a deeper insight on any planets they might harbor.

2.2 Hydrodynamics, heating, and cooling

The two-dimensional (2D) vertically integrated Navier–Stokes equations in cylindrical {r, ϕ, z} coordinates for a perfect gas with surface density Σ, velocity v and specific midplane internal energy e read $$\begin{array}{ll}\frac{\text{d}\Sigma }{\text{d}t}\hfill & =-\Sigma \nabla \cdot v,\hfill \end{array}$$(1a) $$\begin{array}{ll}\Sigma \frac{\text{d}v}{\text{d}t}\hfill & =-\nabla p+\Sigma g+\nabla \cdot \Sigma ,\hfill \end{array}$$(1b) $$\begin{array}{ll}\frac{\text{d}\Sigma e}{\text{d}t}\hfill & =-\gamma \Sigma e\nabla \cdot v+Q_{\text{visc}}+Q_{\text{irr}}-Q_{\text{cool}},\hfill \end{array}$$(1c)

where γ is the adiabatic index and p = (γ − 1)Σe is the vertically integrated pressure. External source terms (in our case, gravitational forces) are contained in g, and Σ denotes the viscous stress tensor (Tassoul 1978).

For a thin disk of gas rotating on a Keplerian orbit with orbital frequency ${\Omega }_{\text{K}}=\sqrt{GM/r^3}$, at distance r around a star of mass M, the pressure scale height of the gas is $H=c_{\text{s}}{}_{\text{iso}}/{\Omega }_{\text{K}}$. Here, $c_{\text{s}}{}_{\text{iso}}=\sqrt{p/\Sigma }$ is the isothermal sound speed and relates to the adiabatic sound speed c_s as $c_{\text{s}}=\sqrt{\gamma }{c}_{\text{s}}{}_{\text{iso}}$. The aspect ratio of the disk is then h = H∕r. The gravitational constant and the mean molecular weight of the gas are denoted by G and μ, respectively. We adopted a typical solar-composition disk of molecular H–He gas, therefore γ = 7∕5 and μ = 2.35.

In a locally isothermal framework, the energy equation in Eq. (1) is not solved. Instead, a fixed sound speed profile is adopted. This in turn defines a (fixed) radial temperature profile T, as $c_{\text{s}}{}_{\text{iso}}=\sqrt{R_{\text{gas}}T/\mu }$, with R_gas being the gas constant. In the locally isothermal scenario this temperature profile is constant on cylinders and only dependent on radius, but we focus on the midplane temperature in this 2D approximation.

More generally, however, the full system of equations is solved. In this case, we include viscous heating, stellar irradiation, and thermal cooling in the energy equation as follows: $$\begin{array}{ll}{Q}_{\text{visc}}\hfill & =\frac{1}{2\nu \Sigma }\text{Tr}({\Sigma }^2)=\frac{1}{2\nu \Sigma }\left({\sigma }_{rr}^2+2{\sigma }_{r\varphi }^2+{\sigma }_{\varphi \varphi }^2+{\sigma }_{zz}^2\right)\hfill \end{array}$$(2a) $$\begin{array}{ll}{Q}_{\text{irr}}\hfill & =2\frac{L_{*}}{4\pi r^2}(1-\varepsilon )\left(\frac{\text{d}\log H}{\text{d}\log r}-1\right)h\frac{1}{{\tau }_{\text{eff}}},\hfill \end{array}$$(2b) $$\begin{array}{ll}{Q}_{\text{cool}}\hfill & =2{\sigma }_{\text{SB}}\frac{T^4}{{\tau }_{\text{eff}}},\hfill \end{array}$$(2c)

where ν is the kinematic viscosity. For this we used the α-ansatz (Shakura & Sunyaev 1973), which for thin disks reads ν = αc_sH. The Stefan-Boltzmann constant is denoted by σ_SB, and τ_eff is an effective optical depth following Hubeny (1990) $${\tau }_{\text{eff}}=\frac{3\tau }{8}+\frac{\sqrt{3}}{4}+\frac{1}{4\tau },$$(3)

where τ is the vertical optical depth measured from the midplane to the disk surface, $$\tau ={\displaystyle {\int }_0^{\infty }\kappa }\rho \text{d}z\approx c_1{\kappa }_{\text{mid}}{\rho }_{\text{mid}}H.$$(4)

Here, ρ is the volume density and κ_mid(ρ, T) is the Rosseland mean opacity at the disk midplane, which we refer to as simply κ and evaluate in Sect. 2.3 below, either using a constant value or adopting the temperature-dependent opacity law. To match the opacity drop with height, we included a correction factor c₁ = 1∕2 (Müller & Kley 2012). The midplane gas density ρ_mid is related to the surface density such that $\Sigma =\sqrt{2\pi }{\rho }_{\text{mid}}H$ (i.e., corresponding to a Gaussian vertical stratification at hydrostatic equilibrium).

Our cooling and irradiation prescription adapts the model by Menou & Goodman (2004) for a disk with an albedo of ɛ = 1∕2 around a star with luminosity L_*. The factor $\frac{\text{d}\log H}{\text{d}\log r}$ was assumed to be constant and equal to 9∕7 (i.e., disk self-shadowing was not considered). For more details, see the physical setup by Ziampras et al. (2020).

2.3 Numerics

We used the PLUTO code (Mignone et al. 2007) along with the FARGO method (Masset 2000; Mignone et al. 2012) for our simulations. Models with an embedded planet were run on a polar {r, ϕ} grid, logarithmically spaced in the radial direction.

In order to set our initial conditions, we constructed two different disks that mimic the structure of HD 163296 and AS 209. We therefore set M_* and L_* according to Andrews et al. (2018), and decided to embed a single planet with mass M_p = 0.5M_J at r_p = 48 au for HD 163296 and M_p = 0.083M_J at r_p = 99 au for AS 209. We adopted a viscous α following Zhang et al. (2018) (either 10⁻⁴ or 10⁻⁵, constant throughout the disk). For the Rosseland mean opacity we chose either a constant value κ = 0.45 cm² g⁻¹ (adapted from Birnstiel et al. 2018) or a temperature-dependent opacity law by Lin & Papaloizou (1985). The latter dictates that κ = 5 × 10⁻⁴T² for temperatures under 170 K, which is true for the full extent of the simulated disk. A list of the physical and numerical parameters used in our models is given in Table 1.

The initial surface density needs to be prescribed, and the temperature is simply given by Q_visc + Q_irr = Q_cool (see Eq. (1c)) in the absence of a planet (we can ignore compression heating in this scenario because its contribution is negligible). We adopted values inferred in Table 3 by Zhang et al. (2018) for HD 163296 such that the surface density at r_gap = {10, 48, 86} au is Σ(r_gap) = {100, 10, 3} g cm⁻², respectively.We then fit a power law to these three points and find that a fit of $$\Sigma {(r)}^{\text{HD}}\approx 10{\left(\frac{r}{r_{\text{p}}}\right)}^{-3/2}\frac{\text{g}}{{\text{cm}}^2}\approx 3200{\left(\frac{r}{\text{au}}\right)}^{-3/2}\frac{\text{g}}{{\text{cm}}^2}$$(5)

matches very well. For AS 209, we simply adopt a setup similar to that by Zhang et al. (2018), such that $$\Sigma {(r)}^{\text{AS}}\approx 10{\left(\frac{r}{r_{\text{p}}}\right)}^{-1}\frac{\text{g}}{{\text{cm}}^2}.$$(6)

As far as temperature is concerned, the effect of viscous heating Q_visc is functionally negligible within our simulation range for our α value. In the absence of a planet, the temperature profile is therefore set by a balance between irradiation heating and thermal cooling. In Eqs. (2b) and (2c), the optical depths for absorption and emission in Q_irr and Q_cool, respectively, can in principle be different, but as a first approximation we assumed that they are identical and can be factored out when we set Q_irr = Q_cool. The temperature profile then only depends on stellar properties so that $$T(r)=T_{\text{p}}{\left(\frac{r}{r_{\text{p}}}\right)}^{-3/7}\text{K}\propto \frac{h^2}{r}\Rightarrow h(r)=h_{\text{p}}{\left(\frac{r}{r_{\text{p}}}\right)}^{2/7}$$(7)

for our irradiation prescription. By replacing the geometrical factor $\left(\frac{\text{d}\log H}{\text{d}\log r}-1\right)h$ with a constant irradiation angle ϕ (e.g., Dullemond et al. 2018), we recover the familiar formula where T ∝ r^−1∕2 and h ∝ r^1∕4 instead (notused here).

After constructing our initial surface density and temperature profiles, we first ran a one-dimensional radiative simulation without a planet and verified that the profiles indeed correspond to an equilibrium state. We then embedded the planet, and using the power-law profiles given in Eqs. (5)–(7) as initial conditions, executed various radiative as well as locally isothermal simulations.

In all setups, Σ is damped to the initial profiles at the boundaries according to de Val-Borro et al. (2006) over a timescale of 0.3 boundary orbital periods. The planet does not accrete material or migrate through the disk, and traces circular orbits around the star. The target simulation time was 1000 orbits for HD 163296 – which at 48 au translates into roughly 230 kyr – and 5000 orbits for AS 209 (5.4 Myr at 99 au) because these models have a very low viscosity of 10⁻⁵. We used a resolution of 573 × 1118 and 673 × 1312 cells for HD 163296 and AS 209, respectively, which results in square cells with 10–12 cells per scale height in the radial direction at the planet’s location.

Similar to Ziampras et al. (2020), the planet’s presence is described by an additional gravitational force. We accounted for the shift of the system barycenter that is caused by the planet, but neglected backreaction of the disk onto the star and planet.

To prevent singularities near the planet, we introduced a softening length ɛ = 0.6H following Müller et al. (2012), using the local pressure scale height and a correction factor that accounts for the disk thickness.

3 Results

In this section we present the results of our numerical simulations. As stated above, we tested several different models. We focus in this section on the locally isothermal and radiative models with parameters described in the previous section. The remaining models are used for comparison purposes and are discussed in Appendix A.

First, we present our results on HD 163296. We start with focusing on overall disk profiles and differences between the locally isothermal and radiative models in terms of the radial distribution of gas in the disk. We then compare and contrast the structure of spiral arms between the two models in terms of surface density and temperature contrast with respect to the disk background. This is then followed up by results on AS 209 regarding the first point only because this is the most crucial to the discussion.

Fig. 2 Azimuthally averaged surface density, midplane pressure, and aspect ratio (used as a proxy for temperature, see Eq. (7)) of HD 163296 for a locally isothermal (in orange) and a radiative (in green) disk after 1000 orbits with an embedded planet with mass 0.5 MJ at 48 au. Vertical black lines mark the location of the planet. The dashed and dotted red lines in the bottom panel denote the aspect ratio profiles used in the simulations by Zhang et al. (2018), where h ∝ r1∕4 and hp = 0.05 and 0.07, respectively. The fact that our aspect ratio profile is bounded by these two curves allows us to validate our results by comparing them to that study (see Sect. 4).

3.1 HD 163296

The first system we studied is HD 163296. This system features several rings, as shown in Fig. 1.

3.1.1 Disk profiles and secondary gaps

Our main results are plotted in Fig. 2. By comparing our radiative and isothermal runs and computing the aspect ratio as a proxy for temperature(because $h\propto \sqrt{rT}$), we find thatthe azimuthally averaged temperature profiles of the two are practically identical with each other (bottom panel in Fig. 2), and therefore with the profile that an irradiation-dominated disk should have. However, the azimuthally averaged surface density profiles of the two models differ in the inner disk and around the planet within 1000 orbits. The difference lies in the additional depression of the gas surface density at around 28 au for the locally isothermal case, in contrast to the smoother profile in the radiative run. This secondary gap, while shallow, can warrant the existence of a pressure bump in the inner disk, creating a dust trap and therefore a second and potentially visible dust ring. In contrast to the isothermal run, this pressure minimum disappears in the radiative simulation (middle panel in Fig. 2).

A rather small (compared to that at 28 au) but visible difference can also be seen within the gap region. This may be related to nonlineareffects caused by the planet’s gap opening, and discuss this in Sect. 4. Finally, the azimuthally averaged surface density in the outer disk is indistinguishable between the two models.

In Fig. 3 we compare the surface density in the disk to its initial profile. A gap and ring are clearly visible at 28 and 34 au, respectively, in the locally isothermal model. Additionally, spiral arms in the inner disk have slightly higher contrast with the background for that model compared to its radiative counterpart. We carry out a more detailed comparison in the next section.

Fig. 3 Surface density perturbation with respect to the initial profiles for the two models of HD 163296. The locally isothermal run features a secondary gap at roughly 28 au and slightly more prominent inner spirals.

3.1.2 Spiral arm contrast

As mentioned in the work of Miranda & Rafikov (2019), the locally isothermal equation of state overestimates the contrast of structures in protoplanetary disks. We already saw this behavior in Fig. 2, where a secondary gap is visible at 28 au. However, nonaxisymmetric features in the disk such as spiral arms are lost when averaging along annuli. With this in mind, it is useful to track the spiral arms that are excited by the planet and compare the surface density along their crests with its azimuthally averaged values. In doing so, we can estimate their contrast with respect to the disk background.

This comparison was carried out by tracking the trajectory of each individual spiral, logging the surface density Σ_arm along their peaks and then plotting the ratio ${\Sigma }_{\text{arm}}/\overline{\Sigma }$, where $\overline{\Sigma }$ is the azimuthally averaged surface density (see Ziampras et al. 2020 for details). Our results are summarized in Fig. 4 and show that spirals in the locally isothermal model are indeed consistently “stronger” (i.e., have ahigher contrast with respect to the disk), with the effect being more prominent in the inner disk (top panel).

The pitch (or opening) angle β of these spirals can be defined (e.g., Zhu et al. 2015) through $$\tan \beta \equiv \frac{\text{d}r}{r\text{d}\varphi }\approxeq \frac{c_{\text{s}}}{r|\Omega (r)-{\Omega }_{\text{p}}|}=\frac{h}{|1-{(r/r_{\text{p}})}^{3/2}|}.$$(8)

Equation (8) suggests that the pitch angle of spirals should be a function of the aspect ratio (again, a proxy for the temperature) for a given distance r. With this inmind, and given that the two models share an identical radial temperature profile, we expect and observe that the pitch-angle profiles of both primary and secondary spirals match very well for the locally isothermal and radiative model (lower panel).

We also compared the temperature along spiral crests with the initial axisymmetric profile for the radiative model. We found that differences are about 10–15% or smaller, and that spiral heating on the entire disk is negligible (see Fig. 5). This is to be expected in the optically thin irradiation-dominated region of the disk at r ≳ 20 au (e.g., Rafikov 2016).

Fig. 4 Comparison of spiral arm properties between our two models for HD 163296. Top: primary and secondary spiral arm contrast with respect to the disk background. We find that the locally isothermal model shows consistently higher contrast. Bottom: pitch angle of spirals as a function of radius. Because the two models have identical aspect ratios, the overlap of the two curves is to be expected. The vertical black line marks the planet location. The dotted black line corresponds to the analytical formulae in Rafikov (2002), Eq. (44), and Muto et al. (2012), Eq. (1).

Fig. 5 Temperature contrast of spirals with respect to the disk background for the radiative run of HD 163296. The effect of spiral heating on the global disk is negligible, and the temperature along spiral crests can increase by up to 15%.

3.2 AS 209

The goal in modeling this additional source was to see how this disparity between the locally isothermal and adiabatic equations of state would affect the gas surface density profile of a system with more gaps observable in continuum emission: five in AS 209 (at r ∈ {24, 35, 61, 90, 105} au) compared to only one in HD 163296 (at r = 48 au) for our simulation domain (Huang et al. 2018). While reproducing all five gaps might be difficult with an aspect ratio h_p ≈ 0.08, we see substantial differences between the two setups not only in the number of producible gaps, but also in the shape of the primary gap around the planetary orbit.

Our results are plotted in Figs. 6 and 7 and match our findings for HD 163296 in that the radiative model does not agree with the locally isothermal one. A single planet at 99 au reproduces at least three gaps at 47, 84 and 115 au in the locally isothermal model (with a possible fourth gap at 22 au, which is however subject to boundary effects), whereas the radiative model produces a very smooth monotonic profile with no gaps except for one exactly at the planet location (Fig. 6). In addition, the primary gap structure is entirely different between the two models, with the locally isothermal model showing a wider gap region that contains more material at the orbital radius of the planet, however. This can be seen in both figures, and is discussed in detail in Sect. 4.

Fig. 6 Surface density profile for our numerical model of AS 209 for a locally isothermal and a radiative equation of state after 5000 orbits. The differences in the inner disk as well as the gap region are clear between the two models. Vertical dotted lines mark location candidates for dust gaps based on the formula determined by Zhang et al. (2018), and the dashed line marks the planet location.

Fig. 7 Similar to Fig. 3, but for AS 209. The sharper annular and angular structures in the locally isothermal model are clearly visible.

4 Discussion

In Sect. 3 we showed that modeling a protoplanetary disk with a locally isothermal equation of state can under certain circumstances prove incorrect because this assumption might affect the planet–disk interaction process even when the temperatureprofile is largely unaffected by the choice of equation of state. By comparing against a radiative simulation where viscous and irradiation heating and thermal cooling are included, we showed that radiative effects play an important role in the disk evolution even at large distances from the source star.

Our testbeds aimed at reproducing the inner 50 au structure of the system HD 163296 and the gap-rich structure of AS 209. At distancesof 48 au for HD 163296 and 99 au for AS 209 (which is where planets are suspected to be), the disk is strongly dominated by irradiation, whereas viscous and spiral heating effects are negligible. Even so, the cooling timescale atthe range of 20–100 au in these systems is still significant enough to invalidate the locally isothermal assumption.

We can estimate the cooling timescale as $$\frac{\partial \Sigma e}{\partial t}~\frac{\Sigma e}{t_{\text{cool}}}~Q_{\text{cool}}\Rightarrow t_{\text{cool}}\propto \frac{{\tau }_{\text{eff}}\Sigma }{T^3}.$$(9)

We then compare t_cool to the local orbital period P_orb = 2π∕Ω_K in Fig. 8. We find that t_cool corresponds to 2–20% of P_orb in the 20–50 au range in the absence of a planet for HD 163296, and 1–10% in the 40–100 au range for AS 209. These fractions are clearly significant in this context, implying that the angular momentum flux driven by the spiral waves of a planet is described by an adiabatic framework (Miranda & Rafikov 2019).

At the same time, we see that the cooling timescale drops below 0.1% of the orbital period at distances of ~100 au for HD 163296, or ~200 au for AS 209 (see Fig. 8). At these distances we find that the locally isothermal model approaches the radiative one both in terms of radial disk structure and spiral arm contrast, suggesting that the locally isothermal assumption can still be justified when the cooling timescale is short enough. This is consistent with the latest results by Miranda & Rafikov (2020), who suggested that a cooling timescale of about 10⁻³–10⁻² orbital periods is sufficiently short to match locally isothermal disks for low- and high-mass planets, respectively. Regarding this last statement, the different gap structure between equations of state around the planetary orbit in AS 209 can be explained by comparing the planet mass to the thermal mass of the disk (Zhu et al. 2015) $$M_{\text{th}}\equiv \frac{c_{\text{s}}^3}{\text{G}{\Omega }_{\text{p}}}=h_{\text{p}}^3{M}_{*}\approx 1M_J{\left(\frac{h_{\text{p}}}{0.1}\right)}^3\left(\frac{M_{*}}{M_{\odot }}\right).$$(10)

Using this formula, we find that a 0.5 M_J planet in HD 163296 corresponds to 1.3 M_th, and a 0.083 M_J planet in AS 209 to 0.18 M_th. This means that the former system is prone to nonlinear effects due to gap opening, and therefore the surface density profiles around the primary gap are largely the same between the two different equations of state. On the other hand, our numerical model of AS 209 hosts a planet that falls in the linear regime, and therefore a cooling time of about 1% of the orbital period is not short enough to allow the radiative and locally isothermal models to agree with each other. Instead, the two profiles show a substantially different structure in the region that corotates with the planet, and they only overlap at the outskirts of the disk (where planet-driven effects are weak and/or damped anyway).

To verify that the locally isothermal profile can indeed be recovered with a more efficient cooling prescription, we executed several additional radiative simulations of HD 163296 where we limited adiabatic effects (γ = 1.01) or artificially reduced the cooling timescale (using κ = 0.045 cm² g⁻¹), thereby limiting radiative effects. These results are shown in Appendix A.

While this exercise shows that the inner disk is prone to radiative effects in our simulations, this result should be taken with a grain of salt when comparing to observations. Both the optical depth and surface density (which depend on the opacity and initial conditions) are uncertain and not set in stone. Additionally, if a substantial part of the dust grains has grown into larger grains, the opacity (and in turn, the optical depth) could be much lower, reducing the cooling timescale.

Our results on HD 163296 can also be compared to the 2D locally isothermal simulations by Zhang et al. (2018). In Fig. 2 of their study, they show that a planet with M_p = 0.3 M_J (q = 2.9 × 10⁻⁴, in our case q = 2.3 × 10⁻⁴) can open a secondary gap at roughly 0.6 r_p (28 au for our single-planet model of HD 163296) when the aspect ratio h_p is between 0.05 and 0.07 and scales with r^1∕4. Because our computed aspect ratio lies within the profiles generated using those values (see Fig. 2), we expect exactly (or at least) one secondary gap in our locally isothermal simulation. The fact that we indeed observe it suggests that our results agree with Zhang et al. (2018) in the locally isothermal limit. In addition, we observe the secondary gap at the location inferred by the fitting formula suggested in that study. However, we used a steeper surface density profile than they did (Σ₀ ∝ r^−3∕2 as opposed to r⁻¹). We carried out two more simulations with a shallower surface density profile and found similar results (see Appendix A).

By taking into account the above points, we see that the secondary-gap-opening capabilities of a single planet are exaggerated in a locally isothermal disk, such that a secondary gap in gas surface density can form within that framework but not always when radiative effects are properly treated. In the case of HD 163296 this could constrain the time of planet formation because a ring at 34 au is indeed not visible (Huang et al. 2018). This assumes that our estimates of α among other model parameters are viable.

Our simulations do not include a dust component, meaning that we cannot verify whether this secondary gap generated in the locally isothermal simulation is visible. Nevertheless, we expect that large grains can be trapped in the pressure bump formed at 34 au, forming a ring and therefore increasing the contrast between the inner disk and the secondary gap region. A similar argument can be used for the crescent structure within the gap that can be seen in the ALMA observationsof HD 163296. We find that the perturbed gas surface density at the trailing Lagrange point is higher than that at the leading one (see Fig. 9), and we suggest that dust–gas interaction could collect grains at the L4 point. This feature, however, should be transient as the gap region continues to empty as the disk evolves.

The system HD 163296 shows more structure between 50 and 200 au, namely a clear second gap in continuum emission at 86 au and a third, slightly less visible gap at 145 au. Assuming the existence of planets at all three locations, we can justify our results as far as annular structures in the 10–50 au range are concerned because the planet at 48 au will shield the inner disk by opening a gap and therefore halt the propagation of spiral waves by the outer planet(s).

Nevertheless, it would be interesting to include more planets and model the 10–200 au range of this system, assuming that each gap in dust continuum stems from a single planet opening a corresponding gap in gas density. This is further motivated by the kinematic detection of Jupiter-sized planets at 83 and 137 au by Teague et al. (2018). We therefore carried out a simulation with three planets at 48, 86, and 145 au, respectively, and find similar results in the 10–50 au range. These results are shown and discussed in Appendix B.

There are several other sources for which ALMA has provided high-fidelity observational datasets. As far as ring structures are concerned, one system stands out: AS 209 (Andrews et al. 2018) shows at least five rings with as many gaps in dust continuum emission (Huang et al. 2018). Numerical models strongly suggest that one or more growing planets are responsible for their formation (e.g., Fedele et al. 2018).

It has been suggested that a single planet at 99 au might be able to carve most of these gaps (Zhang et al. 2018). We therefore found it useful to examine the importance of radiative effects on a system with such a rich annular structure, and carried out a comparison similar to that for HD 163296 against a locally isothermal simulation, shown in Sect. 3.2. Strikingly, while our locally isothermal model produces several rings throughout the disk, we find that this is not the case in the radiative model, which shows no ring structure except for the unavoidable pressure bump in the outer disk that is caused by the clearing of the planet’s corotating region. Our results suggest that we would need multiple planets to explain the ring–gap structure in such systems.

Fig. 8 Azimuthally averaged cooling timescale in units of the local orbital period of the disk after 1000 and 5000 planetary orbits for HD 163296 and AS 209, respectively. Differences between locally isothermal and radiative models become clear in the inner disk, where the cooling timescale corresponds to a non-negligible fraction of the orbital period (roughly 1–10% near the planet, and comparable to the orbital period farther in). In the outer disk, cooling is sufficiently fast so that the two models overlap as far as azimuthally averaged profiles are concerned.

Fig. 9 Azimuthal cut at r = rp = 48 au. The surface density trailing behind the planet (left side, ϕ < 0) is slightly higher than that at the leading Lagrange point (ϕ > 0). The dashed black lines mark the planet’s hill radius and contain 20 cells. The dotted lines mark the location of the L4 (left) and L5 (right) Lagrange points.

5 Conclusions

Our goal was to examine the importance of radiative effects in disk evolution with regard to the gap-opening capabilities of a single planet. To this end, we carried out 2D numerical simulations of planets embedded in disks that resemble two systems that have recently been imaged in high angular resolution: HD 163296 and AS 209. We then compared locally isothermal simulations, where a fixed radial sound speed profile is prescribed, to setups where radiative effects are self-consistently treated.

We found that locally isothermal models exaggerate the contrast of planet-generated features with respect to radiative models and therefore overestimate the planet’s ability to carve a secondary gap in its disk. This is consistent with previous results and implies that a single planet cannot always explain the existence of multiple gaps.

We also found that the contrast of spiral arms launched by a planet is artificially sharpened within a locally isothermal framework. While this phenomenon is weak or negligible in the optically thinner outer disk, it becomes more significant for inner spirals. Regardless, spiral shock heating is negligible, with spirals having a low temperature contrast with the background and the disk being sufficiently optically thin.

Finally, by running a simulation of HD 163296 over an extended range that contained three planets, we found that our results in the limited range of 10–60 au (i.e., inside and around the innermost planet) remain unchanged. We also showed that an interplay between planet mass and the cooling timescale can lead to slight differences between a locally isothermal and an adiabatic equation of state even at 145 au, although viscous evolution might render such differences negligible.

In conclusion, the locally isothermal assumption proves to be dangerous even at the range of tens of au regarding planet–disk interaction and should therefore be in general avoided in favor of an adiabatic equation of state with a prescription for radiative cooling in the disk. By estimating the cooling timescale t_cool, the usage of such an assumption in a regime where cooling occurs very rapidly compared to the orbital period P_orb might be justified. This corresponds to t_cool at least shorter than 1% of P_orb for massive planets, or 0.1% for low-mass planets, in good agreement with Miranda & Rafikov (2020).

Acknowledgements

C.P.D. and W.K. acknowledge funding from the DFG research group FOR 2634 “Planet Formation Witnesses and Probes: Transition Disks” under grant DU 414/22-1 and KL 650/29-1, 650/30-1. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 37/935-1 FUGG. All plots in this paper were made with the Python library matplotlib (Hunter 2007).

Appendix A Approaching the locally isothermal limit in HD 163296

To verify whether the locally isothermal results for HD 163296 concerning the primary and secondary gaps can be recovered by appropriately tweaking radiative effects, we tested several cases where we either restrained changes in internal energy by setting γ = 1.01, or amplifiedthe contribution of irradiation and cooling by lowering the opacity to κ = 0.045cm² g⁻¹. A detailed description of the model parameters used here is given in Table 1, Cols. 2 and 3.

Our results, plotted in Fig. A.1, show that locally isothermal conditions can be emulated toan extent by manipulating the contribution of the energy equation altogether through the adiabatic index, or controlling the cooling timescale through the opacity (Eq. (9) suggests that t_cool ∝ τ ∝ κ). It might be possible to completely recover the locally isothermal limit by further constraining γ → 1 or artificially lowering the opacity even further.

Fig. A.1 Azimuthally averaged surface density profiles after 1000 orbits for several comparison tests, attempting to recover the locally isothermal limit by constraining radiative effects. The orange and green curves correspond to our fiducial locally isothermal and radiative models, and the change being tested in each panel is plotted with a black dashed line. An adiabatic equation of state with γ = 1.01 or κ = 0.045 cm2/g (instead of 1.4 and 0.45 cm2/g, respectively) roughly reconstructs the gap structure of the locally isothermal model and could allow the formation of a shallower secondary gap. A model with lower viscosity (α = 10−5 instead of 10−4) also shows a secondary gap, but differences in the width of the primary gap and in the outer disk are also visible. The opacity model by Lin & Papaloizou (1985) yields no observable difference.

In addition, to test the planet’s secondary-gap opening capabilities under conditions that could potentially lead to the reemergence of the secondary gap, as well as to verify our opacity model of choice, we executed a few more pairs of locally isothermal and radiative simulations. We set α = 10⁻⁵ in the first pair to facilitate the gap opening process. We find that this does lead to the surface density contrast (and the pressure bump) reappearing in the inner disk, but at the same time, it also results in a wider gap and a different outer disk profile (see Fig. A.1). This value of α is extremely low and gap opening is known to be easier when the viscosity is low (e.g., Crida et al. 2006), therefore we argue that this adds nothing new to our results.

Fig. A.2 Azimuthally averaged surface density and midplane pressure profiles after 1000 orbits for two comparison tests where disks with different surface density profiles are simulated. Right panels: a shallower surface density profile (Σ0 ∝ r−1 instead of ∝ r−3∕2) is used.

In the second pair, we compared our constant-opacity models and the analytical opacity model by Lin & Papaloizou (1985), which dictates a relation κ ∝ T² for T < 170 K (i.e., within our simulation domain). We find that the choice of opacity model has little to no effect as far as the secondary-gap opening capabilities of the planet are concerned (Fig. A.1). This makes sense because irradiation heating and thermal cooling dominate the energy equation, and equating these two together factors the opacity out given our prescription in Eq. (2) while still yielding a similar value for the cooling timescale at the 20–40 au region.

Finally,in the third pair we prescribed a shallower initial surface density profile $\Sigma (r)\approx 10{(r/r_{\text{p}})}^{-1}$ in an attempt to reduce the cooling timescale in the inner disk while preserving the conditions in the planet’s vicinity. In doing so, we can also compare against the results from locally isothermal simulations by Zhang et al. (2018). The results are summarized in Fig. A.2 and paint a picture similar to that for the steeper, Σ(r) ∝ r^−3∕2 profile. However, because the cooling timescale is shorter for the shallower Σ profile (by a factor of 30 and 40% at the location of the secondary pressure bump and gap, respectively), the disk acts more “locally isothermally” and the secondary density bump is still visible after 1000 orbits but a pressure minimum does not form in the midplane. At the same time, pressure torques near the secondary gap edge are weaker for the shallower profile; this effect assists the gap-opening process. A detailed description of the model parameters used in these models is given in Table 1, Cols. 4–6.

Appendix B Modeling HD 163296 with three planets

As mentioned in Sect. 4, we conducted an additional simulation of HD 163296 with three planets at 48, 86, and 145 au to examine the structure of the system on a larger scale, more easily comparable to the DSHARP observation (see Fig. 1, left panel). To do so, we employed a fourth-order Runge–Kutta N-body integrator (Thun & Kley 2018) that allows the planets to interact gravitationally and accounts for the noninertial term that arises due to centering our system around r = 0 instead of the barycenter of the system.

In our model, all three planets have the same mass of 0.5 M_J. Because the thermal mass of the disk scales purely with the aspect ratio at a planet’s location, we expect a damping of nonlinear effects as we slowly transition to the linear regime upon moving farther out in the disk. More specifically, the ratio M_p ∕M_th is 1.3, 0.8, and0.5 for the three planets in ascending distance from the star. Because the cooling timescale at 86 and 145 au is 0.33 and 0.06% of the local orbital period, respectively, it is possible that we can still observe deviations between the locally isothermal and radiative models around the planets’ corotating regions because we now study the low-mass regime.

Fig. B.1 Azimuthally averaged surface density profile after 1000 orbits at 48 au (417 at 86 au, 190 at 145 au) for the three-planet model (see Table B.1). The colored dashed lines refer to the base model with a single planet at 48 au. Vertical black lines show the location of planets in the disk.

The resulting surface density profile is plotted in Fig. B.1. As expected, the two models show better overall overlap withincreasing distances from the star. Nevertheless, small deviations are still visible within the corotating regions of the outer planets. We rationalize this outcome both through the cooling timescale argument and by considering how the interaction between a planet and the disturbances caused by other planets (e.g., spirals) could affect their corotating regions. It should be noted, however, that the viscous timescale at 145 au is roughly three times longer than that at 48 au, so it is possible that differences between equations of state can emerge on much longer timescales, which might not be reasonable for such a young system.

Finally, by comparing the 10–60 au range of the three-planet simulation against that of our base model with a single planet (lower panel of Fig. B.1) we find that the multi-planet system shows stronger perturbations in surface density around the secondary gap at 28 au, but still fails to form a pressure maximum at 34 au in the radiative model. This helps support the relevance of our base model, which targeted a limited range of HD 163296.

References

All Tables

All Figures

Fig. 1 ALMA continuum emission observations of the two systems we used as testbeds for our numerical models, as shown in Huang et al. (2018). Solid and dotted white arcs mark the location of bright rings and dark gaps, respectively, along with their distance from the host star in au. It is likely that a planet is responsible for one or more of the gaps observed in either system. Left: HD 163296, a 12.6 Myr old likely isolated system that features at least two clearly visible gaps at 48 and 86 au. Our models focus on the 20–60 au range, with a planet carving a gap at 48 au. Right: AS 209, a 1 Myr old system in the Ophiuchus region that is rich in annular structures, boasting five rings between 20 and 130 au. This source has been modeled by Zhang et al. (2018), who found that a fit with a single planet at 99 au could reconstruct all five of these rings.

In the text

Fig. 2 Azimuthally averaged surface density, midplane pressure, and aspect ratio (used as a proxy for temperature, see Eq. (7)) of HD 163296 for a locally isothermal (in orange) and a radiative (in green) disk after 1000 orbits with an embedded planet with mass 0.5 MJ at 48 au. Vertical black lines mark the location of the planet. The dashed and dotted red lines in the bottom panel denote the aspect ratio profiles used in the simulations by Zhang et al. (2018), where h ∝ r1∕4 and hp = 0.05 and 0.07, respectively. The fact that our aspect ratio profile is bounded by these two curves allows us to validate our results by comparing them to that study (see Sect. 4).

In the text

	Fig. 3 Surface density perturbation with respect to the initial profiles for the two models of HD 163296. The locally isothermal run features a secondary gap at roughly 28 au and slightly more prominent inner spirals.
In the text

Fig. 4 Comparison of spiral arm properties between our two models for HD 163296. Top: primary and secondary spiral arm contrast with respect to the disk background. We find that the locally isothermal model shows consistently higher contrast. Bottom: pitch angle of spirals as a function of radius. Because the two models have identical aspect ratios, the overlap of the two curves is to be expected. The vertical black line marks the planet location. The dotted black line corresponds to the analytical formulae in Rafikov (2002), Eq. (44), and Muto et al. (2012), Eq. (1).

In the text

	Fig. 5 Temperature contrast of spirals with respect to the disk background for the radiative run of HD 163296. The effect of spiral heating on the global disk is negligible, and the temperature along spiral crests can increase by up to 15%.
In the text

Fig. 6 Surface density profile for our numerical model of AS 209 for a locally isothermal and a radiative equation of state after 5000 orbits. The differences in the inner disk as well as the gap region are clear between the two models. Vertical dotted lines mark location candidates for dust gaps based on the formula determined by Zhang et al. (2018), and the dashed line marks the planet location.

In the text

	Fig. 7 Similar to Fig. 3, but for AS 209. The sharper annular and angular structures in the locally isothermal model are clearly visible.
In the text

Fig. 8 Azimuthally averaged cooling timescale in units of the local orbital period of the disk after 1000 and 5000 planetary orbits for HD 163296 and AS 209, respectively. Differences between locally isothermal and radiative models become clear in the inner disk, where the cooling timescale corresponds to a non-negligible fraction of the orbital period (roughly 1–10% near the planet, and comparable to the orbital period farther in). In the outer disk, cooling is sufficiently fast so that the two models overlap as far as azimuthally averaged profiles are concerned.

In the text

	Fig. 9 Azimuthal cut at r = rp = 48 au. The surface density trailing behind the planet (left side, ϕ < 0) is slightly higher than that at the leading Lagrange point (ϕ > 0). The dashed black lines mark the planet’s hill radius and contain 20 cells. The dotted lines mark the location of the L4 (left) and L5 (right) Lagrange points.
In the text

Fig. A.1 Azimuthally averaged surface density profiles after 1000 orbits for several comparison tests, attempting to recover the locally isothermal limit by constraining radiative effects. The orange and green curves correspond to our fiducial locally isothermal and radiative models, and the change being tested in each panel is plotted with a black dashed line. An adiabatic equation of state with γ = 1.01 or κ = 0.045 cm2/g (instead of 1.4 and 0.45 cm2/g, respectively) roughly reconstructs the gap structure of the locally isothermal model and could allow the formation of a shallower secondary gap. A model with lower viscosity (α = 10−5 instead of 10−4) also shows a secondary gap, but differences in the width of the primary gap and in the outer disk are also visible. The opacity model by Lin & Papaloizou (1985) yields no observable difference.

In the text

	Fig. A.2 Azimuthally averaged surface density and midplane pressure profiles after 1000 orbits for two comparison tests where disks with different surface density profiles are simulated. Right panels: a shallower surface density profile (Σ0 ∝ r−1 instead of ∝ r−3∕2) is used.
In the text

	Fig. B.1 Azimuthally averaged surface density profile after 1000 orbits at 48 au (417 at 86 au, 190 at 145 au) for the three-planet model (see Table B.1). The colored dashed lines refer to the base model with a single planet at 48 au. Vertical black lines show the location of planets in the disk.
In the text