A global two-layer radiative transfer model for axisymmetric, shadowed protoplanetary disks

To properly treat shadows, we do not use an approximate irradiation surface defined in terms of the vertical optical depth as used in figure 1. We employ the exact irradiation surface defined as where the optical depth to the direct starlight is unity (D’Alessio et al., 2001). At position $(r,z)$ in a disk, the Planck-mean optical depth to the direct starlight is given by

where $\chi_{*}$ is the Planck-mean extinction opacity per gas mass for the starlight, $\rho$ is the gas mass density, and $s$ denotes the path length of the stellar ray from the star to the position considered. Thus, the height of the irradiation surface, $z_{s}(r)$, at arbitrary cylindrical radius $r$ is defined by the relation $\tau_{*}(r,z_{s}(r))=1$. The disk’s vertical structure is specified in section 3.4.

At large radial distances where the central star can be seen as a point source, stellar rays travel along lines of approximately constant $z/r$. For such regions, it is useful to map the irradiation surface in the $r$–$z/r$ plane. Figure 3 schematically shows the irradiation surface of the gapped disk shown in figure 1(b). As already mentioned in section 2, the irradiation surface outside the shadowed gap approximately matches the surface at which the vertical visual optical depth is $\sim H/r$. Inside the shadowed gap, the irradiation surface as seen in the $r$–$z/r$ plane is represented by a horizontal line. Below this line, the material inward of the gap blocks the direct starlight (see the lower panel of figure 3).

For given $z_{s}(r)$, one can calculate $\mu_{*}=\sin\alpha_{*}$ as (Kusaka et al., 1970; Ruden & Pollack, 1991)

In the right-hand side of equation (9), the first term accounts for the finite stellar size, while the sum of the second and third terms is related to the grazing angle of the ray from the central point source. At large radial distances where the first term is negligible, $\mu_{*}$ vanishes in shadowed regions with radially constant $z_{s}/r$ (see the lower panel of figure 3).

We define the thermal photosphere as the surface on which the optical depth to reprocessed starlight from the irradiation surface reaches unity (see figure 3). Strictly speaking, the location of the thermal photosphere depends on the incident angle of the reprocessed radiation. To avoid this complexity, we approximate the thermal photosphere by the surface where the vertical optical depth for downward infrared radiation exceeds $1/2$. This choice is based on the fact that the angle-averaged optical depth of an optically thin disk to isotropic radiation is twice the vertical optical depth (Nakamoto & Nakagawa, 1994). Writing the Planck-mean extinction opacity for infrared radiation as $\chi_{\rm P}$, the corresponding vertical optical depth is given by

Using this, the height of the thermal photosphere, $z_{\rm IR}(r)$, at radial position $r$ can be approximately defined by the relation $\tau_{\rm P}(r,z_{\rm IR}(r))=1/2$. The region well below the thermal photosphere is also optically thick to its own thermal emission.

For disks with $\tau_{\rm P}<1/2$ at the midplane, we set $z_{\rm IR}=0$ for convenience. However, the midplane of such an optically thin disk is not a real thermal photosphere, absorbing only a small fraction of the incoming infrared radiation. The low absorptivity and emissivity of optically thin disks are taken into account in section 3.5.

We split the irradiation surface into thin concentric ring elements of radius $r^{\prime}$, radial extent ${dr^{\prime}}$, and height $z_{s}(r^{\prime})$ (figure 4), Each ring emits infrared radiation with a luminosity proportional to the received stellar flux. To derive an analytic expression for the downward radiation flux from each ring element, we further divide the rings azimuthally into rectangular segments of azimuthal width $r^{\prime}d\phi^{\prime}$. The area of each segment is $dS=\sqrt{dr^{\prime 2}+dz_{s}^{2}}\,r^{\prime}d\phi^{\prime}=\sqrt{1+(dz_{s}/dr^{\prime})^{2}}\,r^{\prime}dr^{\prime}d\phi^{\prime}$. The luminosity of the downward radiation from each segment is

