# Spirals, rings, and vortices shaped by shadows in protoplanetary discs: from radiative hydrodynamical simulations to observable signatures

Alexandros Ziampras, Cornelis P. Dullemond, Tilman Birnstiel, Myriam Benisty, Richard P. Nelson

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

## To cite this version:

Alexandros Ziampras, Cornelis P. Dullemond, Tilman Birnstiel, Myriam Benisty, Richard P. Nelson. Spirals, rings, and vortices shaped by shadows in protoplanetary discs: from radiative hydrodynamical simulations to observable signatures. Monthly Notices of the Royal Astronomical Society, 2025, 540, pp.1185-1201. 10.1093/mnras/stf785. insu-05099854

HAL Id: insu-05099854

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

Submitted on 5 Jun 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spirals, rings, and vortices shaped by shadows in protoplanetary discs: from radiative hydrodynamical simulations to observable signatures

Alexandros Ziampras<sup>1,2,3★</sup> Cornelis P. Dullemond<sup>2</sup>, Tilman Birnstiel<sup>2,5</sup> Myriam Benisty<sup>3,6</sup>
and Richard P. Nelson<sup>1</sup>

Accepted 2025 April 24. Received 2025 April 21; in original form 2024 October 16

## ABSTRACT

Numerous protoplanetary discs exhibit shadows in scattered light observations. These shadows are typically cast by misaligned inner discs and are associated with observable structures in the outer disc, such as bright arcs and spirals. Investigating the dynamics of the shadowed outer disc is therefore essential in understanding the formation and evolution of these structures. We carry out two-dimensional radiation hydrodynamics simulations that include radiative diffusion and dust-gas dynamics to study the formation of substructure in shadowed discs. We find that spiral arms are launched at shadow edges, permeating the entire disc. The local dissipation of these spirals results in an angular momentum flux, opening multiple gaps, and leading to a series of concentric, regularly spaced rings. We find that ring formation is favoured in weakly turbulent discs where dust growth is taking place. These conditions are met for typical Class II discs, in which bright rings should form well within a fraction of their lifetime (~0.1–0.2 Myr). For hotter discs, gap opening is more efficient, such that the gap edges quickly collapse into vortices that can appear as bright arcs in continuum emission before decaying into rings or merging into massive, long-lived structures. Synthetic observations show that these structures should be observable in scattered light and millimetre continuum emission, providing a new way to probe the presence of substructure in protoplanetary discs. Our results suggest that the formation of rings and gaps is a common process in shadowed discs and can explain the rich radial substructure observed in several protoplanetary discs.

**Key words:** accretion, accretion discs – hydrodynamics – radiation: dynamics – methods: numerical.

## 1 INTRODUCTION

The last decade of observations has revealed that shadows are a common feature in protoplanetary discs around young stars. These shadows, typically theorized to be cast by misaligned inner discs or warps (e.g. Marino, Perez & Casassus 2015; Muro-Arena et al. 2020; Bohn et al. 2022), can span a wide azimuthal extent of the outer disc in scattered light observations when the inner disc is only tilted by a few degrees, and are associated with a variety of observed structures such as arcs and spirals. Such examples include HD 139614 (Muro-Arena et al. 2020), HD 143006 (Benisty et al. 2018), and Wray 15–788 (Bohn et al. 2019). In the case of highly inclined inner discs, the shadow can be more localized and appear as a dark lane, blocking incident stellar radiation in a narrow azimuthal range. Examples of narrow shadows include HD 100453 (Benisty et al. 2017), HD 135344B (Stolker et al. 2017), DoAr 44 (Avenhaus

et al. 2018; Casassus et al. 2018), and RX J1604.3–2130A (Pinilla et al. 2018).

The presence of one or more shadows or warps in a protoplanetary disc has been shown to significantly influence both its thermal and chemical structure (Young et al. 2021) as well as its dynamical evolution (e.g. Juhász & Facchini 2017; Kimmig & Dullemond 2024). Gas in the shadowed region can cool quite efficiently, especially so in the optically thin outer disc (e.g. Casassus et al. 2019; Nealon, Price & Pinte 2020). The resulting temperature contrast between the shadowed and illuminated regions translates to a radial force imbalance, as pressure support is weaker inside the shadow, which manifests as spirals in the disc (Montesinos et al. 2016; Montesinos & Cuello 2018; Zhang & Zhu 2024). These spirals provide an additional interpretation to spiral structures besides embedded planets (e.g. Zhu et al. 2015) and turbulence due to the gravitational instability (e.g. Gammie 2001; Béthune, Latter & Kley 2021).

At the same time, such vigorous spiral activity would provide an efficient mechanism to redistribute angular momentum within the disc. The non-linear damping of spiral arms can drive a local angular

<footnote n="1"><sup>1</sup>Astronomy Unit, School of Physics and Astronomy, Queen Mary University of London, London E1 4NS, UK</footnote>

<footnote n="2"><sup>2</sup>Ludwig-Maximilians-Universität München, Universitäts-Sternwarte, Scheinerstr 1, D-81679 München, Germany</footnote>

<footnote n="3"><sup>3</sup>Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="4"><sup>4</sup>Institute for Theoretical Astrophysics, Center for Astronomy (ZAH), Heidelberg University, Albert Ueberle Str 2, D-69120 Heidelberg, Germany</footnote>

<footnote n="5"><sup>5</sup>Exzellenzcluster Origins, Boltzmannstr 2, D-85748 Garching, Germany</footnote>

<footnote n="6"><sup>6</sup>Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, F-06304 Nice Cedex 4, France</footnote>

<footnote>*E-mail: a.ziampras@lmu.de</footnote>
momentum flux (AMF) into the disc due to them steepening into shocks (e.g. Rafikov 2002) at the interface between the shadowed and illuminated regions, or due to radiative cooling (Miranda & Rafikov 2020a; Zhang & Zhu 2020; Ziampras, Nelson & Rafikov 2023). This process essentially acts to remove material from the location where damping occurred, ultimately leading to the formation of a gap if the AMF is both sustained long enough at a given radial location and strong enough to overcome the gap-refilling effects of turbulence and local pressure gradients (Crida, Morbidelli & Masset 2006). Due to the long time-scales involved in the gap-opening process, however, the possibility of gap opening due to the spiral arms excited in shadowed discs has not been explored in detail (but see Su & Bai 2024, whose work was published during the preparation of this paper).

In this work, we aim to investigate the formation of radial structure in the form of rings and gaps due to the spiral wakes excited in a protoplanetary disc shadowed by a misaligned inner disc. We approach this problem with radiation hydrodynamics simulations that include a self-consistent treatment of cooling and radiative diffusion, coupled to a dust component that accounts for the dynamics of both submicron and millimetre grains. In this way, we can study both the thermal and dynamical evolution of the gas structure in the presence of a shadow, as well as the observability of any resulting substructure in scattered light or millimetre continuum observations.

In Section 2, we introduce our physical and numerical framework and describe our model of a shadowed protoplanetary disc. We analyse our results from a suite of hydrodynamical simulations in Sections 3 and 4, and present a set of synthetic observations in Section 5. We discuss our findings in Section 6, and finally conclude in Section 7.

## 2 PHYSICS AND NUMERICS

In this section, we introduce our physical framework and describe the approach we take to model a shadowed protoplanetary disc. We also outline our numerical set-up and list the parameters we use in our simulations.

### 2.1 Radiation hydrodynamics

We solve the Navier-Stokes equations (e.g. Tassoul 1978) in cylindrical coordinates {$R, \phi$} for a vertically integrated distribution of perfect gas with adiabatic index $\gamma = 7/5$ and mean molecular weight $\mu = 2.35$ orbiting around a star with mass $M_*$ and luminosity $L_*$. We also include a distribution of dust grains with bulk density $\bar{\rho}_d$ and size $a_d$, modelled as a pressureless fluid. Defining the gas and dust surface density $\Sigma_g$ and $\Sigma_d$, the gas and dust velocities $u_g$ and $u_d$, and the vertically integrated gas pressure $P$, the equations read

$$ \frac{d\Sigma_g}{dt} = -\Sigma_g \nabla \cdot u_g, \qquad (1a) $$

$$ \frac{d\Sigma_d}{dt} = -\Sigma_d \nabla \cdot u_d - \nabla \cdot j, \quad j = -\Sigma \frac{\nu}{Sc} \frac{1 + 4\text{St}^2}{(1 + \text{St}^2)^2} \nabla \left( \frac{\Sigma_d}{\Sigma} \right), \qquad (1b) $$

$$ \Sigma_g \frac{du_g}{dt} = -\nabla P - \Sigma_g \nabla \Phi_* + \nabla \cdot \bar{\sigma} - \Sigma_d \frac{u_g - u_d}{St} \Omega_K, \qquad (1c) $$

$$ \Sigma_d \frac{du_d}{dt} = -\Sigma_d \nabla \Phi_* - \Sigma_d \frac{u_d - u_g}{St} \Omega_K, \qquad (1d) $$

$$ \frac{de}{dt} = -\gamma e \nabla \cdot u_g + Q_{\text{visc}} + Q_{\text{irr}} + Q_{\text{cool}} + Q_{\text{rad}}. \qquad (1e) $$

Here, $\Sigma = \Sigma_g + \Sigma_d$ is the total surface density, $e = P/(\gamma - 1)$ is the internal energy density, $\Phi_* = -GM_*/R$ is the gravitational potential of the central star of mass $M_*$ at distance $R$, $\bar{\sigma}$ is the viscous stress tensor, and $\Omega_K = \sqrt{GM_*/R^3}$ is the Keplerian angular velocity, with $G$ being the gravitational constant. We can then define the gas temperature $T = \mu c_s^2/R$ with $c_s = \sqrt{P/\Sigma_g}$ the isothermal sound speed, which relates to the pressure scale height $H = c_s/\Omega_K$. The disc aspect ratio is then $h = H/R$.

In equation (1b), $j$ represents the dust diffusion flux (Morfill & Voelk 1984; Youdin & Lithwick 2007) for a kinematic viscosity $\nu$ and assuming that the Schmidt number is $Sc = 1$. Finally, $St$ denotes the Stokes number or dimensionless stopping time of the dust grains, which are subject to Epstein drag (e.g. Armitage 2009) with

$$ \text{St} = \frac{\pi a_d \bar{\rho}_d}{2 \Sigma_g}. \qquad (2) $$

The terms $Q_{\text{visc}}$, $Q_{\text{irr}}$, $Q_{\text{cool}}$, and $Q_{\text{rad}}$ represent viscous heating, irradiation heating, surface cooling, and in-plane radiative diffusion, respectively:

$$ Q_{\text{visc}} = \frac{1}{2\nu\Sigma_g} \text{Tr}(\bar{\sigma}^2) \approx \frac{9}{4}\nu\Sigma_g\Omega_K^2, \quad \nu = \alpha\sqrt{\gamma c_s H}, \qquad (3a) $$

$$ Q_{\text{irr}} = 2 \frac{L_*}{4\pi R^2} (1 - \epsilon) \frac{\theta}{\tau_{\text{eff}}}, \quad \theta = R \frac{dh}{dR} \approx \frac{2h}{7}, \qquad (3b) $$