The radiation from each surface segment strikes the thermal photosphere. For simplicity, we here assume $dz_{\rm IR}/dr\ll 1$ and approximate the thermal photosphere as a plane locally parallel to the midplane. Because of this approximation, our model neglects the shadowing of the reprocessed radiation field on the thermal photosphere. For typical protoplanetary disks with $z_{\rm IR}\ll r$, the approximation made here can be justified unless $z_{\rm IR}$ varies on a radial lengthscale $\ll r$. Because we assume axisymmetry, the azimuthal coordinate of an arbitrary segment of the thermal photosphere can be taken to be zero without loss of generality. The angle between the reprocessed starlight and the thermal photosphere segment can then be written by

where

is the height of the irradiation surface segment relative to the photosphere segment and

is the distance between the two segments.

Assuming that the downward reprocessed starlight is isotropic, and approximating the irradiation surface segment as a point source, the magnitude of the downward reprocessed starlight flux $d{\bm{F}}_{\rm rep,\downarrow}$ (see figure 4) can be written as

where the factor $2\pi$ comes from the solid angle of the lower hemisphere. Because the reprocessed starlight strikes at angle $\alpha_{\rm rep}$, its flux incident on the thermal photosphere is $\sin\alpha_{\rm rep}|d{\bm{F}}_{\rm rep,\downarrow}|$. Integrating this over $\phi^{\prime}$, we obtain

where $E(k)=\int_{0}^{\pi/2}(1-k^{2}\sin^{2}x)^{1/2}dx$ is the complete elliptic integral of the second kind with $k=2\sqrt{rr^{\prime}/((r^{\prime}+r)^{2}+h^{2})}$. Finally, by integrating equation (16) over $r^{\prime}$, we obtain the net downward vertical flux $F_{\rm rep,\downarrow}(r)$ of the reprocessed starlight from the entire irradiation surface to the disk interior at radius $r$,

Here, the integration is performed over regions of $\alpha_{\rm rep}>0$ $(h>0)$³³3The quantity $h$ can become negative because it compares $z_{s}$ and $z_{\rm IR}$ at different radial locations. For instance, $z_{s}$ at small $r^{\prime}$ can be smaller than $z_{\rm IR}$ at large $r$.. Note that $f_{\downarrow}$, $\mu_{*}$, $F_{*}$, and $h$ generally depend on $r^{\prime}$ and hence should be inside the integral.

Equation (LABEL:eq:Fs) is a generalization of the reprocessed starlight flux in the local two-layer model, equation (1). To see this, we rewrite equation (LABEL:eq:Fs) as

where the function $W(r,r^{\prime})$ is

for $h>0$ and $W(r,r^{\prime})=0$ otherwise. Equation (18) can be viewed as the radial average of equation (1) weighted by $W(r,r^{\prime})$. As shown in appendix B, equation (18) recovers equation (1) in the limit of small $h$.

The weighting function $W$ represents the radially nonlocal nature of the reprocessed radiation. As explained below, this function is peaked around $r^{\prime}=r$ and has a width of $\sim h$. Because the elliptic integral $E$ is bounded in the narrow range of $1\leq E\leq\pi/2$ and because $\sqrt{1+(dz_{s}/dr^{\prime})^{2}}\sim 1$, the profile of $W$ is essentially determined by the factor $hr^{\prime}/[((r^{\prime}-r)^{2}+h^{2})\sqrt{(r^{\prime}+r)^{2}+h^{2}}]$. Since $h\ll r$, the profile around the peak is approximately Lorentzian $\propto 1/((r-r^{\prime})^{2}+h^{2})$ with the half width at half maximum of $h$. Figure 5 illustrates the functional form of $W$ for $z_{s}(r^{\prime})=0.15r^{\prime}$ and $z_{\rm IR}(r)=0.05r$, i.e., $h=0.15r^{\prime}-0.05r$. In this particular example, $W$ vanishes at $r^{\prime}/r<1/3$, at which $h$ becomes negative.

So far, we have not specified the disk’s vertical structure. In the present study, we avoid detailed modeling of the vertical structure by approximating the disk as vertically isothermal and hydrostatic. These two approximations yield the vertical distribution of $\rho$ of the simple analytic form

where $\Sigma$ is the gas surface density and $H=c_{s}/\Omega_{\rm K}$ is the gas scale height, with $c_{s}$ and $\Omega_{\rm K}$ being the isothermal sound speed and Keplerian frequency, respectively. The isothermal sound speed is related to the disk interior temperature as $c_{s}=\sqrt{k_{\rm B}T/m_{g}}$, where $k_{\rm B}$ is the Boltzmann constant and $m_{g}$ is the mean mass of the gas molecules.

We comment on the validity and limitations of the vertical isothermal and hydrostatic approximations. For passively irradiated disks in radiative equilibrium, the vertical isothermal approximation is valid well below the irradiation surface. However, this approximation breaks down across the irradiation surface, above which the temperature is higher because of direct stellar irradiation (Calvet et al., 1991; Chiang & Goldreich, 1997; D’Alessio et al., 1998). A more reasonable approximation would be to assign different temperatures to the interior ($z<z_{s}$) and warmer surface region ($z>z_{s}$) as described in appendix B of Watanabe & Lin (2008). However, such two-temperature modeling would result in a vertical density profile dependent on $z_{s}$, requiring iterative determination of $z_{s}$ and $\rho(z)$. We defer this complication to future work. Regarding the equilibrium disk structure, our vertical isothermal approximation should suffice to determine $z_{s}$ because the warmer but optically thin gas well above the irradiation surface should have little effect on the location of the irradiation surface (see Watanabe & Lin, 2008). Nevertheless, the vertical variation of the temperature, and more importantly of the thermal relaxation time, could have a more critical impact on the time evolution of the disk temperature. We discuss this point in more detail in section 3.5.

The vertical hydrostatic approximation is valid as long as the disk’s thermal relaxation time $t_{\rm th}$ (see equation (27) for definition) is much longer than the orbital timescale $t_{\rm K}=2\pi/\Omega_{\rm K}$. In protoplanetary disks, the condition $t_{\rm th}\gg t_{\rm K}$ is typically fulfilled at $r\ll 100~{}\rm au$ (Dullemond, 2000; Watanabe & Lin, 2008; Wu & Lithwick, 2021, see also section 5.2) but can break down farther out. As discussed by Wu & Lithwick (2021), the disk’s hydrodynamic response may influence its thermal and dynamical evolution. In this study, we concentrate on the relatively inner disk region where the vertical hydrostatic assumption is applicable.

The opacities used in our model depend on the dust-to-gas mass ratio and size distribution of the opacity-dominating dust grains. One can account for dust settling by using vertically varying opacities.

Because passively irradiated disks can be thermally unstable, we treat the disk interior temperature as intrinsically time-dependent. For simplicity, we neglect radial heat advection by accreting gas (see, e.g., Cannizzo, 1993), radial diffusion of the disk’s own thermal radiation (e.g., Latter & Balbus, 2012; Owen & Armitage, 2014), and accretion heating. We plan to include these effects in future work.

Neglecting the above mentioned effects, the time evolution of the interior temperature $T$ obeys the following energy equation,

where $\gamma$ is the adiabatic index (taken to be $1.4$ throughout this paper), $T_{\rm ex}$ is the radiation temperature of the parent molecular cloud (Ueda et al., 2021), and ${\cal C}$ is a dimensionless factor correcting for the effects of the disk’s infrared optical thickness and albedo.

Equation (21) is the vertically integrated equation of total energy conservation (Watanabe et al., 1990; Watanabe & Lin, 2008). Its left-hand side stands for the rate of change in the total energy per unit disk area. The temperature on the left-hand side originally stands for the density-weighted vertical average of the temperature (Watanabe et al., 1990). Following Watanabe & Lin (2008), we have applied the vertical isothermal approximation and represented the vertically averaged temperature with the single interior temperature. In reality, the rate of change in temperature in optically thick disks depends on the depth from the disk surface, with shallower and optically thinner regions having shorter cooling timescales. Recently, Pavlyuchenkov et al. (2022) have pointed out that the vertical variation of the thermal relaxation timescales could weaken the thermal wave instability around the midplane of optically thick disks. This effect is not included in our current modeling.

The right-hand side of Equation (21) accounts for radiative heating and cooling on both sides ($z<0$ and $z>0$) of the disk. We take the correction factor ${\cal C}$ to be

where

is the ratio between the Planck-mean absorption opacity $\kappa_{\rm P}$ (see equation (39) for definition) and Rosseland extinction opacity $\chi_{\rm R}$,

is the effective absorption optical depth to the midplane accounting for multiple scattering (Rybicki & Lightman, 1979), and