$$ Q_{\text{cool}} = -2 \frac{\sigma_{\text{SB}} T^4}{\tau_{\text{eff}}}, \quad \tau_{\text{eff}} = \frac{3\tau_R}{8} + \frac{\sqrt{3}}{4} + \frac{1}{4\tau_p}, \quad \tau_{\text{K,P}} = \frac{\kappa_{\text{K,P}} \Sigma_g}{2}, \qquad (3c) $$

$$ Q_{\text{rad}} = \sqrt{2\pi H \nabla \cdot \left( \lambda \frac{4\sigma_{\text{SB}}}{\kappa_{\text{R}}\rho_{\text{mid}}} \nabla T^4 \right)}, \quad \rho_{\text{mid}} = \frac{1}{\sqrt{2\pi}} \frac{\Sigma_g}{H}. \qquad (3d) $$

In the above, we follow the $\alpha$-viscosity prescription of Shakura & Sunyaev (1973), the irradiation model of Menou & Goodman (2004) for a star with luminosity $L_*$ assuming a disc albedo $\epsilon = 1/2$, an effective optical depth $\tau_{\text{eff}}$ that depends on the Rosseland and Planck mean opacities $\kappa_R$ and $\kappa_P$ according to Hubeny (1990), and the flux-limited diffusion (FLD) closure of Levermore & Pomraning (1981) for the radiative diffusion term. We use the prescription by Kley (1989) for the flux limiter $\lambda$. Finally, $\sigma_{\text{SB}}$ denotes the Stefan-Boltzmann constant. For more details on the individual radiative terms, we refer the reader to Ziampras et al. (2023).

### 2.2 Dust model, parameters, and initial conditions

We consider a solar-type star with $M_* = 1 M_\odot$ and $L_* = 1 L_\odot$ and a disc with surface density profiles $\Sigma_g = 1000 R_{au}^{-1} \text{ g cm}^{-2}$ and $\Sigma_d = 0.01 \Sigma_g$ for the gas and dust, respectively. This translates to a disc mass of $\sim 0.07 M_\odot$ assuming a disc that extends between 0.1 and 100 au. For the dust grains, we assume a two-population model similar to Ziampras et al. (2025), where the total dust mass is distributed between ‘small’ and ‘big’ grains with $a_d^{\text{small}} = 0.1 \mu\text{m}$ and $a_d^{\text{big}} = 1 \text{ mm}$. For convenience, we define a coagulation fraction $X$ and write

$$ X = \Sigma_d^{\text{big}} / \Sigma_d, \quad \Sigma_d = \Sigma_d^{\text{small}} + \Sigma_d^{\text{big}}. \qquad (4) $$

In this way, $X=0$ and $X=1$ correspond to discs containing only small or big grains, respectively. The small grains are assumed to be perfectly coupled to the gas, and at any given time their surface density is given by $\Sigma_d^{\text{small}} = (1 - X_0)\epsilon_0 \Sigma_g$, where $X_0$ and $\epsilon_0 = 0.01$ are the coagulation fraction and dust-to-gas ratio at $t=0$. The big grains are initialized as $\Sigma_d^{\text{big}} = X_0\epsilon_0 \Sigma_g$, and evolved according to equations (1b) and (1d).
<table><caption>**Table 1.** List of hydrodynamical models in this study. Values for X refer to t = 0, and values for $h_0$ refer to t = 0 and R = $R_0$ = 30 au.</caption><thead><tr><th>X</th><th>α</th><th>h<sub>0</sub></th><th>Comment</th></tr></thead><tbody><tr><td>0.9</td><td>10<sup>-5</sup></td><td>0.05</td><td>Fiducial model</td></tr><tr><td>0.1</td><td>10<sup>-5</sup></td><td>0.05</td><td>Fewer millimetre grains (10 per cent by mass)</td></tr><tr><td>0.99</td><td>10<sup>-5</sup></td><td>0.05</td><td>More millimetre grains (99 per cent by mass)</td></tr><tr><td>0.9</td><td>10<sup>-4</sup></td><td>0.05</td><td>Moderately viscous</td></tr><tr><td>0.9</td><td>10<sup>-3</sup></td><td>0.05</td><td>Highly viscous</td></tr><tr><td>0.9</td><td>10<sup>-5</sup></td><td>0.07</td><td>Moderately hotter disc (L<sub>*</sub> = 10.6 L<sub>☉</sub>)</td></tr><tr><td>0.9</td><td>10<sup>-5</sup></td><td>0.1</td><td>Much hotter disc (L<sub>*</sub> = 128.7 L<sub>☉</sub>)</td></tr></tbody></table>

The Rosseland and Planck mean opacity of each grain population is given by

$$ \kappa_{\mathrm{R}}^{\mathrm{small}} \approx 0.27 T_{\mathrm{K}}^{1.6} \mathrm{cm}^2 \mathrm{g}^{-1}, \quad \kappa_{\mathrm{P}}^{\mathrm{small}} \approx 0.41 T_{\mathrm{K}}^{1.6} \mathrm{cm}^2 \mathrm{g}^{-1}, \quad (5) $$

$$ \kappa_{\mathrm{R}}^{\mathrm{big}} \approx \kappa_{\mathrm{P}}^{\mathrm{big}} \approx 4.5 \mathrm{cm}^2 \mathrm{g}^{-1}, \quad (6) $$

following the calculations by Ziampras et al. (2025) with the code OPTOOL (Dominik, Min & Tazaki 2021), using the distribution of hollow spheres approach (Min, Hovenier & de Koter 2005). The dust grains are composed of 87 per cent amorphous pyroxenes and 13 per cent amorphous carbon (Zubko et al. 1996) with a porosity of 25 per cent, resulting in $\bar{\rho}_d = 2.08 \mathrm{g} \mathrm{cm}^{-3}$. The total opacity per gram of gas to be used in equations (3b)–(3d) is then

$$ \kappa_{\mathrm{R,P}} = \frac{\Sigma_{\mathrm{d}}^{\mathrm{small}} \kappa_{\mathrm{R,P}}^{\mathrm{small}} + \Sigma_{\mathrm{d}}^{\mathrm{big}} \kappa_{\mathrm{R,P}}^{\mathrm{big}}}{\Sigma_{\mathrm{g}}} = \varepsilon \left[ (1-X)\kappa_{\mathrm{R,P}}^{\mathrm{small}} + X\kappa_{\mathrm{R,P}}^{\mathrm{big}} \right]. \quad (7) $$

As shown in Ziampras et al. (2025), this approach works well even for marginally optically thick regions where the Rosseland mean opacity becomes important.

For a small $\alpha = 10^{-5}$ and our disc model, equation (1e) yields in steady state $h \approx 0.019 R_{\mathrm{au}}^{2/7}$, or $T \approx 90 R_{\mathrm{au}}^{-3/7}$ K (see e.g. Chiang & Goldreich 1997), regardless of our choice of X due to heating from viscous dissipation being negligible. Translating the above to a reference radius $R_0 = 30$ au, our initial conditions are then

$$ \Sigma_{\mathrm{g,0}} = 33.3 \mathrm{g} \mathrm{cm}^{-2} \left(\frac{R}{R_0}\right)^{-1}, \quad \Sigma_{\mathrm{d,0}}^{\mathrm{big}} = X_0 \varepsilon_0 \Sigma_{\mathrm{g,0}}, \quad (8a) $$

$$ h_0 = 0.05 \left(\frac{R}{R_0}\right)^{2/7} \Rightarrow T_0 = 20.9 \mathrm{K} \left(\frac{R}{R_0}\right)^{-3/7}. \quad (8b) $$

In our models, we choose a fiducial value of X = 0.9 (i.e. 90 per cent of the dust mass is in big grains) but also consider X = {0.1, 0.99} to investigate the effect of dust growth on both radiative processes and the observability of any resulting substructure. We also explore the effects of diffusion by varying $\alpha = \{10^{-5}, 10^{-4}, 10^{-3}\}$, with $10^{-5}$ being our fiducial value.

Finally, we consider three different disc aspect ratios by varying $h(R_0) = \{0.05, 0.07, 0.1\}$, with 0.05 being our fiducial value. To achieve this, we simply adjust the stellar luminosity $L_*$ in equation (3b) to $\{1, 10.6, 128.7\} L_{\odot}$ for $h = \{0.05, 0.07, 0.1\}$, respectively. While this is of course unrealistic, it allows us to investigate both the effect of a higher sound speed on the formation of spirals and the importance of radiative cooling on any resulting substructure, as the cooling time-scale (loosely) scales with $t_{\mathrm{cool}} \propto e/|Q_{\mathrm{cool}}| \propto T^{-3} \propto h^{-6}$.

A list of all models carried out in this work can be found in Table 1.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic of the disc configuration with slices along the x (top) and y (bottom) axes. The inner disc is tilted about the x-axis by 30° with respect to the x-y plane and casts a shadow on the main disc.</figcaption></figure>

## 2.3 Shadows due to a misaligned inner disc

Similar to Montesinos et al. (2016) and Cuello et al. (2019), we consider a small, misaligned inner disc casting a shadow on the ‘main’ disc that we actually evolve in our simulations. We assume that the inner disc extends in the range $R \in [0.2, 4]$ au and is tilted by 30° with respect to the x-y plane about the x-axis, while the main disc extends beyond 5 au. The two discs otherwise share the same physical properties (see Section 2.2), and we apply an exponential taper in the radial direction for both discs for a smooth transition between the two. A detailed description of our two-disc model is provided in Appendix A.

We construct this model in PYTHON using a 3D spherical grid {$r, \theta, \phi$} logarithmically spaced in $r$ with $r \in [0.1, 300]$ au, $\theta \in [0, \pi]$, $\phi \in [-\pi, \pi]$, and $N_r \times N_\theta \times N_\phi = 512 \times 768 \times 512$ cells. A plot of the configuration is shown in Fig. 1.

To measure the width and depth of the shadow cast by the inner disc, we then trace rays from the central star along the radial direction and measure the irradiation heating term $S$:

$$ S \approx \frac{L_* e^{-\tau_{\mathrm{in}}(r)}}{4\pi\sigma^2} \left(1 - e^{-\kappa_{\mathrm{P}}^{\mathrm{small}} \rho_d^{\mathrm{small}} \Delta r}\right), \quad \tau_{\mathrm{irr}}(r) = \int_{0.1 \mathrm{au}}^{r} \kappa_{\mathrm{P}}^{\mathrm{small}} \rho_d^{\mathrm{small}} dr. \quad (9) $$

Here, $\Delta r$ is the radial cell width, and $\rho_d^{\mathrm{small}}$ is the volume density in small dust, assumed to be perfectly coupled to the gas and therefore given by $\rho_d^{\mathrm{small}}(r, \theta) = \varepsilon(1-X)\rho_g$, where $\rho_g$ is the gas density (see equation A4). We set $X = 0.9$ and $\varepsilon = 0.01$.