is the Rosseland-mean vertical optical depth to the midplane. All the mean opacities appearing here are for the disk thermal radiation. In equation (22), the factor involving $\tau_{\rm eff,mid}$ corrects for the infrared emissivity and absorptivity of both optically thick and thin disks, derived by approximating a disk with a uniform slab (see equation (53) in appendix D). The factor $(1+3\tau_{\rm R,mid}/4)^{-1}$ corrects for the vertical radiative diffusion flux in optically thick ($\tau_{\rm R,mid}\gg 1$) disks being inversely proportional to $\tau_{\rm R,mid}$ (Wu & Lithwick, 2021).

We consider a radial computational domain spanning $r=r_{\rm min}$ to $r=r_{\rm max}$ and discretize it into logarithmically spaced cells. Each radial cell has the inner and outer boundaries at $r_{j-1/2}=r_{\rm min}(r_{\rm max}/r_{\rm min})^{(j-1)/J}$ and $r_{j+1/2}=r_{\rm min}(r_{\rm max}/r_{\rm min})^{j/J}$, respectively, and the logarithmic center at $r_{j}=\sqrt{r_{j-1/2}r_{j+1/2}}=r_{\rm min}(r_{\rm max}/r_{\rm min})^{(j-1/2)/J}$, where $j=1,2,\dots J$ labels the cells and $J$ is the number of the cells. One should take the cell width should to be sufficiently smaller than the width $\sim h$ of the weighting function $W$.

The irradiation surface is found with a ray-tracing approach. We use rays emanating from the coordinate origin $(r,z)=(0,0)$ at angle $\theta$ with respect to the $r$-axis. All simulations presented in this paper adopt 180 linearly spaced $\theta$ grids spanning $\theta=0$ to $\theta=\theta_{\rm max}=\pi/12$, with a resolution of $d\theta=\theta_{\rm max}/180=1.45\times 10^{-3}$. We have also tried simulations with $d\theta=\theta_{\rm max}/360$ and confirmed that the higher angular resolution gives no appreciable change in the simulation results.

When searching for the irradiation surface, one must assume the starlight optical depth to the inner computational boundary. The starlight optical depth $\tau_{*}$ (equation (8)) to arbitrary radial position $r$ along a ray with angle $\theta$ can be written as

where $\tau_{*,\rm in}$ denotes the optical depth to the inner computational boundary. The second term in the right hand side of equation (26) uses $z=r\tan\theta$ and $ds=dr/\cos\theta$. The problem here is that $\tau_{*,\rm in}$ is intrinsically unknown because it is determined by the disk structure outside the computational domain. However, as demonstrated in appendix C, an unreasonable choice of $\tau_{*,\rm in}$ can cause unwanted artifacts in the resulting temperature distribution near the inner boundary. Our prescription for $\tau_{*,\rm in}$ is presented in appendix C.

The energy equation (21) is solved with a first-order forward differencing scheme. The time step must be smaller than the thermal relaxation timescale for the vertically averaged temperature (Watanabe & Lin, 2008),

Note that $t_{\rm th}$ depends on the correction factor ${\cal C}$ introduced in section 3.5.

Because the radial computational domain is limited to $r_{\rm min}<r^{\prime}<r_{\rm max}$, the radial integration in equation (LABEL:eq:Fs) cannot be extended to $r^{\prime}<r_{\rm min}$ and $r^{\prime}>r_{\rm max}$. However, neglecting the reprocessed starlight from outside the computational domain would result in an underestimate of the temperature near the computational boundaries. Near the inner computational boundary, an underestimated temperature amplifies the numerical wiggle in the radial disk structure (see appendix C). To mitigate the underestimation of the temperature, we account for the flux outside the computational domain in an approximate way. Specifically, we write $F_{\rm rep,\downarrow}$ as

where $F_{\rm rep,\downarrow,\rm domain}$ represents the flux from within the computational domain, whereas $F_{\rm rep,\downarrow,in}$ and $F_{\rm rep,\downarrow,out}$ represent the additional fluxes from $r^{\prime}<r_{\rm min}$ and $r^{\prime}>r_{\rm max}$, respectively. Our prescription for these additional fluxes is described in appendix C.

To suppress numerical instabilities arising from grid-scale fluctuations of the reprocessed starlight flux, we replace $F_{\rm rep,\downarrow}(r_{j})$ at $j=2,3,\dots,J-1$ with the geometric means of the raw fluxes at $j-1$ and $j+1$, $\sqrt{F_{\rm rep,\downarrow}(r_{j-1})F_{\rm rep,\downarrow}(r_{j+1})}$. At $j=1$ and $J$, we use $\sqrt{F_{\rm rep,\downarrow}(r_{1})F_{\rm rep,\downarrow}(r_{2})}$ and $\sqrt{F_{\rm rep,\downarrow}(r_{J-1})F_{\rm rep,\downarrow}(r_{J})}$ instead. We expect that this grid-scale smoothing would become unnecessary once we include the radial diffusion of the disk’s thermal emission in equation (21) in future work.

Following Jang-Condell & Turner (2012), we consider axisymmetric gapped disks with the radial gas surface density profile given by

where $19({r}/{10~{}\rm au})^{-0.9}~{}\rm g~{}cm^{-2}$ is the surface density profile with no planet, $a=10~{}\rm au$ is the planet’s orbital radius, and $d$ and $w$ represent the depth and width of the planet-carving gap, respectively. They considered one case with no planet and two cases with a planet of mass $M_{p}=70$ or $200M_{\oplus}$ orbiting at 10 au from the central star. The parameters $(d,w)$ were taken to be $(0.56,0.11a)$ and $(0.84,0.17a)$ for $M_{p}=$ 70 and 200$M_{\oplus}$, respectively. The top row of figure 6 shows the gas surface density profiles for three disk models. The central star was assumed to have mass $M_{*}=M_{\odot}$, radius $R_{*}=2.6R_{\odot}$, surface temperature $T_{*}=4280~{}\rm K$.

Jang-Condell & Turner (2012) provided the temperature structure of the three disk models ($M_{p}=0$, 70, and $200M_{\oplus}$) using two radiative transfer models. One is a model conceptually similar to ours, relying on the locally one-dimensional analytic solution for a plane-parallel disk (Jang-Condell & Sasselov, 2003, 2004; Jang-Condell, 2008, 2009). The other is the Monte Carlo model developed by Turner et al. (2012), which directly follows the emission, absorption, and scattering of photon packets using the frequency-dependent opacities provided by Jang-Condell (2009). The two approaches returned similar results as shown in figure 4 of Jang-Condell & Turner (2012), so we only use their Monte Carlo results in the following comparison. The dashed lines in the second row of figure 6 show the temperature profiles from the Monte Carlo calculations by Jang-Condell & Turner (2012). We note that their Monte Carlo calculations only cover 3–20.8 au.

We adopt the mean opacities for stellar and disk thermal radiation computed from the frequency-dependent opacities of Jang-Condell (2009). For disk thermal radiation, we evaluate the mean opacities at a fixed temperature of $T=50~{}\rm K$. The adopted mean opcities are $\chi_{*}=11.3~{}\rm cm^{2}~{}g^{-1}$, $\kappa_{*}=1.36~{}\rm cm^{2}~{}g^{-1}$, $\chi_{\rm P}=1.92~{}\rm cm^{2}~{}g^{-1}$, $\chi_{\rm R}=1.98~{}\rm cm^{2}~{}g^{-1}$, and , $\kappa_{\rm P}=0.975~{}\rm cm^{2}~{}g^{-1}$, which yield $\epsilon_{*}=0.12$, $q=5.7$, $f_{\downarrow}={0.24}$ (see figure 2), and $\epsilon=0.49$. The disk model considered here is optically thick to its own thermal emission, with $\tau_{\rm R}\sim 10$ at $r\sim 10~{}\rm au$. The thermal relaxation time at $r\sim 10~{}\rm au$ is $t_{\rm th}\sim k_{\rm B}\Sigma/(m_{g}\sigma_{\rm SB}T^{3}\tau_{\rm R})\sim 10^{2}~{}\rm yr$.

Unlike the radiative equilibrium Monte Carlo calculations of Jang-Condell & Turner (2012), our global two-layer calculations are fully time-dependent. The time step in our calculations is fixed to be 1 yr, which is shorter than the thermal relaxation time in the computational domain. The initial temperature profile is set to be $T=120(r/1~{}\rm au)^{-3/7}~{}\rm K$. This choice is arbitrary, but is not far from the steady state temperature profile for the model with no planet.