In Fig. 2, we show a heat map of $S$ along a surface tilted by $4h_{5au} \approx 7^\circ$ with respect to the x-y plane about the x-axis, roughly tracing the irradiation surface of the main disc in order to showcase the geometry of the shadow. From this figure, it becomes clear that inside the shadow, the irradiation heating is reduced to practically zero, while outside the shadow, the heating is unaffected. We can then define a mask that we apply to the irradiation heating term $Q_{\mathrm{irr}}$.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Heat map of the irradiation heating term $S$ through equation (9) along a surface tilted by 7° with respect to the x-y plane, tracing the irradiation surface of the main disc. The shadow cast by the inner disc is clearly visible as a dark wedge. Dashed lines mark the shadow edges, fit by eye.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Azimuthal profile of the irradiation heating term $S$ at $z=0$ and $R=10$ au. The mask $f_{sh}$ used in our disc models, defined in equation (10), is shown with a dashed line.</figcaption></figure>

in equation (3c) to simulate the shadow cast by the inner disc. This mask is compared to the heating term $S$ at the mid-plane in Fig. 3 and is given by

$$f_{sh} = \begin{cases} d_{sh}, & |x| \le 0.7, \\ d_{sh} + (1 - d_{sh}) \sin\left(\frac{|\phi - \phi_e|}{0.3 w_{sh}} \frac{\pi}{2}\right)^2, & 0.7 < |x| < 1, \\ 1, & |x| \ge 1, \end{cases} \quad (10)$$

with $x = (\phi - \phi_0)/w_{sh}$ and $\phi_e = \phi_0 + \text{sgn}(x)/0.7$. The shadow depth and half-width are set to $d_{sh} = 10^{-5}$ and $w_{sh} = 0.24$, respectively. We model two diametrically opposed shadows at $\phi_0 = 0$ and $\phi_0 = \pi$. The mask is then applied to the irradiation heating term $Q_{irr}$ in equation (3c) as

$$Q_{irr} \rightarrow Q_{irr} f_{sh}. \quad (11)$$

We note that we are not interested in the absolute value of $S$ but rather its azimuthal dependence, which we then use to mask the main disc in our simulations. Since the disc model has already been

established in Section 2.2 (i.e. we do not actually calculate the disc thermal structure through ray tracing), our results are only weakly sensitive to the opacity model. For this reason, we use the same mask for all models instead of recomputing it for different $X$. To isolate the effects of different $h$ in the dynamics of the disc, we also use the same mask for all $h$, even though the shadow would realistically be wider for a geometrically thicker disc.

## 2.4 Numerical set-up

We solve the 2D, vertically integrated set of equations (1a)–(1e) using the numerical hydrodynamics code PLUTO (Mignone et al. 2007), combined with the dust module described in Ziampras et al. (2025), the FLD module in Ziampras et al. (2020), and a dust diffusion module following Weber et al. (2019). To relax time-step limitations, we use the FARGO algorithm (Masset 2000), implemented in PLUTO by Mignone et al. (2012). We use a second-order time- and spatially accurate scheme with the HLLC (Harten–Lax–van Leer contact) Riemann solver (Toro, Spruce & Speares 1994) and the flux limiter of Van Leer (1974). Viscosity is implemented using the super-time-stepping method of Alexiades, Amiez & Gremaud (1996).

Our grid extends radially with a logarithmic spacing in the range $R \in [6, 75]$ au $= [0.2, 2.5]$ $R_0$ with a reference radius of $R_0 = 30$ au and covers the full azimuthal extent $\phi \in [0, 2\pi]$ with $N_r \times N_\phi = 576 \times 1536$ cells. This grid resolution translates to approximately 12 cells per gas scale height in either direction at $R = R_0$ for $h_0 = 0.05$. Our disc is initialized according to equations (8a) and (8b), with $u_R = 0$ for both gas and dust and

$$u_{\phi,g} = \Omega_K R \sqrt{1 + \frac{d \log P}{d \log R} h^2}, \quad u_{\phi,d} = \Omega_K R, \quad (12)$$

to correct for pressure support.

We impose wave-damping boundary conditions at both radial boundaries on all quantities following de Val-Borro et al. (2006) with a damping time-scale of 0.1 boundary orbits. At the boundary walls, all quantities are fixed to their initial values. We then evolve the system for 1000 orbits at $R_0 (= 1000 P_0)$, which corresponds to approximately 164 kyr ($P_0 \approx 164$ yr).

# 3 FIDUCIAL MODEL

In this section, we present the results of our fiducial numerical simulation with $X = 0.9$, $\alpha = 10^{-5}$, and $h(R_0) = 0.05$. This set of parameters comprises a reasonably realistic disc model with a significant fraction of the dust mass in big grains after substantial dust growth has taken place (see e.g. Birnstiel 2024), a small viscosity motivated by the lack of magnetically driven turbulence (e.g. Bai & Stone 2013), and a temperature profile that corresponds to a passively irradiated disc (Chiang & Goldreich 1997).

## 3.1 Disc structure

We start by showing heat maps of the perturbed temperature and surface density of the disc after 1000 $P_0$ of Fig. 4. The edges of the two shadows are marked with dotted lines in the perturbed temperature (left panel), with the gas cooling rapidly once it enters a shadow and then heating back up to the irradiation temperature upon exiting. The temperature troughs are azimuthally offset with respect to the centre of each shadow as the gas responds to heating and cooling over a finite time-scale.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The perturbed gas temperature (left) and gas surface density (middle) as well as the dust surface density (right) after 1000 $P_0$. Two dark wedges are visible in the temperature structure, corresponding to the attenuated heating due to the ‘shadows’ cast by an inner disc. A series of concentric rings is visible in all quantities, forming clearly visible features in the dust surface density. The shadow edges are marked with dashed lines in the left panel. The disc rotates counterclockwise.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Azimuthal profile of the perturbed gas temperature at different radii at $t = 1000 P_0$. The temperature perturbation is not symmetric about the shadow exits and becomes smaller in the slowly cooling inner disc. Shaded regions mark the two shadows. Gas moves to the right.</figcaption></figure>

The azimuthal profile of temperature at various radii is shown in Fig. 5, further revealing that the temperature perturbation is not symmetric about the exit of the shadow. This happens due to radiative cooling being proportional to $T^4$, and results in the gas cooling more rapidly than it heats up. We can also see that the gas does not quite reach the irradiation temperature before entering the next shadow, leading to a net cooling of the disc by a few per cent compared to the initial state. This effect is more pronounced in the inner disc, where the cooling time-scale is longer due to the higher gas density and opacity. The slower cooling also results in a smaller temperature contrast between the shadowed and illuminated regions. Specifically, as the cooling time-scale increases, the azimuthal temperature perturbation becomes smaller and approaches a value of ≈4 per cent, reflecting the fact that the disc is irradiated less on average due to the pair of shadows.

Interestingly, while expecting a spiral structure in the gas surface density (middle panel of Fig. 4) as the result of the radial force imbalance due to the shadow, we instead find a series of concentric rings. Even though the contrast of these rings is of the order of

5 per cent in the gas surface density, the local modulation of the pressure gradient is enough to significantly enhance the dust-to-gas ratio at the peak of each ring over time, as shown in the dust surface density (right panel).

We speculate that the formation of rings is inherently linked to the launching of spiral waves at the shadow edges, due to the lower temperature in the shadowed region resulting in a radial force imbalance (see equation 12). If these spiral waves can dissipate into the disc, the resulting AMF can drive the formation of gaps in the gas surface density. This, in turn, leads to a smaller radial pressure gradient between two gaps, potentially trapping – or at least slowing down – dust grains. The result is a series of rings in the dust continuum emission.

Fig. 6 shows the perturbed gas surface density with respect to the initial state $\Sigma_{g,0}$ (top panels) and the azimuthally averaged state $\bar{\Sigma}_g$ (bottom panels) at different snapshots. At early times, an $m = 2$ spiral pattern can be seen in both top and bottom panels, owing to the radial force imbalance in the shadow perturbing the initially axisymmetric disc. This spiral pattern persists in the bottom panels, indicating that spirals are continuously generated throughout the disc. In the top panels, however, we observe a transition from spiral- to ring-like structures, with the rings becoming more pronounced over time. The transition seems to take place after about 120 $P_0$ (second column), where we can see what looks like spirals breaking up before merging into rings.

This behaviour supports our hypothesis that the constant generation of spiral waves due to the shadows drives an AMF into the disc once the waves dissipate, ultimately leading to the formation of rings. We investigate further by computing the AMF in the disc, defined as (e.g. Miranda & Rafikov 2020a)

$$F_J = \oint_{\phi} \delta F_J d\phi, \quad \delta F_J(R, \phi) = \bar{\Sigma} R^2 (u_R - \bar{u}_R)(u_\phi - \bar{u}_\phi), \quad (13)$$

where all quantities refer to the gas and bars denote azimuthal averaging. We then plot azimuthally averaged profiles of the perturbed gas and dust surface densities as well as the AMF $F_J$ as a function of time in Fig. 7. This figure confirms that an AMF is indeed present
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The perturbed gas surface density with respect to the initial state $\Sigma_{g,0}$ (top panels) and the azimuthally averaged state $\bar{\Sigma}_g$ (bottom panels) at different snapshots. The disc initially shows spiral-like structures, but transitions to the rings visible in Fig. 4 after about $120 P_0$. The undulations seem to be more pronounced at the exits of the shadows and increase in amplitude over time.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Azimuthally averaged profiles of the perturbed gas and dust surface densities (panels a and b), the AMF $F_J$ (panel c), and the cooling time-scale $\beta$ (panel d) as a function of radius and for several snapshots. The AMF is approximately constant over time, leading to deeper gaps and sharper features as the disc evolves. The disc is roughly separated into a featureless inner region and a ringed outer region at $R \approx 24$ au, which corresponds to $\beta \approx \beta_{crit}$ (see equation 15). The AMF is arbitrarily normalized to $F_{J,0} = 10^{-10} M_{\odot} \text{au}^2 \text{yr}^{-2}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Time evolution of the azimuthally averaged perturbed gas and dust surface densities. Gaps deepen and dust accumulates in rings roughly linearly over time.</figcaption></figure>

in the disc, with the troughs in the gas surface density corresponding to peaks in the AMF. Dust features then form at the radial locations where the radial pressure gradient is minimized, slightly offset from the peaks in gas surface density due to the radial temperature gradient.

The fact that $F_J$ remains roughly constant over time also suggests that the gap depth will approximately increase linearly with time, resulting in sharper features as the disc evolves. We verify this statement by plotting the azimuthally averaged perturbed gas and dust densities as a function of time in Fig. 8, where we show that the undulations in both quantities increase in amplitude over time.

We also plot a heat map of the two-dimensional quantity $\delta F_J$ averaged over 200 snapshots in the range $t \in [100, 300] P_0$ in Fig. 9.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Heat map of the quantity $\delta F_J$ averaged over 200 snapshots in the range $t \in [100, 300]$ $P_0$. A high AMF is present at the exits of the shadows, indicating the deposition of angular momentum into the disc and leading to the sustained formation of gaps. Similar to Fig. 7, $F_{J,0} = 10^{-10}$ M$_\odot$ au$^2$ yr$^{-2}$.</figcaption></figure>