Our two-layer calculations use a computational domain covering $r_{\rm min}=1$ au to $r_{\rm max}=$ 100 au, divided into 200 logarithmically spaced radial cells. Our computational domain is wider than in the Monte Carlo calculations by Jang-Condell & Turner (2012), and therefore our results are likely less affected by the computational boundaries and by any outer edge to the disk. We neglect radiation from the parent molecular cloud by setting $T_{\rm ex}=0$ in equation (21). The additional reprocessed starlight fluxes from outside the computational domain, $F_{\rm rep,\downarrow,in}$ and $F_{\rm rep,\downarrow,out}$, are included.

Figure 7 shows the time evolution of the temperature $T$ for three Jang-Condell & Turner (2012) disk models obtained from our global two-layer calculations. In the $M_{p}=0$ and $70M_{\oplus}$ models, the temperature profile relaxes into a steady state. During relaxation, the temperature exhibits a damped oscillation with a period of $\approx 700~{}\rm yr$. This oscillation period is comparable to the thermal relaxation time of the disk, which is $\sim 200$–400 yr for the model with no planet. This implies that the observed oscillation is not an artifact but reflects the system’s thermal relaxation. In the $M_{p}=200M_{\oplus}$ disk model, an oscillation of a similar period persists at a constant amplitude. The oscillation propagates toward the star, indicating that it is likely due to the thermal wave instability (see section 5 for more details about the long-term behavior of this instability). In the following, we describe the results for each disk model in more detail.

We adopt the disk model used by Ueda et al. (2021). The gas surface density profile is given by $\Sigma=1700(r/1~{}{\rm au})^{-3/2}\exp(-r/100~{}\rm au)~{}g~{}cm^{-2}$. The disk opacities are assumed to scale with the dust-to-gas mass ratio $f_{\rm d2g}$ of opacity-dominating dust grains. For the fiducial value of $f_{\rm d2g}=0.01$, the mean opacities are $\chi_{*}=\kappa_{*}=8~{}\rm cm^{2}~{}g^{-1}$ and $\chi_{\rm R}=\chi_{\rm P}=\kappa_{\rm P}=4~{}\rm cm^{2}~{}g^{-1}$. for dust grains is a free parameter of the model. In this study, we only consider this fiducial model. Because $\epsilon_{*}=1$, we have $f_{\downarrow}=0.5$. External thermal radiation of $T_{\rm ex}=10~{}\rm K$ is included, so that the disk temperature never falls below 10 K.

The radial computational domain ranges between 0.03 to 300 au and is divided into 480 logarithmically spaced cells, resulting in a radial resolution of $dr/r=0.019$. The adopted radial resolution is two times better than in the simulation by Ueda et al. (2021). In section E, we also present simulation runs with lower resolution to study the convergence of our simulation results.

Our simulations differ from those by Ueda et al. (2021) in the treatment of radial radiative transfer and thermal relaxation. To discriminate between the effects of the two different treatments, we also conduct a simulation using the same local radiative transfer model as in the simulations by Ueda et al. (2021), but including the correction for thermal relaxation in optically thick regions, the factor $(1+3\tau_{\rm R,mid}/4)^{-1}$ in equation (22). Their radiative transfer model, originally developed by Watanabe & Lin (2008), approximates the downward reprocessed starlight flux as

where the angled brackets denote a radial average over a radial zone near $r^{\prime}=r$. This averaging is introduced to mimic the radial propagation of reprocessed starlight over distance $\sim z_{s}(r)$. Following Ueda et al. (2021), we use the Gaussian weight function $\exp[-(r^{\prime}-r)^{2}/z_{s}(r)^{2}]$ for the radial averaging in equation (30).

The simulations are carried out over 1 Myr. The simulated time covers $\sim 100t_{\rm th}$ at 10 au and $\sim 10t_{\rm th}$ at 1 au (see section 5.2). In comparison, the simulation reported by Wu & Lithwick (2021) only covers several thermal times at 10 au and less than one thermal time at 1 au. Therefore, our simulation for the first time fully captures thermal waves traveling from the 10 au to 1 au regions, with a correct treatment for thermal relaxation in the optically thick regime. The timestep $\Delta t$ is taken to be 1 yr, which is 10 times shorter than the minimum thermal relaxation time in the disk ($\sim 10~{}\rm yr$; see figure 11 in section 5.2). Using a shorter timestep of $\Delta t=0.25$ yr makes no appreciable change in the simulation results.

All our simulations assume vertical hydrostatic equilibrium (see section 3.4). As shown in section 5.2, the requirement $t_{\rm th}\gg t_{\rm K}$ for the hydrostatic equilibrium is fulfilled over the entire thermally unstable region.

Figure 10 presents snapshots of the temperature distribution at selected times. We find that the temperature structure at $r\approx 0.3$–30 au is unstable and exhibits oscillations (thermal waves) that propagate inward. The thermal waves consist of sharp temperature peaks with width $\sim z_{s}(r_{\rm peak})$, where $r_{\rm peak}$ is the radial position of each peak. Each peak casts a shadow with a radial width $\sim r_{\rm peak}$ (see the middle and bottom panels of figure 10), and each shadow causes a dip in the temperature profile. The irradiation surface height in the unstable region exceeds $4H$, which is higher than in the Jang-Condell & Turner (2012) disk models. This provides further support for the prediction by (Wu & Lithwick, 2021) that disks with larger $z_{s}/H$ are more prone to the thermal wave instability.

At $r\gtrsim 30~{}\rm au$ and $r\lesssim 0.3~{}\rm au$, the thermal wave instability is suppressed for different reasons. In the outer region of $r\gtrsim 30~{}\rm au$, the external radiation flux $\sigma_{\rm SB}T_{\rm ex}^{4}$ is comparable to or even dominates over the reprocessed starlight flux $F_{\rm rep,\downarrow}$, directly suppressing the thermal wave instability (Ueda et al., 2021). In the inner region of $r\lesssim 0.1~{}\rm au$, the finite size of the central star determines the starlight grazing angle (i.e., $\alpha_{*}\approx 4R_{*}/(3\pi r)$; see equation (9)), which stabilizes the thermal wave instability as it is triggered by the variation of the stellar grazing angle with temperature fluctuations. In addition, the thermal timescale in this inner region is too long for thermal waves to develop within 1 Myr. Indeed, figure 11, shows that the thermal relaxation time $t_{\rm th}$ (equation (27)) at $r\lesssim 0.1~{}\rm au$ exceeds 1 Myr. Figure 11 also shows that the condition $t_{\rm th}\gg t_{\rm K}$ for vertical hydrostatic equilibrium is fulfilled in the entire region where the thermal wave instability is observed.

Compared to the thermal waves observed in the simulations by Ueda et al. (2021, see their figures 2 and 3), our thermal waves are less coherent and propagate inward on different timescales at different radial positions. This can be more clearly seen in figure 12, which shows a spacetime plot of the temperature distribution. To highlight the wave components, we here normalize the temperature profile by the time-average $T_{\rm avr}(r)$ over $t=0.1$–1 Myr. One can see that temperature peaks at smaller $r$ migrate inward on longer timescales. The observed migration timescale at each $r$ is crudely consistent with the local thermal timescale $t_{\rm th}$ shown in figure 11. This suggests that the radial variation of the migration speed is a manifestation of the radial variation of $t_{\rm th}$. In contrast, in the simulations by Ueda et al. (2021), $t_{\rm th}$ was nearly independent of $r$ because they did not include the correction for the thermal relaxation rate. This explains why the thermal waves observed by Ueda et al. (2021) propagate inward on a radially uniform timescale. For the Ueda et al. (2021) disk model with $f_{\rm d2g}=0.01$, the thermal relaxation correction must be included because the disk is optically thick ($\tau_{\rm R,\rm mid}>1$) at $r\lesssim 100~{}\rm au$ (see figure 11).

Another interesting feature that was not visible in the simulations by Ueda et al. (2021) is the collision of temperature peaks. As an example, the top left panel of figure 10 shows how two temperature peaks labelled by 1 and 2 collide. The outer peak 2 migrates faster than the inner peak 1, and they merge together at $t\approx 0.40~{}\rm Myr$. We speculate that the radially varying migration timescales of the temperature peaks induce their collisions.

Because the observed temperature peaks are narrow, it is important to examine whether the finite numerical resolution affects our simulation results. In appendix E, we show that a radial resolution of $dr/r\lesssim 0.04$ is enough to resolve the thermal waves.