We find that the AMF is strongest at the exit of either shadow, which also implies that gaps should be slightly deeper at those locations compared to the illuminated regions, a quarter of an orbit later. This exact behaviour can be seen in the top panels of Fig. 6, with perturbations being strongest at the exits of the shadows and weakest after an angular distance of $\pi/2$.

## 3.2 A criterion for the formation of rings

The dependence of the azimuthal temperature perturbation on the cooling time-scale (see Fig. 5) might hint at the reason why we do not observe rings or spirals in the inner disc. In particular, smaller azimuthal temperature perturbations would result in weaker spiral waves and therefore a smaller AMF, in turn leading to less pronounced (if at all present) rings. By combining the contributions to cooling due to thermal emission from the disc surfaces and in-plane radiative diffusion, we define the cooling time-scale $\beta$ following Ziampras et al. (2025) as

$$ \beta = \frac{1}{f+1} \frac{e}{|Q_{\text{cool}}|} \Omega_{\text{K}}, \quad f = 16\pi \frac{\tau_{\text{eff}}}{\tau_{\text{R}}} \frac{\tau_{\text{p}}^2}{6\tau_{\text{p}}^2 + \pi}. \qquad (14) $$

We then plot the azimuthally averaged cooling time-scale $\bar{\beta}$ at $t=0$ in the bottom panel of Fig. 7, showing that it decreases with radius as expected.

To crudely estimate the conditions under which temperature perturbations due to the shadows are too small to efficiently generate spiral waves, we require that the cooling time-scale $\beta$ becomes longer than the time it takes to cross the shadowed regions (see also Casassus et al. 2019). Assuming Keplerian rotation, and since we model two shadows, this critical time-scale is then

$$ \beta_{\text{crit}} \sim 2 \frac{\Delta t_{\text{shadow}}}{\Omega_{\text{K}}} \approx 2 \frac{2w_{\text{sh}}\Omega_{\text{K}}}{\Omega_{\text{K}}} = 4w_{\text{sh}}, \qquad (15) $$

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Top: Radial profile of the perturbed gas surface density at $t = \{100, 30, 15\}$ $P_0$ for $h = \{0.05, 0.07, 0.1\}$, respectively. Dots mark the location of peaks. Bottom: Radial spacing between rings in the gas surface density as a function of radius. The dotted line is based on equation (17) for $f = 2/7$.</figcaption></figure>

or $\beta_{\text{crit}} \approx 0.96$ for our models, which translates to $R \approx 24$ au. We mark this radial location with a vertical dashed line in Fig. 7, showing that it roughly divides the disc into a feature-poor inner region and a ringed outer region, consistent with our expectations.

We note that this is only a rough estimate in the interest of providing a simple explanation for the lack of features in the inner disc. In principle, a weak AMF will exist even for small temperature perturbations due to a shadow, and rings would still form after an arbitrarily long time, which might of course be longer than the lifetime of the disc.

## 3.3 Spacing between rings

Figs 4 and 7 hint at a regular spacing between rings in the gas surface density. We compute this spacing as the radial distance between two adjacent rings normalized to the radius of the interior ring in Fig. 10, confirming that the rings are fairly regularly spaced in radius.

Given that the disc exhibits two large-scale, continuous spirals (see bottom panels of Fig. 6) due to the two shadows, and considering that the AMF peaks at the shadow exits, we can make a rough estimate of the radial spacing between the two rings. To do this, we define the pitch angle $\phi_p$ of the spiral waves as

$$ \tan \phi_p = \frac{1}{R} \frac{\partial R}{\partial \phi} = \frac{1}{R} \frac{\partial R}{\partial t} \frac{\partial t}{\partial \phi} \approx \frac{\sqrt{\gamma c_s}}{R \Omega_K} = \sqrt{\gamma h}, \qquad (16) $$

assuming that the waves travel at the adiabatic sound speed $\sqrt{\gamma c_s}$. With $h = h_0 (R/R_0)^t$, we then integrate equation (16) from one shadow exit to the next to obtain

$$ \int_{R_0}^{R_0+\Delta R} \frac{1}{hR} dR = \sqrt{\gamma} \int_{\phi_0}^{\phi_0+\pi} d\phi \Rightarrow \frac{\Delta R}{R_0} = (1 - \sqrt{\gamma \pi f h_0})^{-1/f} - 1. \qquad (17) $$
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Heat maps similar to Fig. 4 for coagulation fractions $X = \{0.1, 0.9, 0.99\}$ at $t = 1000 P_0$. The $X = 0.1$ model is dominated by small grains and is therefore optically thicker, resulting in rings forming slightly further out due to the cooling criterion in equation (15). The rings are also much fainter in this model due to the smaller fraction of big grains.</figcaption></figure>

We then plot this result in Fig. 10 with a dotted line for $f = 2/7$ and all values of $h$ (discussed further in Section 4.3), showing that this simple model captures the spacing between rings fairly well, aside from the left and right edges of the disc due to radial boundary effects.

# 4 PARAMETER STUDY

Having analysed the fiducial model, we explore here the effects of varying the coagulation fraction $X$, the viscosity $\alpha$, and the aspect ratio $h$.

## 4.1 Coagulation fraction $X$

Fig. 11 shows heat maps of the perturbed gas density and the dust density similar to Fig. 4 for the three different values of the coagulation fraction $X \in \{0.1, 0.9, 0.99\}$ used in our study. Here, we find that the effect of having a significant fraction of the dust mass as small grains ($X = 0.1$) is twofold. For one, it results in rather faint and less sharp features in the dust distribution, as the majority of the dust is well coupled to the gas and does not accumulate on the forming pressure bumps.

At the same time, the higher opacity of the submicron-sized grains renders a larger portion of the inner disc optically thick, resulting in slower cooling and ultimately limiting the formation of rings in the gas (and therefore the dust) structure. To verify this statement, we show a radial profile of the cooling time-scale at $t = 0$ for all three models following equation (14) in Fig. 12. Through this figure, we can also explain why the models with $X \ge 0.9$ look quite similar: the cooling time-scale is quite similar between the two for $R \gtrsim 20$ au.

Nevertheless, for realistic values of $X \gtrsim 0.9$ that account for the dust growth happening in a typical protoplanetary disc by its T Tauri stage (e.g. Birnstiel 2024), we find that cooling is sufficiently fast to facilitate the formation of rings in the bulk of the disc.

## 4.2 Turbulent $\alpha$

The effect of modelling turbulence via a viscous $\alpha$ stress is rather straightforward, acting against the formation of radial structure by diffusing the ring-like gas features and related pressure bumps. It

<figure><img src="image_14.png" /><figcaption>**Figure 12.** The azimuthally averaged, density-weighted cooling time-scale following equation (14) at $t = 0$ for the three models with different $X$. A lower $X$ results in more small grains and an optically thicker, slowly cooling disc. A dot marks the radial location where $\beta \approx \beta_{crit}$ according to equation (15).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Similar to Fig. 11, but for different viscosity parameters $\alpha = \{10^{-5}, 10^{-4}, 10^{-3}\}$. Rings are slightly more diffuse for $\alpha = 10^{-4}$ and hardly visible for $\alpha = 10^{-3}$ when compared to the fiducial model with $\alpha = 10^{-5}$.</figcaption></figure>

is nevertheless useful to check how strong turbulence should be to eliminate (or at least strongly damp) the substructures we find in our fiducial model.

In Fig. 13, we show our results for models with $\alpha = 10^{-4}$ and $10^{-3}$, complementing our fiducial value of $10^{-5}$. While the model with $\alpha = 10^{-4}$ shows slightly more diffuse rings with a weaker contrast compared to our fiducial case, the radial structure is still clearly visible. This is not the case for $\alpha = 10^{-3}$, however, where the viscous torques manage to wash away any developing radial features.

Overall, for typical values of $\alpha \lesssim 10^{-4}$ motivated by hydrodynamical mechanisms such as the vertical shear instability (e.g. Nelson, Gressel & Umurhan 2013; Flock et al. 2020) and the lack of efficient magnetohydrodynamic (MHD) turbulence in the ‘magnetically dead’ zone of the disc (Bai & Stone 2013), we expect that turbulence should not inhibit significantly the formation of rings through the mechanism studied in this paper.

## 4.3 Aspect ratio $h$

Given that the aspect ratio $h$ largely determines the width of the shadow through the scale height of the misaligned inner disc (see Fig. 1), the cooling time-scale as $\beta \propto h^{-6}$ (see equation 14), and the radial spacing between rings (equation 17), it is very interesting to
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Top: The perturbed gas surface density for different aspect ratios, marked in the corner of each panel. The leftmost panel corresponds to $t = 50 P_0$, while the others refer to $t = 1000 P_0$. Bottom: The dust surface density for the same models. The significantly shorter cooling time-scale for larger $h$ drives a rapid formation of deep gaps, which often result in the nearby rings collapsing into vortices due to the RWI.</figcaption></figure>

investigate how the disc structure changes with $h$ even though we hold the shadow width $w_{sh}$ fixed. We show a comparison among models with $h(R_0) = \{0.05, 0.07, 0.1\}$ in Fig. 14.

Regarding the formation of rings, a higher $h$ results in a significantly shorter cooling time-scale, which in turn leads to a higher temperature contrast between the shadowed and illuminated regions (see Fig. 5), and therefore more vigorous spiral waves, a higher AMF into the disc, and ultimately more efficient gap opening. This is clearly visible on the left panels of Fig. 14, which show snapshots at $t = 50 P_0$ and illustrate the rapid formation of structure in both gas and dust for higher $h$. In particular, gap opening is so efficient for $h = 0.1$ that the gap edges are already unstable to the Rossby wave instability (RWI; Lovelace et al. 1999) after only $50 P_0$, resulting in the formation of vortices in the gas surface density and bright clumps in the dust.

Within $1000 P_0$, the discs with $h(R_0) \ge 0.07$ destabilize to the RWI throughout, developing a series of vortices that span the full radial extent of the disc. These vortices decay as they orbit, in part due to viscous diffusion, cooling (Fung & Ono 2021; Rometsch et al. 2021), and dust-gas coupling (Raettig, Klahr & Lyra 2015; Lovascio, Paardekooper & McNally 2022; Ziampras et al. 2025). For $h = 0.1$, several vortices in the outer disc ($R \gtrsim 20$ au) merge into a single large vortex, while most of the vortices that formed in the inner disc eventually dissipate along the azimuthal direction, resulting in the dust structure reverting to a set of rings.

The vortex-rich structure seen in models with $h(R_0) \ge 0.07$ is in stark contrast to the fiducial model with $h = 0.05$, which shows a purely radial structure. We highlight, however, that the gap opening process is still ongoing and will inevitably lead to gap edges steep enough to be unstable to the RWI, even for this model. Whether this will happen within the lifetime of the disc is beyond the scope of this study.

# 5 SYNTHETIC OBSERVATIONS

In this section, we test whether the features found in our numerical models could be observable with current observational facilities such as the Atacama Large Millimeter/submillimeter Array (ALMA) and the Very Large Telescope (VLT). To that end, we show synthetic observations of the disc models presented in the previous section using the radiative transfer code RADMC-3D (Dullemond et al. 2012). We focus on the dust continuum emission at 1.3 mm (231 GHz, corresponding to ALMA Band 6) and the scattered-light emission at 1.65 µm [corresponding to VLT-Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) in H band] to examine the disc structure of the mid-plane and surface layers of the disc, respectively.

## 5.1 RADMC-3D set-up

Since we operate in a vertically integrated framework for our suite of hydrodynamical simulations, we construct a 3D model of the disc assuming a vertically isothermal disc in hydrostatic equilibrium. We follow a similar procedure as in Appendix A, this time replacing the surface density profile of the outer disc with that from our hydrodynamical models within the radial range $R \in [0.2, 2.5] R_0 = [6, 75]$ au. The vertical dust profile for both big and small grains is computed assuming a vertically constant turbulent mixing parameter $\alpha_z$ following Fromang & Nelson (2009):

$$ \rho_d^i = \frac{\Sigma_d^i}{\sqrt{2\pi H}} \sqrt{\frac{St_{\text{mid}}^i}{\alpha_z} + 1} \exp \left\{ -\frac{St_{\text{mid}}^i}{\alpha_z} \left[ \exp\left(\frac{z^2}{2H^2}\right) - 1 \right] - \frac{z^2}{2H^2} \right\}, \quad (18) $$

where $St_{\text{mid}}^i$ is the Stokes number of species with index $i$ evaluated at the mid-plane through equation (2), and $H$ is the pressure scale height of the gas. The wavelength-dependent opacity table for each
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Left: Vertical cross-section of the grid used in RADMC-3D models. Each cell in the figure corresponds to a block of cells in the grid used, with the colour indicating the number of cells in the block. Right: Dust density structure for $h_0 = 0.05$. The bottom row zooms into the box in the top row.</figcaption></figure>

grain species is computed using OPTOOL using the composition in Section 2.2.

Even though the hydrodynamical simulations processed with RADMC-3D correspond to runs with $\alpha = 10^{-5}$, we choose $\alpha_z = 10^{-4}$ for our RADMC-3D models. This choice is motivated by the fact that turbulent diffusion does not affect the ring formation process significantly at $\alpha = 10^{-4}$ (see Fig. 13), while a higher $\alpha_z$ helps recover the temperature structure of the disc without the need for an extremely large number of photons.

To resolve the highly settled layer of millimetre grains in the polar direction for both the inner misaligned disc and the outer disc, we use a grid with

(i) 256 cells in the range $\theta \in [(\pi/2), (\pi/2) + h_0]$,

(ii) 64 cells in the range $\theta \in [(\pi/2) + h_0, \theta_{\text{tilt}} - (h_0/2)]$, $\theta_{\text{tilt}} = (\pi/2) + 30^{\circ}$,

(iii) 128 cells in the range $\theta \in [\theta_{\text{tilt}} - (h_0/2), \theta_{\text{tilt}} + (h_0/2)]$, and

(iv) 16 cells in the range $\theta \in [\theta_{\text{tilt}} + (h_0/2), \pi]$,

with a mirror symmetry about the mid-plane ($\theta = \pi/2$) for a total of $N_{\theta} = 920$ cells. The radial extent covers $r \in [0.1, 110]$ au with $N_r = 512$ logarithmically spaced cells, and the azimuthal domain spans $\phi \in [0, 2\pi]$ with $N_{\phi} = 514$ cells. A cross-section of the grid used in RADMC-3D is shown in Fig. 15. The dust surface density from the hydrodynamical models is first interpolated onto the mid-plane of the 3D grid ($r, \theta = \pi/2, \phi$), and the 3D density structure for each species is then computed using equation (18).

To calculate the dust temperature with RADMC-3D we use the mctherm command with $10^{10}$ photons and a solar-type star ($M_* = 1 M_{\odot}$, $T_* = 5772$ K). To make sure that the frequency range for stellar emission is the same for all values of $h_0$, we change the stellar radius to $R_* = [1, 3.25, 11.34] R_{\odot}$ such that $L_* = [1, 10.6, 128.7] L_{\odot}$ for $h_0 = [0.05, 0.07, 0.1]$, respectively, and treat the star as a point source to avoid unexpected geometrical effects. This is a good approximation for our set-up, as $R \gg R_{\odot}$ even for the inner disc. The vertical structure of the misaligned inner disc is also now computed with the appropriate $h_0$ such that its shadow can be captured in

the dust temperature calculation. Finally, we include the effects of anisotropic scattering using the approximation of Henyey & Greenstein (1941) with the scattering asymmetry parameter $g$ computed using OPTOOL.

It is important to note here that the temperatures computed by RADMC-3D are equilibrium temperatures, in contrast to those used in the dynamic calculations. As a result, RADMC-3D cannot capture the asymmetric azimuthal variation of temperature found in the hydrodynamical models around the shadow. We highlight this difference in Appendix C. Given that the azimuthally averaged mid-plane temperature profiles obtained from RADMC-3D are remarkably similar to those from the hydrodynamical models, we replace the RADMC-3D output of the temperature structure within $\pm H/2$ about the mid-plane with the hydrodynamical temperature structure. In doing so, our models can also capture the azimuthal variations in temperature (see e.g. Fig. 5) when computing synthetic images in millimetre emission. We believe that applying this correction to the disc surface would not necessarily produce more accurate results, as our models are 2D, but instead refer the reader to Zhang & Zhu (2024), who computed scattered light images for their (3D) models.

We then compute images at 1.3 mm and 1.65 µm with the image command for various orientations using $10^8$ photons and 768 pixels on both x and y directions, assuming a distance of 100 pc to the source. For images at 1.65 µm, the flux within a circle of radius 100 milliarcsec (mas) from the star is blocked to emulate the use of a coronagraph. The images are then convolved with a Gaussian beam of 40 mas for both wavelengths to mimic typical observations with ALMA Band 6 and SPHERE, respectively.

## 5.2 Synthetic images

In the top two panels of Fig. 16, we present the synthetic images for our fiducial model ($h_0 = 0.05$, $X = 0.9$, $\alpha = 10^{-5}$, $\alpha_z = 10^{-4}$) after 1000 $P_0$. The left panels show the dust continuum emission at 1.3 mm assuming a face-on disc, while the remaining panels show the scattered-light emission at 1.65 µm for different inclinations and azimuthal rotations of the disc with respect to the observer.

In the second row of Fig. 16, which shows the full disc structure, the dust continuum emission features both a series of bright rings as well as a faint azimuthal variation in brightness due to the shadow cast by the misaligned inner disc. In the remaining panels of the same row, corresponding to the scattered-light images, the shadowed wedge is clearly visible, but the rings are not. The top row of the same figure shows zoomed-in views of the inner disc regions, highlighting the directly illuminated inner rim of the outer disc ($R \lesssim 18$ au). However, the structure of this inner region is blurred due to its angular size being comparable to the beam size.

In the third and fourth row of Fig. 16, we show the synthetic images for models with $h_0 = 0.07$ and $h_0 = 0.1$, respectively, similar to the second row of the same figure. As expected, the vortices that form in both models are clearly visible as bright arcs in millimetre emission, and can still be faintly observed in scattered light, but the rings are too faint to be observable. This is a result of the very weak contrast in µm-sized grains, which largely trace the gas and therefore show $\Delta\Sigma/\Sigma \lesssim 5$ per cent. The shadows grow progressively wider with increasing $h_0$ due to the larger scale height of the inner disc, and a dark lane can even be observed in millimetre emission due to the stark temperature contrast between the shadowed and illuminated regions.

Overall, we find that the features in our hydrodynamical models should be observable with current facilities, and especially so for the models with larger aspect ratios, as they both exhibit more pro-
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Synthetic images for our fiducial model with $h_0 = 0.05$ (top two rows) as well as $h_0 = 0.07$ and $h_0 = 0.1$ (third and fourth row) at 1.3 mm (leftmost column) and 1.65 µm (remaining panels) for different disc inclinations and azimuthal rotations. The top row zooms into the (unresolved) inner region at $R \lesssim 18$ au, marked by a box in each panel of the second row. Rings and shadows are visible throughout the entire disc for all values of $h_0$. For better visibility on the same scale, the fluxes are rescaled by a factor shown in the bottom right corner of each panel. A solid white circle marks the beam size of 40 mas, and a dashed white circle the coronagraph of radius 100 mas. A dashed green circle marks the outer edge of our hydrodynamical domain.</figcaption></figure>

nounced substructure and represent hotter, brighter discs. For the case of our fiducial model for $h = 0.05$, however, the features are fainter and might be harder to detect. For this reason, we further process the (unconvolved) images at 1.3 mm using the SIMIO-CONTINUUM tool (Kurtovic 2024)<fnref n="1" /> to simulate an ALMA observation that more

accurately represents the noise induced by the synthesized ALMA beam. We use a configuration that corresponds to the observation of a source analogous to the system HD 163296 ($d = 101$ pc, $i = 47°$, and position angle (PA) = 133°; see e.g. Huang et al. 2018).

Fig. 17 shows the synthetic ALMA image for our fiducial model with $h_0 = 0.05$, with the left panel corresponding to $t = 1000$ $P_0$ as in all panels shown so far in this section. While the contrast is rather weak, rings and gaps are visible to the naked eye without requiring

<footnote n="1">¹SIMIO-CONTINUUM utilizes the CASA package (CASA Team 2022).</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Synthetic ALMA image at 1.3 mm for our fiducial model (see second row of Fig. 16), configured to mimic an observation of the system HD 163296 using SIMIO-CONTINUUM. Left: $t = 1000 P_0 \approx 0.16$ Myr. While faint rings and gaps are visible by eye. Right: $t = 2000 P_0 \approx 0.32$ Myr. The rings are now clearly visible as gap opening progresses. Colours are rescaled with an arcsinh stretch.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Same as Fig. 17, but for models with $h_0 = 0.07$ (left) and $h_0 = 0.1$ (right) at $t = 1000 P_0$. Vortices are clearly visible, but fainter features such as rings and arcs corotating with the vortices are less discernible by eye.</figcaption></figure>

any post-processing. The shadow cast by the misaligned inner disc is, however, not discernible by eye. To highlight the time-sensitive nature of the gap opening process, we also show the synthetic ALMA image at $t = 2000 P_0$ in the right panel of the same figure, where the rings are now clearly visible as gap opening progresses and more dust accumulates in the resulting pressure bumps.

For completeness, we also show the synthetic ALMA images for the models with $h_0 = 0.07$ and $h_0 = 0.1$ in Fig. 18. In these cases, vortices are clearly visible for both models, but the arcs of dust corotating with the vortices for $h_0 = 0.07$ as well as the rings interior to the vortex for $h_0 = 0.1$ are no longer visible by eye. For $h_0 = 0.1$, a dark lane due to the shadow cast by the misaligned inner disc is clearly visible in the synthetic ALMA image, as the shadow is wider and more pronounced for larger aspect ratios.