Figure 13 shows the spacetime temperature plot from the simulation using the radiative transfer model of Watanabe & Lin (2008). The results are qualitatively similar to those from the global two-layer simulation presented above in that both feature collisions of temperature peaks. However, the Watanabe & Lin (2008) model produces wider temperature peaks than our global two-layer model. This is likely due to the ad hoc radial averaging of the reprocessed starlight flux adopted in the model. Therefore, we conclude that the radiative transfer model of Watanabe & Lin (2008) captures the qualitative nature of the thermal wave instability but should not be used for a quantitative study of the instability.

As already noted by previous studies (Watanabe & Lin, 2008; Ueda et al., 2021; Wu & Lithwick, 2021), the thermal wave instability may have two important implications for dust evolution and planetesimal formation in dust-rich disks. First, the thermal wave instability causes temporal variations in snow line locations. We demonstrate this in figure 14, where we show how the positions where $T=$ 160, 70, and 20 K move with time in the simulation shown in figures 10 and 12. The selected temperatures correspond to the sublimation temperatures of H₂O, CO₂, and CO ices, respectively (Okuzumi et al., 2016). Overall, the plot indicates that the instability causes order-of-unity time variations in the positions of the snow lines. For instance, the H₂O snow line would migrate between $\approx 0.4$–1.5 au on the timescale of $\sim 0.1~{}\rm Myr$. The previous simulation by Ueda et al. (2021) predicted a migration timescale of $\sim 10~{}\rm yr$ for the H₂O snow line. Our results indicate that the simulations by Ueda et al. (2021) underestimated the migration timescale of the H₂O snow line because their simulations did not include the correction for thermal relaxation in the optically thick limit. Our new simulation suggests that the H₂O snow line can in fact move on a timescale comparable to the planet formation timescale, $\sim 0.1$–1 Myr, if the disk’s optical thickness is large enough. We expect that the snow line migration should have important effects on planetesimal formation around the snow line (Ros & Johansen, 2013; Schoonenberg & Ormel, 2017; Dra̧żkowska & Alibert, 2017; Ida et al., 2021; Hyodo et al., 2021), in particular on the composition of forming planetesimals. In future work, we will explore the potential effects by including dust evolution in our global two-layer radiative transfer calculations.

Second, the thermal wave instability might produce pressure maxima that can trap dust particles. Dust particles in a gas disk are known to drift in the direction of increasing gas pressure (Whipple, 1972; Adachi et al., 1976; Weidenschilling, 1977). Specifically, the radial drift velocity of a dust particle in a gas disk is given by

where $t_{\rm stop}$ is the stopping time of the particle, $v_{\rm K}=r\Omega_{\rm K}$ is the local Keplerian velocity, and $P=\rho k_{\rm B}T/m_{g}$ is the gas pressure. The radial drift velocity is proportional to the radial pressure gradient ${\partial\ln P}/{\partial\ln r}={\partial\ln\rho}/{\partial\ln r}+{\partial\ln T}/{\partial\ln r}$, with negative and positive radial gradients leading to inward and outward drift, respectively. A pressure maximum acts as a trap for radially drifting particles because their drift velocity converges there. Watanabe & Lin (2008) already suggested that the thermal wave instability can potentially produce local pressure maxima. Our new simulations relying on a more rigorous radiative transfer model confirm this possibility. Figure 15 shows some snapshots of the radial distribution of $\partial\ln P/\partial\ln r$ at the midplane from our global two-layer simulation presented in Figure 10. We find that some of the temperature peaks observed in figure 10 (e.g., peak 2^′ shown in the right panel of figure 10) indeed reverse the sign of the pressure gradient locally. Less pronounced temperature peaks that do not reverse the pressure gradient may also slow down particle inward drift and thereby promote planetesimal formation. However, a question remains as to whether the steep pressure variations observed here are hydrodynamically stable, because they may violate the Rayleigh criterion for stable disk rotation (e.g., Yang & Menou, 2010). More fundamentally, the actual thermal waves around the midplane of optically thick disks could be less pronounced than predicted from our calculations owing to the vertically varying thermal relaxation time (Pavlyuchenkov et al., 2022). Addressing these issues is beyond the scope of this paper but should be done in future work.