# 6 DISCUSSION

In this section, we discuss the implications of our results in the context of planet formation and observations of protoplanetary discs. We also address additional effects that we have not considered in our study, but can complicate the picture.

## 6.1 Implications for planet formation

Recent years have shown that substructures in protoplanetary discs appear to be ubiquitous (Andrews et al. 2018; Long et al. 2018). They are of crucial importance for planet formation models (Guilera et al. 2020; Morbidelli 2020), as they are trapping pebbles and thereby facilitating the growth of planetesimals towards planets (Lau et al. 2022; Jiang & Ormel 2023). Planets, however, are one of the most popular mechanisms to explain the observed substructures as they are indicated by kinematics (e.g. Pinte et al. 2018; Teague et al. 2018; Teague, Bae & Bergin 2019) or in exceptional cases even directly detected (Keppler et al. 2018). Planets and substructures may, in fact, trigger each other, as demonstrated in Lau et al. (2024). It is therefore crucial to understand which effects can lead to the formation and evolution of substructures in discs and how to distinguish between the various mechanisms.

In this work, the shadows cast by an inclined inner disc lead to the formation of multiple rings. This may play a significant role in this process as it can seed an outer disc with structures that kick-start planet formation through dust accumulation. However, with numerous rings, each may only trap a limited amount of dust, potentially hindering efficient planet formation if the material is too spread out. The impact of these shadow-induced rings likely evolves over time, as the disc cools and the dust population changes. The demonstrated formation of vortices as the final stage of the ring evolution can further accumulate dust also along the azimuthal direction. This time, evolution could help explain why we observe relatively few substructures in younger Class 0/I discs, but far more ringed structures in later stages of disc evolution.

## 6.2 Effect of a moving shadow

In our models, we assume that the shadow is static in time, or equivalently that the misaligned inner disc is not precessing. This choice was made for the sake of simplicity, although misaligned inner discs are expected to precess (e.g. Facchini, Juhász & Lodato 2018; Kuffmeier et al. 2021), with a common explanation being an inclined planet (e.g. Nealon et al. 2018) or binary companion (e.g. Facchini et al. 2018). A shadow due to a warp in the inner region of the disc (e.g. Papaloizou & Pringle 1983) would also precess (Kimmig & Dullemond 2024). If the warp propagates too fast to be eliminated by viscous damping, it could result in a broken inner disc (Doğan et al. 2015), which would then precess as well.

In the case of a precessing shadow, Montesinos & Cuello (2018) showed that the emerging spiral structure resembles planet-like signatures, with the radial location of the spiral-launching point being a function of the precession direction and frequency of the shadow. A radial profile of the cooling time-scale complicates the picture, as forming radial substructure would require an appropriate balance between the cooling, shadow precession, and orbital/shadow crossing time-scales.

## 6.3 Comparison to related studies

Cuello et al. (2019) were the first to investigate the trapping of dust particles in spirals driven by shadows in a framework similar to ours, but their models did not feature rings or gaps. Instead, the spiral structures presented in their work are much brighter in 1.3 mm emission, contrary to our models, where the spirals are quickly overshadowed by the forming rings or vortices. All of the above can be attributed to several different reasons.
For one, the omission of in-plane radiative diffusion in Cuello et al. (2019) renders the disc less radiatively efficient, as the in-plane cooling channel (which can be much more efficient in the presence of temperature gradients) is unavailable (Miranda & Rafikov 2020b; Ziampras et al. 2023). This results in a longer cooling time-scale and therefore slower development of the gap opening process. In combination with their significantly shorter integration time-scale of $10^4$ yr – 10 orbits at 100 au, their reference radius – it is not surprising that the only prominent feature in their models is the spiral structure, which develops very early on.

Furthermore, the brightness of the spirals in millimetre emission, which primarily depends on the local dust-to-gas ratio, heavily depends on the background dust model. Cuello et al. (2019) used a dust size distribution with $a_d^{\text{max}} = 10$ cm, or $St^{\text{max}} \sim 2$ at 100 au, such that a substantial fraction of dust grains lies in the regime $St \approx 0.1–1$ and can very efficiently accumulate around the spiral arms. The fact that the dust particles were treated using post-processing can further exaggerate this effect, as dust-gas dynamics are missed. In contrast, our models feature millimetre grains with $a_d^{\text{big}} = 1$ mm or $St \approx 0.01$ at $R_0 = 30$ au, which do not experience such efficient accumulation around spirals but will nevertheless be trapped in the resulting pressure bumps.

Nevertheless, the studies by both Montesinos et al. (2016) and Cuello et al. (2019) laid the groundwork in probing the observational signatures of spirals induced in shadowed discs, and in particular, our results relate to Cuello et al. (2019), where the central focus was the dynamics of the dust component. Our results, while seemingly qualitatively different, are in fact fully consistent with theirs once the differences in parameter space and model complexity are taken into consideration.

During the preparation of this paper, we became aware of the work by Su & Bai (2024, hereafter SB24), who also studied the dynamical evolution and formation of substructure in shadowed protoplanetary discs. SB24 opted for a more controlled, statistical approach by using a simplified, constant cooling time-scale $\beta$ per model and carrying out a large set of simulations for different $\alpha$, $w_{\text{sh}}$, $h$, and $\beta$. They then commented on the existence of different ‘branches’ in the parameter space, with the formation of rings being more likely for discs with $\beta \sim 1$ and a weak temperature contrast in the shadow $|\Delta T|/T \lesssim 0.15$, and vortices forming for $\beta = 10^{-3}$.

Our results are entirely consistent with their findings and provide further insight into the physical mechanisms behind the formation of rings and vortices. For a realistic temperature contrast $|\Delta T|/T \sim 0.15$ for $h_0 = 0.05$, all models are firmly in the non-linear regime where spirals can deposit angular momentum into the disc. It then becomes a question of whether the AMF is sufficiently high to open gaps, a process mitigated by a high $\alpha$ or limited by the simulation runtime (or, more realistically, the disc lifetime). For sufficiently low $\alpha$, the gap edges are eventually unstable to the RWI, triggering the formation of vortices. In that regard, the three branches (‘spirals’, ‘rings’, and ‘vortices’) presented in SB24 can be understood as different snapshots of the same process.

We also note that their quoted value of $\beta \sim 1$ needed to form rings is entirely consistent with our criterion in equation (15), which for their models evaluates to $\beta_{\text{crit}} \approx 4\sigma_{\phi} \approx 1$. Given that our models already show vortices by the end of the simulation runtime for $h_0 = 0.07$, where $\beta \lesssim 0.1$ at $R \gtrsim R_0$, it makes sense that the disc destabilizes to the RWI nearly immediately for $\beta = 10^{-3}$ in their models and further suggests that a cooling time-scale of $\beta \lesssim 0.1$ should be sufficient to trigger the RWI well within the disc lifetime under the constant thermal forcing of a shadow.

In addition, recently Zhang & Zhu (2024) showed with three-dimensional radiative hydrodynamics simulations that the spiral-driving mechanism due to a shadow is robust and can excite substantial angular momentum transport as well as bright features in scattered-light observations. While their work focused on the spiral formation process and its observable signatures, it is worth pointing out that their models also featured gaps and vortices (see fig. 2 therein) due to the same mechanism studied in this work. Combined with the recent linear theory calculations by Zhu, Zhang & Johnson (2025), it is clear that shadow-induced spirals can very efficiently transport angular momentum and drive substructure in protoplanetary discs.

All in all, our findings are in very good agreement with previous and current studies, with a key difference being the inclusion of dust hydrodynamics. In this context, the latter enables more consistent dust-gas interaction both via momentum exchange and the determination of the opacity, as well as synthetic observations of millimetre emission with ALMA. Further work is needed to understand dust-gas dynamics in a three-dimensional model, given that the excitation of spirals can drive substantial vertical motion (Zhang & Zhu 2024), which can set the dust scale height and therefore overall cooling efficiency.

## 6.4 Observational prospects

In Section 5, we showed that the rings formed due to the shadowing by a misaligned inner disc should be visible with modern observational facilities such as ALMA and VLT-SPHERE. At the same time, however, the shadow itself might not be detectable for discs with a larger aspect ratio (i.e. hotter). It is therefore possible that several observed discs may exhibit substructure due to a shadow, but the shadow itself might not be visible. Nevertheless, systems with both a visible shadow and substructure have been observed, making them ideal targets for further study. Such systems include HD 100453 (Benisty et al. 2017), HD 142527 (Avenhaus et al. 2017), and HD 143006 (Andrews et al. 2018; Benisty et al. 2018).

In addition, in Fig. 10, we showed that the radial spacing between rings is highly regular and depends solely on the aspect ratio of the disc, a quantity that can be constrained with observations of gas tracers at different altitudes from the disc mid-plane (e.g. Law et al. 2021). This regular spacing could be used as a diagnostic tool to infer whether the observed rings are due to a shadow or another mechanism, such as magnetic fields (Béthune, Lesur & Ferreira 2017; Riols, Lesur & Menard 2020) or undetected exoplanets (e.g. Zhang et al. 2018), as both of these mechanisms result in irregularly spaced rings in millimetre emission. In particular, discs such as around AS 209 (Guzmán et al. 2018) or SO 1274 (Huang et al. 2024) display many narrowly separated rings, which may be linked to the mechanism explored in this work.

Our results also suggest that the formation of substructure is encouraged for short cooling time-scales, which is generally true for hotter discs. With that in mind, it is possible that the formation of rings or vortices due to a shadow is more common in discs around Herbig Ae/Be stars.

Nevertheless, a misalignment of more than $\sim 10^\circ$ is required for the shadow to manifest as a sharp, narrow dark lane rather than an extended shadow. The latter would most likely induce warps at the disc surface rather than strong spirals and rings at the mid-plane, the dynamical consequences of which might be more difficult to observe.
## 6.5 Limitations of this work

In this study, we have made several simplifying assumptions. For one, we have worked in a vertically integrated (2D) framework, which neglects both the propagation of waves in the vertical direction and the three-dimensional structure of a shadow cast by a misaligned disc. As a result, the dynamical evolution of possible warps (e.g. Kimmig & Dullemond 2024) is not captured in our models. Our 2D framework further assumes that the dust is in settling-mixing equilibrium, which might not be a good approximation around the shadowed regions. A fully 3D hydrodynamical model would be required to address these issues, but would also be much more computationally expensive. This will be the focus of follow-up work.

Furthermore, in our models, we assume that the dust and gas grains are perfectly thermally coupled such that cooling is limited by the optical properties of dust grains rather than their collisional coupling to the gas. This assumption can easily break down in the outer disc, and especially so for larger grains (e.g. Dullemond et al. 2022; Muley, Melon Fuksman & Klahr 2023). In this case, the cooling time-scale would be longer than what is shown in Fig. 12 in the outer disc, possibly inhibiting the formation of rings for $R \gtrsim 70$ au. We nevertheless carry out a comparison between our radiation hydrodynamical models and a dynamical multifrequency radiative transfer model in Appendix B, showing that a large enough temperature contrast between the shadowed and illuminated regions is still present when accounting for the finite thermal coupling time-scale between gas and dust species.

Finally, we have neglected magnetic fields, and in particular the role of non-ideal MHD effects, which have been shown to result in the formation of rings (Riols et al. 2020). It is possible that this process could compound with the one studied in this paper, such that the MHD-driven gap opening would accelerate the formation of gaps seeded by the spirals launched at the shadow edges.

## 7 SUMMARY

In this study, we investigated the dynamical evolution of a protoplanetary disc with shadows cast by a misaligned inner disc using high-resolution 2D hydrodynamical simulations. Our models feature a self-consistent treatment of radiation transport, including radiative diffusion, and a dust component with dust-gas thermal and dynamical coupling.

We found that the thermal forcing due to the shadowing effect induces a series of spiral waves that deposit angular momentum into the disc at the shadow edges, leading to the formation of concentric gaps that trap millimetre grains efficiently, in turn forming bright rings in the dust surface density. The gap depth increases linearly with time, as spirals drive a sustained AMF into the disc, with the process being efficient as long as the local cooling time-scale is approximately shorter than the time it takes for gas to cross the shadow. Further analysis showed that the rings are spaced regularly in the radial direction, with their separation depending solely on the disc aspect ratio.

We then carried out a parameter study and found that the formation of rings in the disc is a robust process that is largely insensitive to the coagulation fraction $X$ and the viscous parameter $\alpha$ for reasonable levels of dust growth and turbulence, respectively, as long as the disc is sufficiently optically thin to allow for efficient cooling. The aspect ratio $h$ (as a proxy for the temperature), on the other hand, plays a crucial role in the formation of rings as the cooling time-scale depends sensitively on $h$. Models with a large $h$ (i.e. hotter discs) exhibit vigorous spirals and efficient gap opening, leading to Rossby

wave unstable gap edges and the formation of numerous vortices. The latter can appear as arcs, merge into clumps, or diffuse into rings depending on the local thermodynamics and dust-gas interaction.

Following up on our hydrodynamical models, we computed a suite of images of synthetic observations using the radiative transfer code RADMC-3D to check whether the features found in our models could be observable with current facilities such as ALMA and VLT-SPHERE. We found that both the rings and the shadow are visible in dust continuum emission at 1.3 mm for our fiducial model, and a dark lane can be seen in scattered-light images at 1.65 µm. The vortices found in models with higher aspect ratios are clearly visible as well, albeit with less pronounced shadows from the misaligned inner disc. Simulated ALMA images show that both the rings and vortices should be detectable with current facilities if present, with the characteristic regular spacing between rings being a diagnostic tool to infer the presence of a shadow.

Our results suggest that shadows can drive substructure in protoplanetary discs in the form of rings, arcs, and vortices, all of which should be observable with current facilities. Dust trapping at the resulting pressure bumps can facilitate the formation of planetesimals via the streaming instability, providing a possible solution to the long-standing chicken-and-egg problem of planet formation.

## ACKNOWLEDGEMENTS

The authors thank the anonymous referee for a careful and constructive report that helped identify errors related to the RADMC-3D models and greatly improved the quality of the paper. AZ would like to thank Xuening Bai, Ondřej Chrenko, Dhruv Muley, Giovanni Rosotti, and Prakrutı Sudarshan for helpful suggestions and fruitful conversations. This research utilized Queen Mary’s Apocrita HPC facility, supported by QMUL Research-IT (http://doi.org/10.5281/zenodo.438045). This work was performed using the DiRAC Data Intensive Service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. AZ and RPN were supported by STFC grants ST/T000341/1 and ST/X000931/1. TB acknowledges funding from the European Union under the European Union’s Horizon Europe Research and Innovation Programme 101124282 (EARLYBIRD) and funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 325594231, and Germany’s Excellence Strategy – EXC-2094 – 390783311. MB has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Framework Programme (PROTOPLANETS, grant agreement no. 101002188). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. All plots in this paper were made with the PYTHON library MATPLOTLIB (Hunter 2007). Typesetting was expedited with the use of GitHub Copilot, but without the use of any AI-generated text.

## DATA AVAILABILITY

Data from our numerical models are available upon reasonable request to the corresponding author.
REFERENCES

Alexiades V., Amiez G., Gremaud P.-A., 1996, Commun. Numer. Methods
Eng., 12, 31

Andrews S. M. et al., 2018, ApJ, 869, L41

Armitage P. J., 2009, Astrophysics of Planet Formation. Cambridge Univ.
Press, Cambridge

Avenhaus H. et al., 2017, AJ, 154, 33

Avenhaus H. et al., 2018, ApJ, 863, 44

Bai X.-N., Stone J. M., 2013, ApJ, 769, 76

Benisty M. et al., 2017, A&A, 597, A42

Benisty M. et al., 2018, A&A, 619, A171

Béthune W., Lesur G., Ferreira J., 2017, A&A, 600, A75

Béthune W., Latter H., Kley W., 2021, A&A, 650, A49

Birnstiel T., 2024, ARA&A, 62, 157

Bohn A. J. et al., 2019, A&A, 624, A87

Bohn A. J. et al., 2022, A&A, 658, A183

CASA Team, 2022, PASP, 134, 114501

Casassus S. et al., 2018, MNRAS, 477, 5104

Casassus S., Pérez S., Osses A., Marino S., 2019, MNRAS, 486, L58

Chiang E. I., Goldreich P., 1997, ApJ, 490, 368

Crida A., Morbidelli A., Masset F., 2006, Icarus, 181, 587

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

de Val-Borro M. et al., 2006, MNRAS, 370, 529

Doğan S., Nixon C., King A., Price D. J., 2015, MNRAS, 449, 1251

Dominik C., Min M., Tazaki R., 2021, Astrophysics Source Code Library,
record ascl:2104.010

Dullemond C. P., van Zadelhoff G. J., Natta A., 2002, A&A, 389, 464

Dullemond C. P., Juhasz A., Pohl A., Sereshti F., Shetty R., Peters T.,
Commercon B., Flock M., 2012, Astrophysics Source Code Library,
record ascl:1202.015

Dullemond C. P., Ziampras A., Ostertag D., Dominik C., 2022, A&A, 668, A105

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

Flock M., Turner N. J., Nelson R. P., Lyra W., Manger N., Klahr H., 2020,
ApJ, 897, 155

Fromang S., Nelson R. P., 2009, A&A, 496, 597

Fung J., Ono T., 2021, ApJ, 922, 13

Gammie C. F., 2001, ApJ, 553, 174

Guilera O. M., Sándor Z., Ronco M. P., Venturini J., Miller Bertolami M. M.,
2020, A&A, 642, A140

Guzmán V. V. et al., 2018, ApJ, 869, L48

Henyey L. G., Greenstein J. L., 1941, ApJ, 93, 70

Huang J. et al., 2018, ApJ, 869, L42

Huang J., Ansdell M., Birnstiel T., Czekala I., Long F., Williams J., Zhang
S., Zhu Z., 2024, ApJ, 976, 132

Hubeny I., 1990, ApJ, 351, 632

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Jiang H., Ormel C. W., 2023, MNRAS, 518, 3877

Juhász A., Facchini S., 2017, MNRAS, 466, 4053

Keppler M. et al., 2018, A&A, 617, A44

Kimmig C. N., Dullemond C. P., 2024, A&A, 689, A45

Kley W., 1989, A&A, 208, 98

Kuffmeier M., Dullemond C. P., Reissl S., Goicovic F. G., 2021, A&A, 656, A161

Kurtovic N., 2024, J. Open Source Softw., 9, 4942

Lau T. C. H., Drażkowska J., Stammler S. M., Birnstiel T., Dullemond C. P.,
2022, A&A, 668, A170

Lau T. C. H., Birnstiel T., Drażkowska J., Stammler S. M., 2024, A&A, 688, A22

Law C. J. et al., 2021, ApJS, 257, 4

Levermore C. D., Pomraning G. C., 1981, ApJ, 248, 321

Long F. et al., 2018, ApJ, 869, 17

Lovascio F., Paardekooper S.-J., McNally C., 2022, MNRAS, 516, 1635

Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805

Marino S., Perez S., Casassus S., 2015, ApJ, 798, L44

Masset F., 2000, A&AS, 141, 165

Menou K., Goodman J., 2004, ApJ, 606, 520

Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C.,
Ferrari A., 2007, ApJS, 170, 228

Mignone A., Flock M., Stute M., Kolb S. M., Muscianisi G., 2012, A&A,
545, A152

Min M., Hovenier J. W., de Koter A., 2005, A&A, 432, 909

Miranda R., Rafikov R. R., 2020a, ApJ, 892, 65

Miranda R., Rafikov R. R., 2020b, ApJ, 904, 121

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

Montesinos M., Perez S., Casassus S., Marino S., Cuadra J., Christiaens V.,
2016, ApJ, 823, L8

Morbidelli A., 2020, A&A, 638, A1

Morfill G. E., Voelk H. J., 1984, ApJ, 287, 371

Muley D., Melon Fuksman J. D., Klahr H., 2023, A&A, 678, A162

Muro-Arena G. A. et al., 2020, A&A, 635, A121

Nealon R., Dipierro G., Alexander R., Martin R. G., Nixon C., 2018, MNRAS,
481, 20

Nealon R., Price D. J., Pinte C., 2020, MNRAS, 493, L143

Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610

Papaloizou J. C. B., Pringle J. E., 1983, MNRAS, 202, 1181

Pinilla P. et al., 2018, ApJ, 868, 85

Pinte C. et al., 2018, ApJ, 860, L13

Raettig N., Klahr H., Lyra W., 2015, ApJ, 804, 35

Rafikov R. R., 2002, ApJ, 569, 997

Riols A., Lesur G., Menard F., 2020, A&A, 639, A95

Rometsch T., Ziampras A., Kley W., Béthune W., 2021, A&A, 656, A130

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

Stolker T. et al., 2017, ApJ, 849, 143

Su Z., Bai X.-N., 2024, ApJ, 975, 126 (SB24)

Tassoul J.-L., 1978, Theory of Rotating Stars. Princeton Univ. Press, Prince-
ton, NJ

Teague R., Bae J., Bergin E. A., Birnstiel T., Foreman-Mackey D., 2018, ApJ,
860, L12

Teague R., Bae J., Bergin E. A., 2019, Nature, 574, 378

Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25

Van Leer B., 1974, J. Comput. Phys., 14, 361

Weber P., Pérez S., Benítez-Llambay P., Gressel O., Casassus S., Krapp L.,
2019, ApJ, 884, 178

Youdin A. N., Lithwick Y., 2007, Icarus, 192, 588

Young A. K., Alexander R., Walsh C., Nealon R., Booth A., Pinte C., 2021,
MNRAS, 505, 4821

Zhang S., Zhu Z., 2020, MNRAS, 493, 2287

Zhang S., Zhu Z., 2024, ApJL, 974, L38

Zhang S. et al., 2018, ApJ, 869, L47

Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ, 813, 88

Zhu Z., Zhang S., Johnson T. M., 2025, ApJ, 980, 259

Ziampras A., Ataiee S., Kley W., Dullemond C. P., Baruteau C., 2020, A&A,
633, A29

Ziampras A., Nelson R. P., Rafikov R. R., 2023, MNRAS, 524, 3930

Ziampras A., Sudarshan P., Dullemond C. P., Flock M., Berta V., Nelson R.
P., Mignone A., 2025, MNRAS, 536, 3322

Zubko V. G., Mennella V., Colangeli L., Bussoletti E., 1996, MNRAS, 282, 1321

APPENDIX A: TWO-DISC MODEL

To set up our two-disc model, we first define the cylindrical
radius $R = \sqrt{x^2 + y^2}$ with respect to a Cartesian coordinate system
$\{x, y, z\}$. The mid-plane is then defined by $z = 0$, and the disc is
assumed to be vertically isothermal with an aspect ratio $h(R) =
H(R)/R$ and a radial surface density profile $\Sigma(R)$, tapered between
$R_{\text{in}}$ and $R_{\text{out}}$ as

$$
\Sigma(R) = \Sigma_0 \left(\frac{R}{R_{\text{in}}}\right)^{-1} \frac{1}{1 + \exp\left(\frac{R_{\text{in}}-R}{0.01 R_{\text{in}}}\right)} \frac{1}{1 + \exp\left(\frac{R-R_{\text{out}}}{0.01 R_{\text{out}}}\right)}. \quad \text{(A1)}
$$
The gas volume density is then given by (e.g. Nelson et al. 2013)

$$ \rho_g(R, z) = \rho_g^{\text{mid}} e^{-\frac{1}{h^2}\left(1-\frac{R}{\sqrt{R^2+z^2}}\right)}, \quad \rho_g^{\text{mid}}(R) = \frac{1}{\sqrt{2\pi}} \frac{\Sigma}{H}. \quad (A2) $$

The disc can then be arbitrarily rotated about the y-axis by an angle $\theta_{\text{rot}}$ with respect to the x-y plane by applying the transformation

$$ x' = x \cos \theta_{\text{rot}} + z \sin \theta_{\text{rot}}, \quad y' = y, \quad z' = -x \sin \theta_{\text{rot}} + z \cos \theta_{\text{rot}}, \quad (A3) $$

and finally $R' = \sqrt{x'^2 + y'^2}$.

We first define a fixed grid in spherical coordinates $\{r, \theta, \phi\}$, through which we compute $\{x, y, z\}$ at all points. We then rotate the Cartesian coordinate system by $\theta_{\text{rot}}$ to obtain $\{x', y', z'\}$, define $R'$, and finally compute the gas density $\rho_g(R', z')$ through the above equations. We set $R_{\text{in}} = 0.2$ au, $R_{\text{out}} = 4$ au, and $\theta_{\text{rot}} = 30^{\circ}$ for the inner disc, and $R_{\text{in}} = 5$ au, $R_{\text{out}} = 300$ au, and $\theta_{\text{rot}} = 0^{\circ}$ for the main disc. The total gas density is then given by the sum of the contributions of both discs. Finally, by assuming that the small grains are perfectly coupled to the gas, we write the dust density $\rho_d^{\text{small}}$ as

$$ \rho_d^{\text{small}} = \varepsilon(1-X)\rho_g, \quad (A4) $$

where $X = 0.9$ is the coagulation fraction defined in equation (4), and $\varepsilon = 0.01$ is the dust-to-gas ratio. This approach allows us to model the two-disc system on a single, spherical, regular grid, on which ray tracing from the central star is straightforward.

## APPENDIX B: COMPARISON WITH DYNAMIC RADIATIVE TRANSFER

To verify how accurate the time-dependent heating/cooling prescription of this paper is, we carry out a 1D vertical time-dependent radiative transfer calculation, for a column of gas and dust in approximate vertical hydrostatic equilibrium orbiting around the star and periodically passing through a shadow. In this appendix, we describe the model, albeit in abbreviated form, and refer to Dullemond, van Zadelhoff & Natta (2002) and Dullemond et al. (2022) for details.

The model is set up on a vertical grid of 90 grid points in the coordinate $z$, with $z = 0$ being the mid-plane of the disc. The grid spans between $z = 0$ and $z = 0.4 R$, with $R$ being the cylindrical radial distance to the central star which we set to $R = 60$ au. For simplicity, we assume mirror symmetry in the mid-plane. We are aware of the fact that for disc misalignments unequal to $90^{\circ}$ (or $0^{\circ}$), the shadows of the inner disc onto the outer disc will not be mirror symmetric. However, the dynamic models of this paper do not account for this top/bottom asymmetry, so we keep the 1D model symmetric.

We first compute a steady-state model. The methods used are very similar to those described by Dullemond et al. (2002). The only exception is the convergence algorithm with which the radiative transfer is computed: we use the moment method with fixed mean opacities only for obtaining a first guess. Then we use the simple Lambda Iteration method with Ng acceleration to obtain the accurate results. The reason why we can afford the Lambda Iteration scheme is that we focus on the outer disc regions, which have much lower optical depth than the very inner disc regions for which the moment method was developed.

The parameters of our model are $M_{\star} = 1 M_{\odot}$, $L_{\star} = 1 L_{\odot}$ with a solar-type spectrum, $\Sigma_g = 16.67 \text{ g cm}^{-2}$, a dust growth parameter of $X = 0.9$, standard DIANA opacities from OPTOOL with material density $2.08 \text{ g cm}^{-3}$ for small grains with radius $0.1 \text{ µm}$ and pebbles with radius $1 \text{ mm}$. We use an incidence irradiation angle of

<figure><img src="image_21.png" /><figcaption>**Figure B1.** Mid-plane gas temperature evolution using the dynamic radiative transfer model described in this section (black) and models following the set-up in Section 2. The two models agree well in terms of the peak-to-peak temperature variation. The difference in phase owes both to the partial thermal dust-gas decoupling and the lack of radiative diffusion in the 1D model.</figcaption></figure>

0.0361 rad, yielding an initial approximate mid-plane temperature of 18.2 K, which yields an aspect ratio of $h \simeq 0.066$. The vertical distribution of the gas is then taken to be Gaussian with that pressure scale height. We assume that both the small grains and the pebbles are vertically well mixed. For the wavelength grid, we use 100 gridpoints logarithmically spaced between $0.1 \text{ µm}$ and $2 \text{ mm}$. At each wavelength, the radiation field has an angular grid of 20 points in $\mu = \cos \theta$ in the range $-1 \le \mu \le 1$. After the steady-state radiative transfer has converged, the mid-plane temperature is 16.72 K for both the gas and the small grains, and 15.9 K for the pebbles.

Next, we use this steady-state model as the initial condition for a time-dependent computation. As we comove along the orbit, our 1D column passes through the shadow twice per orbit. The shadow depth is $10^{-6}$, the shadow width is $\Delta\phi = 0.48$ rad, corresponding to a pass-through time of 35.5 yr. In the time-dependent model, the radiative transfer is assumed to be instant, as the light-crossing time for one scale height is only 2250 s (i.e. a factor of half a million times shorter than the crossing time of the gas through the shadow). However, the gas, dust, and pebbles are not in thermal equilibrium with the radiation field. The time dependence of the model is thus the gain/loss of heat from these three components. The physics is described in this paper and in the appendices of Dullemond et al. (2022). We include the dust-gas coupling, and we assume that the gas, having no continuum opacity, cannot efficiently couple to the radiation field. So the radiative heating and cooling go solely through the small grains and the pebbles. At each location, the radiative heating rate of the small dust grains (d) or pebbles (p) is

$$ q_{d/p}^{+}(z) = 4\pi\rho_{d/p}(z) \int_{0}^{\infty} \kappa_{v,d/p}^{\text{abs}} J_{v}(z) \, dv, \quad (B1) $$

in units of energy per time per volume, where $J_v(z)$ is the angular mean intensity of the radiation field:

$$ J_v(z) = \frac{1}{2} \int_{-1}^{+1} I_{\mu,\nu}(z) \, d\mu, \quad (B2) $$
where $I_{\mu,\nu}(z)$ is the radiative intensity as a function of $z$, computed from the radiative transfer part of the algorithm. The radiative cooling rate is

$$q_{d/p}(z) = 4\pi\rho_{d/p}(z) \int_{0}^{\infty} \kappa_{\nu,d/p}^{\text{abs}} B_{\nu}(T_{d/p}(z)) \, d\nu, \quad (B3)$$

where $B_{\nu}$ is the Planck function and $T_{d/p}(z)$ is the temperature of the dust or pebble at a given time. In addition, we have the dust-gas (or pebble-gas) thermal coupling $q_{d/p,g}(z)$ (Dullemond et al. 2022). When $q_{d/p}^{+}(z)$, $q_{d/p}^{-}(z)$, and $q_{d/p,g}(z)$ are not in equilibrium, the temperatures of the three components (dust, pebbles, and gas) change, as described in Dullemond et al. (2022).

We can now compare the resulting mid-plane temperatures as a function of time (or phase) along the orbit, and compare to those computed using the approximate method of this paper. For the latter, we run two dedicated simulations with $L_{*} = 1.3 L_{\odot}$ in order to match the temperature of 16.72 K at $R = 60$ au, with the first following the exact set-up in Section 2 and the second assuming that big grains are thermally decoupled (i.e. $\kappa_{R,P}^{\text{big}} = 0$). The result is shown in Fig. B1, where it can be seen that the two methods match reasonably well. The temperature of the pebbles is always a bit lower than that of the gas and the small grains, but that is to be expected because the opacity of the big grains is flatter than that of the small grains and thus tends to lead to lower temperatures.

## APPENDIX C: AZIMUTHAL T VARIATIONS IN RADMC-3D

As mentioned in Section 5.1, RADMC-3D is a hydrostatic code and therefore does not capture the azimuthal advection of thermal energy. This results in artificially narrower temperature fluctuations compared to our hydrodynamical simulations. To highlight this effect we show the azimuthal profiles of $T_{\text{small}} \approx T_{\text{gas}}$ from the RADMC-3D output of our fiducial model at $t = 1000 P_{0}$ and $R = 40$ au in Fig. C1 and compare it to the PLUTO output at the same timestamp and radius.

As seen in that figure, deviations from said background are found only around the shadowed region in RADMC-3D, whereas they span the full azimuthal range in PLUTO. The RADMC-3D model also shows

weaker deviations compared to the background near the mid-plane compared to PLUTO. Nevertheless, the two codes agree very well at the mid-plane in terms of the background temperature.

<figure><img src="image_22.png" /><figcaption>Figure C1. Azimuthal slice at $R = 40$ au of the gas and small dust temperatures, normalized to the azimuthal median for RADMC-3D (blue and orange) and to the initial temperature at that radius for PLUTO (black). The latter shows a wider spread of temperatures due to the azimuthal advection of thermal energy. RADMC-3D and PLUTO show excellent agreement in background temperature at the mid-plane.</figcaption></figure>

This paper has been typeset from a TeX/Latex file prepared by the author.