# Dynamical Consequence of Shadows Cast to the Outer Protoplanetary Disks: I.
Two-dimensional Simulations

Zehao Su School of Physics and Astronomy, Beijing Normal University, Beijing
100875, China<EMAIL_ADDRESS>Institute for Frontier in Astronomy and
Astrophysics, Beijing Normal University, Beijing 102206, China;  Xue-Ning Bai
Institute for Advanced Study, Tsinghua University, Beijing 100084, China;
<EMAIL_ADDRESS>Department of Astronomy, Tsinghua University, Beijing
100084, China

###### Abstract

There has been increasing evidence of shadows from scattered light
observations of outer protoplanetary disks (PPDs) cast from the (unresolved)
disk inner region, while in the meantime these disks present substructures of
various kinds in the submillimeter. As stellar irradiation is the primary
heating source for the outer PPDs, the presence of such shadows thus suggest
inhomogeneous heating of the outer disk in azimuth, leading to a “thermal
forcing” with dynamical consequences. We conduct a suite of idealized 2D disk
simulations of the outer disk with azimuthally-varying cooling prescription to
mimic the effect of shadows, generally assuming the shadow is static or
slowly-rotating. The linear response to such shadows is two-armed spirals with
the same pattern speed as the shadow. Towards the nonlinear regime, we find
that shadows can potentially lead to the formation of a variety of types of
substructures including rings, spirals and crescents, depending on viscosity,
cooling time, etc. We have conducted systematic and statistical
characterization of the simulation suite, and as thermal forcing from the
shadow strengthens, the dominant form of shadow-induced disk substructures
change from spirals to rings, and eventually to crescents/vortices. Our
results highlight the importance of properly modeling the dynamical impact of
inhomogeneous stellar irradiation, while call for more detailed modeling
incorporating more realistic disk physics.

Protoplanetary disks-shadows-hydrodynamics substructures

##  1 Introduction

Thanks to the advent of new observational facilities and instruments to
conduct spatially resolved observations of protoplanetary disks (PPDs), it has
now been well established that disk substructures are ubiquitous (e.g. Bae et
al., 2023). In the millimeter/sub-millimeter wavelengths, the Atacama Large
Millimeter Array (ALMA) has revealed the richness of disk substructures that
are primarily in the form of rings and gaps in addition to other features such
as spirals and crescents at different radii (e.g. ALMA Partnership et al.,
2015; Monnier et al., 2017; Avenhaus et al., 2018; Isella et al., 2018; Huang
et al., 2020; Gratton et al., 2019; Andrews, 2020). These observations reflect
the thermal emission from mm-sized dusts around the disk midplane, which are
biased tracers of the gas density profiles due to finite aerodynamic coupling
between gas and dust. In the optical and near infrared (NIR), high-contrast
imaging with extreme adaptive optics (e.g., VLT/SHERE, VLT/CRIRES, GPI) reveal
even richer and more complex features (e.g. Benisty et al., 2015; Pinilla et
al., 2015; Pohl et al., 2017; van Boekel et al., 2017; Garufi et al., 2018;
Benisty et al., 2023). Emission at optical/NIR mainly result from the
starlight scattered by micron-sized dust (better coupled to the gas) suspended
in the disk, which are better tracers of the disk surface layers. As a result,
features seen in scattered light do not necessarily have direct correspondence
to substructures recognized by ALMA (e.g. van der Marel et al., 2016; Uyama et
al., 2018; Pérez et al., 2018; Muro-Arena et al., 2018).

At least partly contributing to the complexity in features seen in scattered
light is the presence of shadows, typically defined as low-intensity regions
that are confined to specific azimuthal angles (Benisty et al., 2023). They
must be cast from the (unresolved) disk inner region, and can be mainly
classified into two types: broad extended shadows in azimuthal directions
(e.g. Muro-Arena et al., 2020) and narrow shadow lanes with only a few degrees
(e.g. Ginski et al., 2021). Considerable effort has been devoted to
understanding the origin of shadows because the morphology and temporal
variation of shadows can provide indirect information about the disk’s inner
regions. The most common case for shadow casting is the presence of a
misaligned/warped inner disk. For instance, TW Hya shows a moving shadow
pattern that could suggests a precessing inner disk (Debes et al., 2017),
shadows in HD 143006 can be reproduced using a $30^{\circ}$ misaligned inner
disk (Benisty et al., 2018), fast time variations of shadows in RX
J1604.3-2130A may come from dust very close to the central star in an
irregular misaligned inner disk (Pinilla et al., 2015), narrow shadow lanes in
SU Aur that possibly suggest misalignment caused by late-time interactions
with infalling materials (Kuffmeier et al., 2021), shadows in HD 139614 can be
explained by the combination of a misaligned inner ring and disk (Muro-Arena
et al., 2020), and the flux ratio switches sides in brightest nebula of
IRAS40302 can be achieved by applying a tilted inner disk (Villenave et al.,
2023). Even for disks with nearly aligned inner regions, subtle shadowing
effects can still be recognized (Monnier et al., 2017). In addition to the
misalignment of the inner disk regions, variations in the scale height of the
inner disk atmosphere could also be responsible for generating shadows, such
as in HD 163296 (Rich et al., 2019, 2020; Varga et al., 2021; GRAVITY
Collaboration et al., 2021).

Most effort aiming to understand disk shadows so far has focused on the
modeling the (inner) disk morphology to explain the shadow features using
radiative transfer calculations (e.g. Casassus et al., 2018; Benisty et al.,
2018; Nealon et al., 2019; Muro-Arena et al., 2020). On the other hand, we
note that as stellar irradiation is the primary source of heating in the bulk
of the (outer) PPDs, the presence of shadows must also give rise to dynamical
consequences by “thermal forcing”: the disk gas experiences (quasi-)periodic
cooling and heating as it enters and exits the shadow, which hardly settles to
thermal equilibrium and constantly exerts modest or even strong pressure
perturbations to the neighboring fluid. This effect was first explored in
Montesinos et al. (2016), who conducted 2D hydrodynamic simulations that take
into account both stellar irradiation and periodic forcing of shadows with an
opening angle of $28^{\circ}$ in the context of the transition disk HD 142527.
They found that azimuthal pressure gradients generated by shadows can trigger
$m=2$ spirals, which are enhanced by self-gravity and give rise to observable
quasi-steady spiral signals.

In this work, motivated by the diversity of shadowing features seen in
scattered light images and the case study by Montesinos et al. (2016) for the
HD 142527 disk, we aim at a systematic exploration on the dynamical
consequences of shadows cast onto outer PPDs. As an initial effort, we
restrict ourselves to vertically-integrated systems by 2D hydrodynamic
simulations. We follow the evolution of a passive, viscous gaseous disk with a
thermal relaxation prescription towards a target temperature, which is set by
stellar irradiation subject to shadowing. By exploring a large suite of
simulations varying the viscosity, cooling time, and shadow geometry, we find
that shadows can result in the formation of a wide variety of disk
substructures, and perform a statistical analysis of all the substructures
generated from our simulations.

This paper is organized as follows. We detail our simulation setup in Section
2, followed by a description representative features of shadow-driven disk
substructures in Section 3. We present statistical analysis of all our
simulations and describe the substructure-forming process from linear to non-
linear regimes in Section 4. Finally, we discuss the caveats and conclude in
Section 6.

##  2 Numerical Methods

###  2.1 Simulation Setup

We solve the vertically integrated viscous hydrodynamic equations using the
grid-based higher-order Godunov code ATHENA++ Stone et al. (2020) in
cylindrical coordinates $(r,\varphi)$. The conservative form of control
equations are:

$\frac{\partial\Sigma}{\partial t}+\nabla\cdot(\Sigma\textbf{v})=0,$ (1)

$\frac{\partial\Sigma\textbf{v}}{\partial
t}+\nabla\cdot(\Sigma\textbf{v}\textbf{v}+P\mathcal{I}+\mathcal{T}_{vis})=-\Sigma\nabla\Phi,$
(2)

$\frac{\partial E}{\partial
t}+\nabla\cdot\left[(E+P+\mathcal{T}_{vis})\textbf{v}\right]=-\Sigma\textbf{v}\cdot\nabla\Phi+\Lambda,$
(3)

where $\Sigma$ is the gas surface density in the disk, v is the gas velocity,
$P=\Sigma T$ is the vertically integrated pressure where $T$ is the disk
temperature, $\mathcal{I}$ is the identity tensor, $\Phi$ is the gravitational
potential written as $\Phi=-GM_{*}/r$ where $M_{*}$ is mass of central star,
$E$ is total energy density, and $\Lambda$ gives the cooling source terms. The
viscous stress tensor $\mathcal{T}_{visc}$ in momentum equation reads

$\mathcal{T}_{vis}=-\Sigma\nu\left(\frac{\partial v_{i}}{\partial
x_{j}}+\frac{\partial v_{j}}{\partial x_{i}}-\frac{2}{3}\frac{\partial
v_{k}}{\partial x_{k}}\delta_{ij}\right),$ (4)

with $\nu$ being the kinematic viscosity.

The total energy density $E$ is given by the combination of kinetic energy and
internal energy:

$E=\frac{1}{2}\Sigma v^{2}+\frac{P}{\gamma-1},$ (5)

where $\gamma=7/5$ is the adiabatic index for molecular gas. Note that viscous
heating is automatically included in the energy equation, although it is
generally unimportant in the outer PPDs. The gas temperature $T$ is associated
with the isothermal sound speed as $T=c_{s}^{2}$, which yields the disk scale
height $H=c_{s}/\Omega_{K}$, where $\Omega_{K}=(GM_{*}/r^{3})^{1/2}$ is the
Keplerian angular frequency. The disk aspect ratio is then given by $h=H/r$.
With this, viscosity follows the standard $\alpha$ prescription (Shakura &
Sunyaev, 1973), $\nu=\alpha c_{s}H$. It is worth noticing that viscosity
varies as the disk evolves.

Figure 1: Illustration of our shadow morphology depicted by normalized
temperature in disks with $\sigma_{\phi}=0.236$ and $0.079$ for fast
($\beta=0.001$) and slow ($\beta=1$) cooling processes. With fast cooling, the
low-temperature regions caused by shadows largely reflect the target
temperature we prescribe at the shadow location. With longer cooling
timescale, these regions deviate from the prescribed shadow positions towards
the leading side, and exhibit significantly weaker amplitude.

We choose the initial density profile to be power-law:

$\Sigma_{\rm init}=\Sigma_{0}\left(\frac{r}{r_{0}}\right)^{d}.$ (6)

where $\Sigma_{0}$ is the density at reference radius $r_{0}$, which we take
to be the radius of the inner boundary. We specify the disk initial
temperature as

$T_{\rm init}\equiv
c_{s_{0}}^{2}\left(\frac{r}{r_{0}}\right)^{p}=\left(h_{0}r_{0}\Omega_{K_{0}}\right)^{2}\left(\frac{r}{r_{0}}\right)^{p},$
(7)

where $c_{s_{0}}$, $h_{0}$, $\Omega_{K_{0}}$ are the value of isothermal sound
speed, aspect ratio and Keplerian angular velocity at reference radius
$r_{0}$.

The radial force balance leads to the initial rotation profile calculated by

$v_{\varphi}(r)=\left[\left(p+d\right)c_{s}^{2}+\frac{GM_{*}}{r}\right]^{1/2}\
.$ (8)

With viscosity, the radial velocity is set by the accretion velocity given by

$v_{r}(r)=-\frac{3}{2}\frac{\alpha c_{s}^{2}}{r\Omega_{K}}\ .$ (9)

To ensure steady-state accretion in the initial equilibrium (shadow not
included) with constant $\alpha$, the initial temperature and density profile
should satisfy $d+p=-3/2$.

Our simulations are scale-free, adopting $GM_{*}=\Sigma_{0}=r_{0}=1$ in code
units, with $h_{0}=0.1$. As a result, we have $\Omega_{K_{0}}=c_{s_{0}}=1$.
The computational domain ranges from $r_{\rm in}=1$ to $r_{\rm out}=30$ in
code units, to ensure sufficient dynamical range. We employ a logarithmic grid
in radial direction and a uniform grid in azimuthal direction with
$N_{r}\times N_{\varphi}=512\times 512$, achieving a gird resolution of 15
cells per $H$ in $r$ while keeping cell aspect ratio $\Delta r\approx
0.5r\Delta\phi$.

####  2.1.1 Shadow prescription

In our simulations, we assume an obscuring structure present in disk inner
region, which is outside of our simulation domain (inside the inner boundary).
As an initial study, we tentatively take this obscuring structure as a
slightly misaligned inner disk (inclination angle $\sim h$). In this case,
near half of the outer disk (in azimuth) is illuminated from one hemisphere,
and the opposite side being illuminated from the other hemisphere. In the
vertically-integrated sense, the two sides of the “pie-chart” are heated
largely equally. It is the transition region, the disk can be largely blocked
by the inner disk from both hemisphere, which is mostly affected by the
shadow. The shadow then introduces a thermal forcing to the system, causing
system’s temperature to approach the target temperature $T_{\rm tar}$. For
simplicity, we prescribe this target temperature by

$T_{\rm tar}(r,\phi)=T_{\rm init}(r)\left(1-\epsilon
e^{-\frac{\phi^{2}}{2\sigma_{\phi}^{2}}}\right)\left(1-\epsilon
e^{-\frac{(\phi-\pi)^{2}}{2\sigma_{\phi}^{2}}}\right),$ (10)

where $\epsilon$ reflects the amplitude of the shadow and $\sigma_{\phi}$
characterizes the azimuthal width of the shadow. Although we have argued that
two-sided shadows are the most basic case, we still examine the one-sided case
in C for reference. Also, in most cases in this paper, except for the
simulation mentioned in Section 5.3, the shadow is static with pattern speed
($\Omega_{\rm shadow}$) being zero.

The final temperature structure depends on heating and cooling process, which
is often modeled using the $\beta$ cooling approximation (Gammie, 2001). The
cooling term is given by thermal relaxation towards the target temperature

$\Lambda=-\frac{\Sigma}{\gamma-1}\times\frac{T-T_{\rm tar}}{t_{\rm cool}}\ ,$
(11)

where the cooling timescale is specified by the dimensionless parameter
$\beta$:

$t_{\rm cool}=\beta\Omega_{K}^{-1}\ .$ (12)

It describes the disk’s thermodynamic timescale, which can range from $\sim
10^{-3}$ (approaching the isothermal limit) to at least $\sim 10$ (approaching
the adiabatic limit) in our simulations. Figure 1 shows the expected
temperature structure for four representative shadow prescriptions, with
different shadow widths ($15^{\circ}$ and $45^{\circ}$) and cooling times,
calculated by following fluid elements undergoing heating and cooling on
circular orbits. We have fixed shadow amplitude of $\epsilon=0.8$. With fast
cooling ($\beta=0.001$), we see that the shadow aligns with its expected
position, and the temperature at its center is approximately $0.2T_{0}$ as
desired. When cooling is inefficient, the observed shadow center on the
leading side from its expected location, and the lowest temperature becomes
well above $0.2T_{0}$.

It should be noted that our shadow and cooling prescriptions are highly
simplified and are not necessarily always physical (for instance, a flat disk
with $p=-1$ would not be irradiated). We emphasize that the goal of this work
is not to precisely model any particular system, but to explore the general
phenomenology in a qualitative manner.

####  2.1.2 Boundary Conditions

We use modified outflow boundary conditions, where hydrodynamic variables are
copied from the last grid zone assuming $\Sigma\propto r^{d}$, $P\propto
r^{d+p}$, $v_{\phi}\propto r^{-1/2}$, with $v_{r}$ unchanged—except we set
$v_{r}=0$ in case of inflow. To further dampen unphysical waves, we adopt
wave-killing functions in the form described by de Val-Borro et al. (2006):

$\frac{dx}{dt}=-\frac{x-x_{0}}{\tau_{\rm damp}}R(r)$ (13)

where $x$ represents any fluid quantities (e.g. $\Sigma$, v, etc.). The
damping timescale $\tau_{\rm damp}$ is defined as $\tau_{\rm
damp}=\eta\Omega_{K}^{-1}$, where $\eta$ is the damping rate and is set to 1
for all of our simulations. The function $R(r)$ is a parabolic function
expressed as:

$R(r)=\left(\frac{r-r_{\rm damp}}{L_{\rm damp}}\right)^{2},\ {\rm
for}\left|r-r_{\rm damp}\right|<L_{\rm damp}\ ,$ (14)

where $r_{\rm damp}$ is the boundary of the damping area which we take to be
$2.08$ and $26.57$ in inner and outer part of our computation domain,
respectively, and $L_{\rm damp}$ is the length of wave killing zone.

###  2.2 Simulation Runs

In order to comprehensively investigate the dynamical effects of shadows in
PPDs, we conducted a wide parameter scan. Five main parameters are included in
our simulations: dimensionless cooling timescale $\beta$ ranging from
$10^{-3}$ to $10$, viscosity coefficient $\alpha$ ranging from $0$ to
$10^{-2}$, shadow amplitude coefficient $\epsilon=0.5$ and $\epsilon=0.8$,
shadow width $\sigma_{\phi}=0.236$ and $\sigma_{\phi}=0.079$, and the
temperature slope $p=-1$ (flat case) and $p=-0.5$ (flared case). In the
viscous simulations, they translate to density gradient $d=-0.5$ and $d=-1$ to
ensure steady state accretion. In most simulations, the shadows do not rotate,
and we fix the disk aspect ratio $h_{0}=0.1$ at $r=1$, thus $h=0.1$ is
constant in most $p=-1$ (flat) cases. All of these simulations will be
discussed in Section 4. In Sections 5.3 and 5.4, we also briefly explore
simulations with rotating shadows and vary $h_{0}$ from $0.05$ to $0.15$.

To further comment on our choice of parameters, we first note that in outer
disk conditions, we generally expect $\beta\lesssim 1$ (e.g. Lin & Youdin,
2015; Pfeil & Klahr, 2019), though the finite thermal coupling between dust
and gas may significantly enhance the effective $\beta$ (e.g. Bae et al.,
2021). In inviscid simulations, we further examine the influence of the
density profile ($d=-0.5$ and $d=-1$, which affects thermal forcing) since
this parameter is no longer free when viscosity is included (dependent on
$p$). Note that in the new paradigm of wind-driven accretion, the disk is more
laminar and the surface density profile can be more arbitrary (e.g. Bai, 2016;
Suzuki et al., 2016; Tabone et al., 2022). Although we do not incorporate
wind-driven accretion, this exploration also serves the purpose to partly
mimic “windy” disk conditions. On the choice of shadow amplitudes, note that
given the $T^{4}$ dependence, the two choices correspond to the shadowed
region receiving about $0.2^{4}\approx 0.002$ and $0.5^{4}\approx 0.06$ of the
stellar irradiation compared to the non-shadowed regions.

In all of our simulations, the total run time is chosen to be $T=20000P_{0}$,
where $P_{0}=2\pi/\Omega_{K_{0}}$ is the orbital period at the inner boundary.
This is significantly longer than the timescales for substructure formation,
which we find to be within $5000P_{0}$ for most cases. In only a few cases
(especially $\beta\sim 10,\alpha\sim 0$), even on the timescale of
$20000P_{0}$, we cannot unambiguously identify the dominant form of disk
substructure. However, we can infer their evolution trend from a statistical
point of view.

To facilitate comparison of the various simulations discussed in the following
sections, we provide a list of all our runs and their parameters in Table 1.
Our naming convention is structured as follows. We use “L” for runs in the
linear regime ($\epsilon=0.001$) and “NL” for runs in the nonlinear regime
($\epsilon=0.5,0.8$). The labels “hs”, “hm” and “hl” indicate runs with
$h_{0}=0.05$, $h_{0}=0.1$, and $h_{0}=0.15$, respectively, with $h_{0}=0.1$ as
fiducial. To specify the dominant substructure, we use “S” for spiral-
dominant, “R” for ring-dominant, and “V” for vortex-dominant. Shadow
precession speeds are denoted as “NR” for nonrotating (fiducial), “FR” for
fast rotating, “MR” for moderately rotating, and “SR” for slow rotating. For
simulations dedicated to parameter searches discussed in Section 4, we use the
label “S-h-all,” as we do not discuss individual runs for these simulations.

Table 1: Summary of All Highlighted Simulations.1

Run | $\sigma_{\phi}$ | $\epsilon$ | $\alpha$ | $\beta$ | p | $h_{0}$ | $\Omega_{\rm shadow}$ |
---|---|---|---|---|---|---|---|---
Representative runs (Section 3) |
NL-hm-S-NR | 0.236 | 0.8 | $10^{-3}$ | 10 | -1.0 | 0.1 | 0 |
NL-hm-R-NR | 0.236 | 0.5 | $10^{-4}$ | 1 | -1.0 | 0.1 | 0 |
NL-hm-V-NR | 0.236 | 0.5 | $0$ | $10^{-3}$ | -1.0 | 0.1 | 0 |
Statistical runs (Section 4) |
S-h-all2 | $(0.236,0.079)$ | $(0.5,0.8)$ | $(0,10^{-4},10^{-3},10^{-2})$ | $(10^{-3},10^{-2},10^{-1},1,10)$ | $(-1.0,-0.5)$ | 0.1 | 0 |
Linear run (Section 5.1) |
L-hm-S-NR | 0.236 | 0.001 | 0 | $10^{-3}$ | -1.0 | 0.1 | 0 |
Rotating shadow runs (Section 5.3) |
L-hm-S-FR3 | 0.236 | 0.001 | 0 | $10^{-3}$ | -1.0 | 0.1 | $\Omega_{0}$ |
L-hm-S-MR | 0.236 | 0.001 | 0 | $10^{-3}$ | -1.0 | 0.1 | $0.03\Omega_{0}$ |
L-hm-S-SR | 0.236 | 0.001 | 0 | $10^{-3}$ | -1.0 | 0.1 | $0.003\Omega_{0}$ |
Aspect ratio test runs (Section 5.4) |
NL-hs-S-NR | 0.236 | 0.8 | $10^{-3}$ | 10 | -1.0 | 0.05 | 0 |
NL-hs-R-NR | 0.236 | 0.5 | $10^{-4}$ | 1 | -1.0 | 0.05 | 0 |
NL-hs-V-NR | 0.236 | 0.5 | 0 | $10^{-3}$ | -1.0 | 0.05 | 0 |
NL-hl-S-NR | 0.236 | 0.8 | $10^{-3}$ | 10 | -1.0 | 0.15 | 0 |
NL-hl-R-NR | 0.236 | 0.5 | $10^{-4}$ | 1 | -1.0 | 0.15 | 0 |
NL-hl-V-NR | 0.236 | 0.5 | $0$ | $10^{-3}$ | -1.0 | 0.15 | 0 |

11footnotetext: Simulations mentioned in Appendix A and Appendix C are not
included.

22footnotetext: The unified name for parameter sruvey simulations, with
parameters being all combinations of those listed in this row, totaling 160
simulations.

33footnotetext: The resolution of this run is set to $N=2048$.

Note. — $\sigma_{\phi}$: shadow width parameter; $\epsilon$: shadow amplitude
parameter; $\alpha$: viscosity parameter; $\beta$: cooling rate parameter; p:
temperature slope; $h_{0}$: disk aspect ratio at inner boundary; $\Omega_{\rm
shadow}$: shadow procession angular frequency; $\Omega_{0}$: Keplerain angular
velocity at $r=1$. All runs, except for run L-hm-S-FR, use a resolution of
$N=512$.

###  2.3 Diagnostics of Substructures

As we will demonstrate, our simulations generate a variety of substructures of
all types. In this section, we provide the diagnostics we employ to identify
and characterize these substructures. To minimize the influence of the
boundaries and wave-killing, we restrict the analysis domain to be
$r\in[3,21]$.

Vortices exhibit as anti-cyclonic flows with pressure maxima at the center
that can potentially be strong dust traps. They are identified as regions with
negative vorticity, which is defined as $\nabla\times\delta\textbf{v}$ with
$\delta\textbf{v}$ being the difference between current fluid velocity and
background fluid velocity. We quantify individual vortices based on their mean
vorticity (normalized by background Keplerian angular velocity), density
contrast, spacing, and aspect ratio. In doing so, we first choose the vortex
boundary to be where the density is $10\%$ of the density at the vortex center
after subtracting background, while ensuring that the vorticity remains below
zero. This is motivated from the analytical work of Lyra & Lin (2013) while
being robust to the influence of density waves. In our simulations, vortices
are constantly generated and destroyed; only the largest vortices are chosen
(usually can survive for at least 100 local orbits). We measure the density
contrast by comparing the average density in vortex with the average density
at the same radius. The spacing of vortices is calculated by the radial
distances between neighboring vortices, which are normalized by the local
scale height at the midpoint radius between the two vortices. As vortices can
be highly time variable, all quantities are calculated and averaged over
several snapshots (see in Section 4).

In ring-forming disks, we measure the density contrast, width, spacing and
eccentricity of the rings. We identify the rings by first fitting the
background density as a power law, and consider peaks/troughs above/below the
fitted profile as ring peaks/gap centers. The boundaries of the rings are
identified as the radius at the midpoint between peak and valley densities,
with ring width being the distance between the two boundaries for each ring.
Density contrast is calculated by comparing the density between peaks and
boundaries. The final ring width is obtained by averaging the widths of all
identified rings, and each ring width is normalized to the local scale height
of the disk. Ring spacing is measured as the radial distances between the
boundaries of two neighboring rings, normalized in a way similar to that for
vortices, and averaged over several snapshots. In the above, we have treated
the rings as axsymmetric by working with 1D profiles, whereas in practice we
have found that the rings can be eccentric. For identified rings, we further
track the maximum density in 2D data and measure their eccentricity by fitting
an ellipse. Incomplete rings at the boundary of the analysis domain are
excluded from the statistics.

For spirals, we quantify their density contrast, number of spiral arms,
pattern speed and pitch angle. The density contrast is obtained by comparing
density of the spiral spine and the fitted background density at same radii.
In our simulations, we obtain the spiral phase at each radius using Fourier
decomposition and the pitch angle is obtained by fitting the phase angle with
a logarithmic function $\varphi=m(\tan\alpha_{p})^{-1}\ln r+\phi_{0}$, where
$\alpha_{p}$ is the pitch angle and $m$ is number of spiral arms. The constant
$\phi_{0}$ is further employed to measure the pattern speed of the spirals.

##  3 Representative Results

In this section, we present representative outcomes of shadow-driven
substructures at fixed disk aspect ratio $h$ before giving more comprehensive
statistical results. The three representative runs, denoted as“NL-hm-S-NR,”
“NL-hm-R-NR,” and “NL-hm-V-NR,” can be found in Table 1. We show snapshots of
major fluid quantities of interest (i.e.
$\Sigma,T,\textbf{v}_{r},\textbf{v}_{\phi},\nabla\times\textbf{v}$) from our
simulations, and discuss the results below.

###  3.1 Spirals

Figure 2: Density (a), temperature (b), radial velocity (c) and azimuthal
velocity (d) evolution in spiral-forming disk
($\sigma_{\phi}=0.236,\epsilon=0.8,\alpha=10^{-3},\beta=10,p=-1.0$). All
quantities except for $v_{r}$ are normalized by their initial values, and
$v_{r}$ is normalized by initial $v_{\phi}$. The color range for temperature
differs between the formation process of the early stage (first and second
rows) and the relative quasi-steady state (third to fifth rows). The white
solid lines represent the fitted spirals based on density contrast.

Spirals typically form in disks where the shadow is relatively weak, such as
those characterized by slow cooling or weak shadow amplitude. We choose
spirals formed in a disk with the following parameters as an example:
$\sigma_{\phi}=0.236,\epsilon=0.8,\alpha=10^{-3},\beta=10,p=-1.0$ (run NL-hm-
S-NR). As depicted in Figure 2, spirals form relatively quickly (first row in
Figure 2), typically within approximately 20 local orbits, and once formed,
they remain highly stable111The growth in density perturbations observed in
the last two rows of Figure 2 is primarily due to the combined effects of
strong viscous heating and the influence of wave damping zones.. These spirals
are clearly density waves, showing spiral patterns in all diagnostic physical
quantities. The spiral patterns are stationary (i.e., pattern speed is zero),
which is related to the fact that our shadow patterns have zero angular
velocity. Further discussions regarding the relationship between the
properties of the spirals and the other two substructures will be provided in
Section 5.2. In addition, by examining the second column of Figure 2, we see
that with inefficient cooling, the overall temperature is systematically
cooler than the initial temperature by $\sim 15\%$. The azimuthal varies
smoothly through the shadowing regions, with a maximum temperature variation
of about $4\%$.

###  3.2 Rings

Figure 3: Density (a), temperature (b), radial velocity (c) and azimuthal
velocity (d) evolution in ring-forming disk
($\sigma_{\phi}=0.236,\epsilon=0.5,\alpha=10^{-4},\beta=1,p=-1.0$). All
quantities except for $v_{r}$ are normalized by their initial values, and
$v_{r}$ is normalized by initial $v_{\phi}$. The black dashed lines represent
the fourth and fifth fitted rings in this disk.

The conditions for ring formation generally require either slow cooling or a
combination of moderate viscosity and shadow amplitude (for more detailed
information, see Section 4). In Figure 3, we adopt parameters
$\sigma_{\phi}=0.236,\epsilon=0.5,\alpha=10^{-4},\beta=1,p=-1.0$ (run NL-hm-R-
NR) to illustrate the typical formation process and properties of ring
structures.

–Formation. The formation of rings begins with the presence of two-arm spirals
following a transient period (as seen in the first to third rows of Figure 3).
The spirals appear only marginally stable, which later break apart and
reconnect to form concentric rings in surface density (as shown in the fourth
and fifth rows of Figure 3), which takes a relatively long time of $\sim 100$
local orbits. On the other hand, the spiral patterns remain in the velocity
structure even after ring formation, although they are distorted (as opposed
to the spirals discussed in Section 3.1 and could become a distorted ring
patterns in some cases).

–Evolution and main properties. Once formed, the amplitudes of the rings
continue to increase slowly, reaching a steady state over a few hundred local
orbits, where the gas density in rings are about $10\%$ higher than the
background. However, the density within one ring at quasi-steady state is
unevenly distributed, with surface density near the broken/reconnection
location being smaller, which will be further discussed in Section 5.2 and
Appendix A. The typical ring width is approximately twice the local scale
height, and the spacing is regular (about $4H$ between peaks of two
neighbouring rings) across the disk (further discussed in Section 4). We find
the rings to be eccentric (but centered on the star), with the eccentricity
measured to be $e\sim 0.12$. As can be inferred from the third and fourth
columns in Figure 3, the ratio $v_{r}/v_{\phi}$ is approximately $10^{-3}\ll
e$, suggesting that these rings do not directly correspond to the gas moving
in eccentric orbit. Also, the eccentric rings do not precess, analogous to
spiral patterns that remain stationary, thus coroborating the fact that the
rings emerge as the aftermath of spiral patterns. With moderate cooling, the
azimuthal temperature contrast reaches 15$\%$ and may cause azimuthal
brightness variations in observed rings, although we caution for our highly
simplified thermodynamic treatment (further discussed in Section 5.5).

###  3.3 Vortices and Crescents

Figure 4: Density (a), temperature (b), vorticity (c) and radial velocity (d)
evolution in vortex-forming disk
($\sigma_{\phi}=0.236,\epsilon=0.5,\alpha=0,\beta=0.001,p=-1.0$). The color
ranges for density, vorticity and radial velocity differ between the early
spiral formation stage (first and second rows) and the vortex-dominant stage
(third and fourth rows). The simulation duration increases sequentially from
top to bottom. The white frame delineates the boundary of each identified
vortex, and there are no significant structures inside the white frame in the
temperature plot (b) (compare with density (a), vorticity (c) and radial
velocity (d)) due to rapid cooling. A white star within each frame indicates
the center of the vortex. The red numbers highlight the selected vortices.

Crescents can be described as rings that exhibit an azimuthal variation in
intensity (Bae et al., 2023). Physically, the crescents discussed in this
paper are all induced by vortices, thus we use “vortices” and “crescents”
interchangeably. Figure 4 shows an example of shadow-driven formation of
vortices/crescents. This usually occurs with strong shadow amplitude and rapid
cooling, thus strong thermal forcing, and we adopt $\epsilon=0.5,\beta=0.001$
in this example (run NL-hm-V-NR).

–Formation. With rapid cooling, the disk temperature almost instantly relaxes
to the target temperature both within and outside of the shadow region,
resulting in a 50$\%$ variation in azimuthal temperature given our setup. This
leaves two symmetric low-pressure regions that form quickly at the shadow
locations. In the initial stages (first and second rows in Figure 4), it leads
to the appearance of spiral features in surface density. With strong thermal
forcing that constantly perturbing the disk, the system subsequently becomes
more chaotic (third row in Figure 4) where the velocity field undergoes
significant alterations. Although the physical process is not entirely clear,
vortex/crescent formation ensues, as identified in fourth row of Figure 4.
Selected vortices and crescents are marked by white frames in Figure 4.

–Evolution and main properties. Shadow-driven vortices are all anti-cyclonic
in nature, which can be observed either from the negative vorticity (the 3rd
column of Figure 4) or from the change in the sign of radial velocity across
the vortex center (changing from negative to positive when viewed along the
direction of rotation (counterclockwise), as shown in the 4th column of Figure
4). We observe that vortices started small and are continuously generated.
They merge to form larger ones under the influence of differential rotation
within approximately 60 local orbits, ultimately manifesting as relatively
large crescent-shaped structures. In Figure 4, vortices labeled as 4a and 4b
are undergoing a merger into one single vortex. We find that these patterns
largely corotate with the gas, as expected, and their azimuthal locations are
found to be largely random, with no preference to stay in or out of the
shadows. The disk gas remains turbulent and chaotic throughout the evolution
due to strong perturbations from thermal forcing. Velocity deviations from
local Keplerian inside the vortex region are around $0.5c_{s}$. Additionally,
the local level of turbulence, measured in terms of root mean square (rms)
velocity fluctuations averaged in azimuth, is approximately $10\%$ of the
local sound speed. The typical aspect ratio of the vortices/crescents is about
6, with their density contrast being 1.4. The normalized vorticity in this
case is 0.2. Despite of modest to strong level of turbulence, the large
vortices are relatively long-lived, with typical lifetime of at least 300
local orbits.

##  4 Statistics of Substructures

To gain deeper insights into the dynamical consequences of shadows, we
conducted a comprehensive exploration of parameter space. We performed a total
of 160 simulations (run S-h-all), encompassing a wide combination of
parameters. Most results show similarities with one of the aforementioned
three representative cases. We thus primarily summarize the outcomes in a
statistical manner.

For simulations that exhibit the formation of rings and spirals, we only
measure their properties at the end of the simulations when the system has
already reached a steady state. For simulations with vortex/crescent formation
which are generically chaotic, we select four specific snapshots, denoted as
$P_{\rm orb1}=5000P_{0}$, $P_{\rm orb2}=10000P_{0}$, $P_{\rm
orb3}=15000P_{0}$, and $P_{\rm orb4}=20000P_{0}$. The statistical values for
vorticity, density contrast, spacing, and aspect ratio of the vortices are
calculated by averaging the results at these snapshots.

The simplified statistical results are presented in Figure 5, and more
detailed ones are provided in Figures 14 and 15. It is important to emphasize
that panels shaded with red or blue lines are actually undergoing a vortex-
ring transition or a ring-spiral transition state (see discussion in Appendix
A).

Generally speaking, shadows are capable of generating different kinds of
substructures under different parameter settings. Additionally, we found that
the dominant form of shadow-driven substructures changes from spirals to rings
and eventually becomes vortices/crescents as cooling timescales and/or
viscosity decreases. Where exactly the transition occurs depends on other
parameters such as the shadow amplitude, width, and disk aspect ratio, etc.,
and these will be discussed in more detail in the following subsections.

###  4.1 Statistics for Spirals

Two-arm spirals are fundamental substructures in our simulations, dominating
in disks with cooling timescales significantly longer than the dynamical
timescale, high viscosity ($\alpha>10^{-3}$), or very weak shadow amplitude
(see Section 5). Here, we focus on discussing their density contrast, pattern
speed, and pitch angle.

–Density contrast. In general, stronger thermal forcing, higher shadow
amplitude, wider shadow width, etc. leads to stronger density contrast in the
spirals. However, as the spiral-dominated regime generally requires weak
thermal forcing, the spiral amplitudes are typically low (with upper limit
only $1\%$ higher than background density at the same radius).

–Pattern speed. Spirals found in our simulations are density wave patterns
with zero pattern speed, which also results in non-precessing rings. More
generally, the spiral pattern speed exactly matches the shadow’s pattern
speed, which will be further discussed in Section 5.3.

–Pitch angle. The pitch angle is solely affected by the disk aspect ratio.
With weak thermal forcing, we consider the dispersion relation of spiral
density waves in the linear regime under the WKB approximation (Lin & Shu,
1964)

$m^{2}(\Omega_{p}-\Omega)^{2}=k^{2}c_{s}^{2}+\kappa^{2}.$ (15)

Here, $\Omega_{p}$ represents the spiral pattern speed, $k$ is the radial wave
number, and $\kappa\approx\Omega_{K}$ is the epicyclical frequency. The spiral
pitch angle can be estimated by $\alpha_{p}=\partial r/(r\partial\phi)\approx
m/(|k|r)$. With $\Omega_{p}=0$ and $m=2$, we obtain
$\alpha_{p}\sim\frac{2}{\sqrt{3}}h=constant$ for $p=-1$ disks and
$\alpha_{p}\sim\frac{2}{\sqrt{3}}h_{0}r^{0.25}$ for $p=-0.5$ disks. Taking the
disk parameters used in our simulations (with $h_{0}=0.1$) and averaging over
radius gives $\alpha_{p}=6.6^{\circ}$ for disks with $p=-1$ and
$\alpha_{p}=12.1^{\circ}$ for disks with $p=-0.5$. These estimated values
agree well with our simulation results, which we find to be
$7.344_{-0.535}^{+0.607}$∘ and $13.202_{-1.766}^{+2.16}$∘ (see Figure 15),
respectively.

Figure 5: Statistics of shadow-driven substructures based on 80 of 160
simulations with constant disk aspect ratio ($h=0.1$). The figure is divided
into two parts by a dashed line, representing the shadow width being 45
degrees on the left and 15 degrees on the right. Within each part, the left
column corresponds to disks with $\epsilon=0.5$, while the right column
corresponds to disks with $\epsilon=0.8$. Each subfigure in the
$\beta$-$\alpha$ sections represents a specific combination of parameters.
Within each $\beta$-$\alpha$ section, there are three rows representing the
dominant structures in the disk: vortices/crescents, rings, and spirals from
top to bottom. Each type of structure is represented by a different type of
square marker, whose color represent the density contrast of the
substructures. The figure also includes red and blue line shaded areas,
indicating disks undergoing transitions from vortex-ring and ring-spiral
phases, respectively. Figure 6: The average and error of normalized ring
spacing. The black solid line represents simulations with $p=-1.0$, while the
red solid line represents simulations with $p=-0.5$. Points marked with an ’x’
do not have error bars. In this case, there are at most two selected rings
within the disk. Upward arrow represents lower limit points, where there is
only one ring within the disk. Figure 7: The average and error of normalized
vortex/crescent spacing. The black solid line represents simulations with
$p=-1.0$, while the red solid line represents simulations with $p=-0.5$.
Points marked with an ’x’ do not have error bars. In this case, there are at
most two selected vortices/crescents within the disk. Upward arrow represents
lower limit points, where there is only one selected vortex/crescent within
the disk.

###  4.2 Statistics for Rings

In our simulations, rings dominate in disks with cooling timescales comparable
to the dynamical timescale ($\beta\sim 1$) when $\alpha$ is roughly below
$10^{-3}$. For much higher viscosity, rings dominate even when the cooling
rate approaches the isothermal limit ($\beta=10^{-3}$). Typically, this value
is $\alpha=10^{-2}$ for disks with $\sigma_{\phi}=0.236$ and $\alpha=10^{-3}$
for disks with $\sigma_{\phi}=0.079$. Overall, the parameter space for the
dominance of rings is modest thermal forcing, in between the cases that form
vortices/crescents (strong forcing, see next subsection) and spirals (weak
forcing). In fact, we pose that rings can be viewed either as “reconnected
spirals” (stated in Section 3.2), or “failed vortices”, where the latter
connection arises from the finding that vortex-ring transitions often involve
crescents with very large aspect ratios, although the boundary between this
transition is not necessarily clear-cut, and will be further discussed in
Appendix A. Below, we will focus on “normal” rings (not under transition), and
will discuss the density contrast, ring radial width, ring spacing,
eccentricity, and the parameters that have strong influence on them.

–Density contrast. As shown in Figure 5 and 14, gas densities are typically
$1-20\%$ higher than the background density in ring-dominant disks, and ring
density contrast is enhanced by larger shadow amplitude and width. Density
contrast could reach very small values, such as $0.3\%$, in the ring-spiral
transition, and very large values, such as $50\%$, in the vortex-ring
transition.

–Width and spacing. The ring widths in our simulations are usually 2 times the
local scale height, regardless of shadow parameters. Similarly, for almost all
cases, the spacing between neighboring rings is approximately $4H$, as
depicted in Figure 6. There is very small deviations from the mean, indicating
a highly uniform distribution of rings within the disk.

–Eccentricity. As will be stated in Section 5, ring structures are generated
following the “reconnection” of two-armed spirals in the early stages of disk
evolution, causing the ring to become eccentric with zero pattern speed (as
shadows are stationary). More flared disk morphology results in larger spiral
pitch angles, making the spirals less tightly wound. As a result, the rings
formed in this case tend to be more eccentric. Additionally, we find that
viscosity has a strong impact on eccentricity. Typically, ring eccentricity
varies from 0.1 to 0.7 as $\alpha$ increases from $0$ to $10^{-2}$ in our
simulations (see Figure 15 for details). The angle between the ring’s major
axis and the effective shadow center (e.g. $\phi=0^{\circ},180^{\circ}$ when
$\beta=0.001$) is typically between $80^{\circ}$ and $110^{\circ}$.

###  4.3 Statistics for Vortices and Crescents

As we mentioned in Section 3.3 and better seen in Figure 5, vortices/crescents
tend to dominate in disks characterized by fast cooling processes ($\beta<1$),
low viscosity ($\alpha=0,10^{-4}$), high shadow amplitudes ($\epsilon=0.8$),
and wide shadow widths ($\sigma_{\phi}=0.236$). Such parameter settings all
point to strong thermal forcing. Below, we discuss the properties of the
shadow-driven vortices/crescents, focusing on density contrast, spacing and
aspect ratio of vortices/crescents under the influence of these parameters.

–Vorticity and density contrast. The density contrast of substructures is a
crucial factor as it directly influences their detectability. From our
explorations, the density of the crescents are typically $10-50\%$ higher than
the average density at same radius for all vortex-dominated disks. The density
contrast is generally slightly higher for stronger shadow intensity, larger
shadow width, and faster cooling, but the trend is not definitive given the
chaotic nature of the system. The normalized vorticity ranges from 0.1 to 0.6
in vortex-dominated disks, with vorticity around 0.2 in most cases,
potentially reaching up to 0.6 in the most extreme cases (large $\epsilon$ and
$\sigma_{\phi}$). No clear relationship is found between vorticity and density
contrast due to the high turbulence level, which is around $0.1c_{s}$. The
velocity deviations from local Keplerian inside the vortex region ranges from
0.4 to 1.2 $c_{s}$, indicative of strong rotation in the vortices.

–Spacing. The statistical results of the spacing of vortices/crescents are
plotted in Figure 7. In all simulations, the distance between neighboring
vortices/crescents is typically between $2H$ and $4H$. The spacing is less
uniform compared to rings, and is related to the fact that vortex-dominated
disks are usually turbulent. Note that the small error bars in a few cases are
related to very limited number of vortices/crescents (2 or 3); the lower limit
point represents the case where there is only one vortex-induced crescent in
the disk. Similar to the case shown in Figure 4, the azimuthal locations of
the vortices/crescents are largely random with no direct correlation with the
position of the shadows.

–Aspect ratio. The aspect ratio of crescents/vortices is less affected by
different parameters. Typically, in vortex-dominated disks, this value is
about 6. However, for cases close to (for example, $\sigma_{\phi}=0.079$,
$p=-1$, $\epsilon=0.5$, $\alpha=0$, $\beta=0.01$) or undergoing (for example,
$\sigma_{\phi}=0.236$, $p=-0.5$, $\epsilon=0.8$, $\alpha=10^{-4}$, $\beta=1$)
the vortex-ring transition in parameter space, the aspect ratio can be very
large (greater than 12). More detailed results are shown in Figure 15 in
Appendix B.

##  5 Discussion

Figure 8: Density (a), temperature (b), radial velocity (c) and azimuthal
velocity (d) evolution in linear evolution disk
($\sigma_{\phi}=0.236,\epsilon=0.001,\alpha=0,\beta=0.001,p=-1.0$). All
quantities except for $v_{r}$ are normalized by their initial values, and
$v_{r}$ is normalized by initial $v_{\phi}$.  Figure 9: Formation processes of
substructures from linear to non-linear regime. All plots show density
contrast in disks in logarithmic scale.

In this paper, we have conducted simple numerical experiments to study the
dynamical consequence of shadows cast from the inner disk to the outer disk as
a result of thermal forcing. We have restricted ourselves to a small number of
parameters, and the discussion has been largely phenomenological. In this
section, while not going into full detail, we conduct additional studies to
help better understand the origin and trend of shadow-driven substructures,
and briefly discuss their potential implications.

###  5.1 Linear regime

Based on the analysis and discussions in the previous sections, here we
provide further analysis to gain better physical insights on the shadow-driven
substructure formation. As we observe that in all cases, substructure
formation starts from the formation of two-armed spirals under our shadow
prescriptions. This suggests that spirals are the most fundamental form of
shadow-driven substructure, and it can be instructive to look into how spirals
form and evolve under very weak thermal forcing to avoid nonlinear effects. We
thus further conducted a series of 2D inviscid hydrodynamic simulations with
varying perturbation strengths ($\epsilon=0.001$, $\epsilon=0.01$,
$\epsilon=0.1$) while keeping the cooling timescales consistent
($\beta=0.001$). Without viscosity, the simulations are in hydrostatic
equilibrium to start with before thermal forcing is introduced.

In Figure 8, we present the results from the $\epsilon=0.001$ simulation (run
L-hm-S-NR in Table 1). When the shadow is introduced, gas flows into the
shadowed region in a counterclockwise manner. The gas between the shadow
center (pressure minimum) and its rear edge, i.e., between $315^{\circ}$ and
$360^{\circ}$ in the first row of Figure 8, gets accelerated, while the gas
between the shadow center and its leading edge, i.e., between $0^{\circ}$ and
$45^{\circ}$ in the first row of Figure 8, gets decelerated. This leads to gas
piling up near the shadow center, while the neighboring gas is slightly
rarefied, which naturally launch density waves.

As the disk evolves, such density waves wind up due to differential rotation
(see second row of Figure 8). In the meantime, the periodic forcing at the
shadow location continues, keep launching new density waves, leading to
interference. After a few local orbits, the system reaches a relatively steady
pattern of two-arm spirals (see third and forth rows of Figure 8), which
remain stable over long-term. The spirals share the same pattern speed of the
shadows (in this case, zero), and the pitch angle also remains unchanged. We
note that this is very different from planet-induced spirals in that a planet
launches density waves through discrete Lindblad resonances, while as shadows
are cast over a wide range of radii, each radius can excite its own density
waves. In our case, the pattern speed of the shadow is zero, and the only
relevant resonance condition is simply given by $\Omega=\kappa/m$, where
$m=1,2,\cdots$. However, taking $m=2$, we see that with $\kappa\approx\Omega$
for Keplerian disks, no resonance condition is satisfied. In other words, the
two-armed spirals are not driven by Lindblad resonances, but are the effective
eigen-state of thermally-forced oscillations.

###  5.2 Towards the nonlinear regime

We note that even in the linear regime, the spiral patterns are distorted due
to thermal forcing. These can be most easily seen from the velocity
perturbations in the last three columns of Figure 8. They are also present in
the density perturbations where the amplitude of the spirals varies across the
shadow region. The form of the distortion can depend on system parameters,
which is found to be different in Figure 2 where cooling time is significantly
longer. We speculate that such distortions are the source of instability when
thermal forcing enters the nonlinear regime.

Based on our discussions in the previous sections, we summarize the formation
of shadow-driven substructures in Figure 9. Irrespective of whether thermal
forcing is linear or nonlinear, the initial phase of the development is
similar, involving the formation of two-armed spirals, as shown in (a)-(b).
The spirals persist under linear and weakly nonlinear thermal forcing, as seen
in the “linear branch” and “spiral branch” in (c)-(f). The properties of
spirals are similar between the linear and weakly nonlinear regimes, in terms
of pitch angle and pattern speed.

When thermal forcing becomes slightly stronger, the spiral arms undergo a
relatively quiescent transformation by “reconnecting” into eccentric rings
(see Figure 9(g), (h)). The eccentricity of these rings is largely set by the
pitch angles of the original two-arm spirals stage and disk viscosity.
However, when the thermal forcing becomes too strong, the spirals break in a
highly chaotic manner (see Figure 9(i), (j)), leading to the formation of more
localized vortices/crescents.

Figure 10: Normalized density perturbation for simulations with
$\Omega_{\text{shadow}}=1\Omega_{0}$ (10), $0.03\Omega_{0}$ (10), and
$0.003\Omega_{0}$ (10) at the final steady state under the parameters setting
as $\sigma_{\phi}=0.236,\epsilon=0.001,\alpha=0,\beta=0.001,p=-1.0$. The radii
of corotation resonances (CR), inner Lindblad resonances (ILR), and outer
Lindblad resonances (OLR) are shown as yellow, green, and purple dashed lines
respectively.

###  5.3 Rotating shadows

In this paper, we have only discussed the situation when the shadow’s pattern
speed is zero. However, if the misaligned inner disk precesses around the
central star, the shadow cast from the inner region would have a pattern
speed, which then changes the resonance condition discussed in Secton 5.1. To
extend our study to more general conditions, we have conducted additional
simulations with rotating shadows in the linear regime, with three different
shadow pattern speeds, $\Omega_{\text{shadow}}=1\Omega_{0}$ (run L-hm-S-FR),
$0.03\Omega_{0}$ (run L-hm-S-MR), and $0.003\Omega_{0}$ (run L-hm-S-SR). Here,
$\Omega_{0}$ is Keplerain angular velocity at $r=1$. The detailed parameter
settings can be found in Table 1. The density structure from these simulations
in the final states are shown in Figures 10. We measure the pattern speed of
the spirals $\Omega_{p}$ in these situations, and we confirm that in all three
cases, the spirals all have $\Omega_{p}=\Omega_{\text{shadow}}$.

Given the pattern speed, the radii of corotation resonances (CR), inner
Lindblad resonances (ILR), and outer Lindblad resonances (OLR) can be
calculated by $\Omega_{p}=\Omega$, $\Omega_{p}=\Omega\pm\kappa/m$ (with
$m=2$), and shown as yellow, green, and purple dashed lines in Figures 10.
With the WKB dispersion relation 15, the permitted regions for density wave
propagation are outside the Lindblad resonances. In the fast-rotating case
$\Omega_{p}=\Omega_{0}$, density waves are permitted beyond the OLR, and the
spirals are tightly wound towards outer radii with pitch angle $\alpha_{p}\sim
h(\Omega/\Omega_{p})\sim h(r/r_{0})^{-3/2}$. Even with a resolution of
$N=2048$ in Figure 10, it is still insufficient to resolve the spirals across
the entire disk, weakening the spirals at the outer disk by numerical
dissipation. With intermediate $\Omega_{\rm shadow}=0.03\Omega_{0}$, the ILR
and OLR are located at $r=6.5$ and $r=13.5$, respectively. Clearly, there are
well-defined spirals outside the Lindblad resonances, which break inside the
Lindblad resonances. In the slow-rotating case with $\Omega_{\rm
shadow}=0.003\Omega_{0}$, even the ILR is beyond the computational domain, and
the results are largely identical to the stationary case described in Section
5.1.

Given the discussion above, we expect the results presented in this paper
largely applies to regions inside the ILR in slowly-precessing shadows.
Although not the focus of this paper, it is worth noting the significance of
moderately rotating shadows, where the corotation radius lies within the disk
region. Our findings are morphologically similar to those of Montesinos &
Cuello (2018), who demonstrated that the morphology of shadow-driven spirals
notably resembles the planetary wakes caused by embedded planets in the disc
using radiative transfer. For better comparison with planet-induced spirals,
more detailed investigation with more realistic physics (especially dust and
radiative processes) is necessary for the slow-rotating case, especially in
regions between the ILR and OLR.

###  5.4 Dependence on disk aspect ratio

Figure 11:  Normalized density perturbation for runs with different disk
aspect ratios. The parameters are identical to those in the representative
runs (NL-hm runs) in Section 3, except for the disk aspect ratio. Detailed
parameter settings are provided in Table 1. (a): snapshot with $P_{\rm
orb}=20000P_{0}$ for run NL-hs-S-NR
($h_{0}=0.05,\sigma_{\phi}=0.236,\epsilon=0.8,\alpha=10^{-3},\beta=10,p=-1.0$).
(b): snapshot with $P_{\rm orb}=20000P_{0}$ for run NL-hs-R-NR
($h_{0}=0.05,\sigma_{\phi}=0.236,\epsilon=0.5,\alpha=10^{-4},\beta=1,p=-1.0$).
(c): snapshot with $P_{\rm orb}=20000P_{0}$ for run NL-hs-V-NR
($h_{0}=0.05,\sigma_{\phi}=0.236,\epsilon=0.5,\alpha=0,\beta=0.001,p=-1.0$).
(d): snapshot with $P_{\rm orb}=3000P_{0}$ for run NL-hl-S-NR
($h_{0}=0.15,\sigma_{\phi}=0.236,\epsilon=0.8,\alpha=10^{-3},\beta=10,p=-1.0$).
(e): snapshot with $P_{\rm orb}=10000P_{0}$ for run NL-hl-R-NR
($h_{0}=0.15,\sigma_{\phi}=0.236,\epsilon=0.5,\alpha=10^{-4},\beta=1,p=-1.0$).
(f): snapshot with $P_{\rm orb}=1500P_{0}$ for run NL-hl-V-NR
($h_{0}=0.15,\sigma_{\phi}=0.236,\epsilon=0.5,\alpha=0,\beta=0.001,p=-1.0$).

In the preceding discussion, we observed that shadow-driven substructures are
closely tied to thermal forcing, which is influenced not only by the cooling
process but also by disk temperature. Additionally, detailed characteristics
of substructures, such as pitch angle or eccentricity, are affected by the
disk aspect ratio h. Therefore, it is natural to further investigate the
influence of $h_{0}$. We conducted additional simulations with $h_{0}$ ranging
from 0.03 to 0.15, focusing on $h_{0}=0.05$ and $h_{0}=0.15$. These
simulations, denoted as NL-hs-S-NR, NL-hs-R-NR, NL-hs-V-NR and NL-hl-S-NR, NL-
hl-R-NR, NL-hl-V-NR respectively, maintained the same parameters as the
representative runs discussed in Section 3 except for $h_{0}$ (see Table 1).
We note that here “S”, “R”, and “V” do not necessarily indicate dominant form
of substructures but rather serve to guide the reader that these runs only
vary $h_{0}$ compared to representative runs.

In the NL-hs run series ($h_{0}=0.05$), with lower target temperature, we see
that the NL-hs-S-NR (Figure 11) and NL-hs-R-NR (Figure 11) runs maintain
spirals and rings as the primary substructure, respectively. We see the
spirals are more tightly wound and the rings spacing remains uniform except
for being smaller. The changes are exactly in proportion to $h_{0}$, and the
general properties of the rings and spirals are otherwise identical to those
discussed in the NL-hm runs. For the NL-hs-V-NR run, while the vortices are
clearly the dominant, many of the overdensities close a full circle, and we
identify this run as in the vortex-ring transition state.

In the NL-hl run series ($h_{0}=0.15$), with higher target temperature, we see
that all three NL-hl runs retain their spiral, ring and crescent/vortex as the
dominant substructure, respectively. Similarly, the spirals are more open, the
rings are more eccentric, and the vortices are larger and more widely spaced,
as expected.

Overall, we find that varying $h$ slightly alters the boundary where different
forms of substructures dominate, while the general properties for individual
substructures largely remain consistent with what we have found in the
fiducial simulations with $h_{0}=0.1$.

###  5.5 Observational implications

Given the diverse dynamical consequence of shadowing, such disks is expected
to exhibit a variety signatures that are potentially observable. However, it
should be noted that our work serves as a general exploration without detailed
modeling, including radiation transport, dust dynamics, shadow precession
rates (e.g., Pinilla et al., 2015; Stolker et al., 2016; Wolff et al., 2016;
Debes et al., 2017), and realistic shadow morphologies may differ from our
prescription (e.g., Muro-Arena et al., 2020; Debes et al., 2017).
Additionally, as there are a variety of other mechanisms that can drive
substructures (e.g., see reviews by Andrews, 2020; Bae et al., 2023; Benisty
et al., 2023), such as planet-disk interactions, icelines etc. Our shadowed
disk simulations implicitly assumed a smooth disk to start with, and it is
conceivable that the final outcome is set by the interplay between the
existing substructures and shadowing. Besides such dynamical interplay,
substructures themselves can self-shadow (e.g., Zhang et al., 2021), which can
further complicate the situation. Therefore, a systematic observational
comparison with specific sources is beyond the scope of this work. Below, we
mainly discuss general aspects of potential observational implications.

–Spirals. Spirals generated from shadows may not be easily detectable in the
submm continuum or in kinematics, but may be observable in scatter light.
Nearly all spiral-dominant disks correspond to weak thermal forcing, resulting
in only about $0.1\%$ higher gas density than the background. This not only
makes pressure variations across the spirals small that is difficult for dust
trapping, and only sufficiently small particles with a stopping time shorter
than the spiral crossing time (typically requiring the Stokes number much less
than 0.1) can potentially be trapped by the spiral (e.g. Sturm et al., 2020;
Speedie et al., 2022). With the weak spirals, the gas velocity is found to
show very small deviationsfrom Keplerian ($\sim 0.1\%$ $v_{k}$, as opposed to
$\gtrsim 0.5\%$ $v_{k}$ for typical ALMA observations (Pinte et al., 2023).),
making it difficult for kinematic detection. On the other hand, such spirals
may be detectable in scattered light, as suggested by Montesinos et al. (2016)
for the HD 142527 disk, thanks to azimuthal variation of disk scale heights
across the spirals, though three-dimensional simulations are needed for proper
characterization.

–Rings. For full disks, our simulations predict the presence of multiple gas
rings that are uniformly spaced and weakly eccentric. The relatively high
density contrast in our simulations suggests that these rings likely
concentrate dust, making them readily observable in sub-mm wavelengths. While
the resulting dust rings formed are also likely uniformly spaced, whether they
can be eccentric remains uncertain (as the eccentric gas ring is a pattern and
does not reflect real motion), requiring simulations incorporating dust
dynamics. From all simulations, we find that the azimuthal temperature
contrast in the ring-dominant disks are typically greater than $8\%$ and can
reach up to $50\%$ as they approach to vortex-ring transition in disks with
high viscosity and rapid cooling. Such azimuthal temperature variations should
result in azimuthal brightness variations in the mm continuum image, which
however has not been revealed in in real shadowed disks with rings (e.g., HD
143006). This suggests that thermal forcing by shadows in these systems are
likely not as strong as given in our prescriptions, but we caution that
without detailed modeling of shadow morphology, radiation transport and dust
dynamics, we cannot make specific predictions for individual systems. On the
other hand, we comment that both the weakly eccentric ring pattern and low-
level of azimuthal temperature variation, if present, may affect the
interpretation of azimuthal asymmetries seen in multi-ring systems (e.g. Doi &
Kataoka, 2021; Liu et al., 2022). Finally, we note that detection by kinematic
signatures, with velocity disturbances being $\sim 1\%$ of the Keplerian
velocity, is possible but challenging since they are close to ALMA’s detection
limits.

–Crescents. Vortices generate significant velocity perturbations and are
favored sites for dust trapping. Given the adopted turbulent viscosity
parameter $\alpha_{t}\gtrsim 10^{-2}$ in most vortex-dominated simulations,
dust with Stokes number $St>\alpha_{t}\sim 10^{-2}$ is expected to concentrate
inside vortices overcoming turbulent diffusion (Birnstiel et al., 2013), and
can be readily observable in sub-millimeter wavelengths (Zhu et al., 2014).
Previous studies have found that detecting kinematic signatures of vortices
can be possible but challenging (Huang et al., 2018), despite the relatively
large vorticity (typically around $0.2$) and significant velocity deviations
from local Keplerian ($\delta v$ up to $1.2c_{s}$) inside vortex region . It
is expected that sources with modest inclination favors detection but requires
long integration time with ALMA (more than 10h) to achieve the necessary
signal-to-noise ratio.

##  6 Summary and Future Prospects

In this work, we have systematically studied the dynamical consequence of
thermal forcing by shadows cast to the outer protoplanetary disks. With a
large survey of parameters, we have identified a diverse forms of
substructures generated by shadows and studied their trends under different
thermodynamic and viscosity prescriptions. Our results apply in regimes where
the shadow is static or slowly-rotating (prograde), so that the corotation
radius is further than regions of interest. The main findings of our studies
are as follows.

1\. Two-arm spirals with identical pattern speed as the shadow are fundamental
substructures generated by weak thermal forcing ($\epsilon<=0.5$,
$\sigma_{\phi}=0.079$, $\beta>1$) or high viscosity ($\alpha>10^{-3}$). They
represent linear response to thermal forcing, and their pitch angle well
agrees with standard density waves. Both the density contrast (0.1-1$\%$
higher than background) and velocity disturbance up to 0.5$\%$ $v_{K}$) are
small and scale with the strength of thermal forcing.

2\. Disks with moderate thermal forcing are dominated by ring-like
substructures. In this regime (parameter space in between crescent/vortex and
spiral-dominant disks), the gas density contrast reaches 1-20$\%$ above the
background. The rings are uniformly spaced ($\Delta r/H\sim 4H$) and exhibit
pattern eccentricities on the order of $h/r$ or higher which rotates at the
same rate of the shadow.

3\. Crescents/vortices dominate disks under strong thermal forcing
($\epsilon>0.5$, $\beta\lesssim 0.1$, $\sigma_{\phi}=0.236$) and low viscosity
($\alpha<=10^{-4}$). In this case, the density contrast is typically 10-50$\%$
higher than the average density at the same radius. The vortices in our
simulations exhibit relatively large vorticity (ranging from 0.1 to 0.6,
typically around 0.2) and significant velocity deviations from local Keplerian
inside the vortex region (ranging from 0.4 to 1.2 $c_{s}$). Due to the chaotic
nature (local turbulence level is 0.1 $c_{s}$) of the vortex-dominant disk,
these structures are not uniformly spaced, with $\Delta r/H$ between 2 and 4.

4\. Thermodynamics and viscosity significantly influence the formation of
shadow-driven disk substructures. The dominant substructure transitions from
spirals to rings and eventually to vortices as cooling timescales and/or
viscosity decrease.

5\. Owing to the simplicity of our problem setup, it is premature to
definitely assess the observability of such shadow-driven substructures. We
anticipate that the azimuthal brightness contrast in the sub-mm continuum to
offer important constraints on the strength of the thermal forcing, while
detecting in-plane kinematic signatures is likely challenging.

Through our suite of physically-motivated while highly simplified simulations,
we highlight the importance on the dynamical impact of shadows or more
generally, inhomogeneous stellar irradiation, on the gas dynamics of PPDs
through thermal forcing. Given the fact that shadows are often observed in
scattered light images of disks, our results call for proper consideration and
incorporation of such effects for adequate modeling of such systems.

Our simulations can be considered as a starting point to understand the
dynamical effects of shadows on PPDs, yet real systems are likely much more
complex. This leaves several aspects to be considered and tested in the
future. Proper characterizing disk thermodynamics is a pre-requisite to
accurately model thermal forcing from shadows, which requires better modeling
of the shadow geometry, together with self-consistent radiation transport.
Such modeling under typical disk parameters (that are likely nearly optically-
thin) will likely reduce the azimuthal temperature contrast due to in-plane
radiation transport. Incorporation of dust dynamics is essential to obtain
dust response to the shadow-driven substructures. Such simulations are
expected to link the results with specific sources, as we are aware of efforts
underway (Ziampras et al., in preparation). We have also assumed the shadows
are cast to a full disk, whereas shadows are also observed in transition
disks, and it is also pertinent to account for the interplay other physical
mechanisms that cause disk substructures, with additional effect of self-
shadowing.

Finally, all existing studies of shadow-driven disk dynamics are conducted in
2D in the disk plane, whereas the shadow-driven thermal forcing is also
expected to also drive oscillations in the vertical direction (Liu & Bai, in
preparation). Future studies should incorporate 3D effects, which is essential
to further assess the fidelity of 2D simulation results, and make more
realistic observational predictions and comparisons.

Acknowledgements  We thank Yanqin Wu and Shangjia Zhang for useful
discussions, Pinghui Huang for helpful instructions on problem setup, and
Alexandros Ziampras for constructive exchanges. This work is supported by
National Science Foundation of China under grant No. 12233004, 12325304. We
also acknowledge the China Center of Advanced Science and Technology for
hosting the Protoplanetary Disk and Planet Formation Summer School in 2022
when this work was initiated. Numerical simulations are conducted in the Orion
and Sirius clusters at Department of Astronomy, Tsinghua University and
TianHe-1 (A) at National Supercomputer Center in Tianjin, China.

##  Appendix A Transition State

The vortex-ring transition represents the parameter regime where both the
features of vortices/crescents and rings can be observed in the disk. Four
examples of vortex-ring transition are illustrated in Figure 12. They are
recognized as vortex-ring transitions generally based on two reasons: rings
and vortices/crescents are simultaneously present in the disk (Figure 12), or
the basic morphology appears as rings but with significant asymmetry (Figure
12, 12, 12). In Figures 14 and 15, the left side of vortex-ring transition
cases depicts vortex-dominated disks, while the right side illustrates ring-
dominated disks. Further decreases in $\beta$ or $\alpha$ lead to the disk
being completely dominated by vortices/crescents.

The ring-spiral transition represents the parameter regime where both the
features of rings and spirals can be identified in the disk. Four examples of
ring-spiral transitions are shown in Figure 13. They either exhibit regularly
broken rings (Figure 13 and 13) or clearly display both rings and spirals
within the same disks (Figure 13 and 13). These transition regions lie between
ring-dominated disks and spiral-dominated disks in Figure 14 and 15. The disk
becomes dominated by spirals as $\beta$ or $\alpha$ increases.

From the transition states shown in Figure 12 and 13, we can verify that rings
exhibit characteristics of both vortices/crescents and spirals, as discussed
in Section 4.2. Slightly excessive thermal forcing, relative to ring-dominant
disks, can hamper reconnection (Figure 12) mentioned in Section 5.2, leading
the disk into a vortex-ring transition state with strongly asymmetric rings
(Figure 12) or crescents with large aspect ratios (Figure 12). Conversely,
with weak thermal forcing, the breaking of two-armed spirals is partial
(Figure 13), placing the disk into a ring-spiral transition state.

Figure 12: The demonstration of vortex-ring transition situation. The
parameters of these 4 plots are
$\sigma_{\phi}=0.236,p=-0.5,\epsilon=0.5,\alpha=0.0001,\beta=1$;
$\sigma_{\phi}=0.236,p=-0.5,\epsilon=0.8,\alpha=0.01,\beta=0.001$;
$\sigma_{\phi}=0.236,p=-1.0,\epsilon=0.8,\alpha=0.01,\beta=0.001$;
$\sigma_{\phi}=0.079,p=-0.5,\epsilon=0.8,\alpha=0,\beta=0.1$ respectively. All
these cases can be found in Figure 14 and Figure 15 as vortex or ring blocks
covered by red hatch lines.

Figure 13: The demonstration of ring-spiral transition situation. The
parameters of these 4 plots are
$\sigma_{\phi}=0.236,p=-0.5,\epsilon=0.8,\alpha=0,\beta=10$;
$\sigma_{\phi}=0.079,p=-1.0,\epsilon=0.8,\alpha=0.001,\beta=1$;
$\sigma_{\phi}=0.079,p=-0.5,\epsilon=0.5,\alpha=0,\beta=10$;
$\sigma_{\phi}=0.079,p=-0.5,\epsilon=0.5,\alpha=0.001,\beta=1$ respectively.
All these cases can be found in Figure 14 and Figure 15 as ring or spiral
blocks covered by blue hatch lines.

##  Appendix B Simulation statistics

The detailed statistical plot of vorticity and density contrast (Figure 14),
along with other parameters (Figure 15) of substructures, is presented here.
These two figures share the same structures. Each of these figures is divided
into two sections by a dashed line, representing shadow ranges of 45 degrees
($\sigma_{\phi}=0.236$) and 15 degrees ($\sigma_{\phi}=0.079$), respectively.
In the left column, disks with a temperature slope of $-1$ are shown, while
the right column represents disks with a temperature slope of $-0.5$. Each
row, from top to bottom, corresponds to shadow amplitudes of $\epsilon=0.5$
and $\epsilon=0.8$. The $x$-axis of the subfigures represents $\beta$, while
the $y$-axis represents $\alpha$. Within each $\beta$-$\alpha$ section, there
are three rows indicating the dominant structures in the disk:
vortices/crescents, rings, and spirals, each represented by different types of
squares, colored by the relevant physical properties as indicated in the color
bars. The figure also includes red and blue line shaded areas, indicating
disks undergoing transitions from vortex-ring and ring-spiral phases,
respectively. We note that the inviscid ($\alpha=0$) simulations maintain the
same temperature gradient with a slope of $-1$ (indicating that the $p$ value
shown in the title of each subfigure only applies to viscid runs) and vary the
density gradient of $d=-0.5$ and $d=-1$ in the left and right columns,
respectively, which help us exclude the influence of density gradient.

Figure 14: Statistics of substructures’ vorticity and density ratio for 160
simulations. The structure of the figure is described in Appendix B and is
similar to Figure 5. Figure 15: Statistics of vortices’ aspect ratio, rings’
eccentricity and spirals’ pitch angles for 160 simulations. The structure of
the figure is same as Fig14.

##  Appendix C One-sided shadow test

In this Appendix, we briefly examine how the morphology and form of
substructures can be affected by the morphology of the shadow region. As an
experiment, we performed simulations with only the right side of the shadow
shown in Figure 1 present, and the target temperature is taken as

$T_{\rm tar}(r,\phi)=T_{\rm init}(r)\left(1-\epsilon
e^{-\frac{\phi^{2}}{2\sigma_{\phi}^{2}}}\right).$ (C1)

The remaining parameters for the disk and shadow are the same as those in the
representative simulations (NL-hm runs). For detailed parameter settings for
the NL-hm runs, please refer to Table 1. It can be seen from Figure 16 that
the types of dominant substructures have not changed compare with NL-hm runs.
The dominant spiral now has $m=1$, and the rings become asymmetric (with
$m=1$, as opposed to eccentric with $m=2$), while crescents are generated as
usual. These outcomes similarly follow the formation process described in
Section 5.2. These simulations illustrate that besides a morphological change
from $m=2$ to $m=1$, the general trends of shadow-driven substructures are not
sensitive to shadow prescriptions.

Figure 16: Normalized density perturbation for one-sided shadow. Left panel:
Spiral forming disk with parameter taken to be same as run NL-hm-S-NR (Figure
2). Middle panel: Ring forming disk with parameter taken to be same as run NL-
hm-R-NR (Figure 3). Right panel: Vortex forming disk with parameter taken to
be same as run NL-hm-V-NR (Figure 4).

## References

  * ALMA Partnership et al. (2015) ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3, doi: 10.1088/2041-8205/808/1/L3
  * Andrews (2020) Andrews, S. M. 2020, ARA&A, 58, 483, doi: 10.1146/annurev-astro-031220-010302
  * Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44, doi: 10.3847/1538-4357/aab846
  * Bae et al. (2023) Bae, J., Isella, A., Zhu, Z., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 423, doi: 10.48550/arXiv.2210.13314
  * Bae et al. (2021) Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56, doi: 10.3847/1538-4357/abe45e
  * Bai (2016) Bai, X.-N. 2016, ApJ, 821, 80, doi: 10.3847/0004-637X/821/2/80
  * Benisty et al. (2015) Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6, doi: 10.1051/0004-6361/201526011
  * Benisty et al. (2018) Benisty, M., Juhász, A., Facchini, S., et al. 2018, A&A, 619, A171, doi: 10.1051/0004-6361/201833913
  * Benisty et al. (2023) Benisty, M., Dominik, C., Follette, K., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 605, doi: 10.48550/arXiv.2203.09991
  * Birnstiel et al. (2013) Birnstiel, T., Dullemond, C. P., & Pinilla, P. 2013, A&A, 550, L8, doi: 10.1051/0004-6361/201220847
  * Casassus et al. (2018) Casassus, S., Avenhaus, H., Pérez, S., et al. 2018, MNRAS, 477, 5104, doi: 10.1093/mnras/sty894
  * de Val-Borro et al. (2006) de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529, doi: 10.1111/j.1365-2966.2006.10488.x
  * Debes et al. (2017) Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205, doi: 10.3847/1538-4357/835/2/205
  * Doi & Kataoka (2021) Doi, K., & Kataoka, A. 2021, ApJ, 912, 164, doi: 10.3847/1538-4357/abe5a6
  * Gammie (2001) Gammie, C. F. 2001, ApJ, 553, 174, doi: 10.1086/320631
  * Garufi et al. (2018) Garufi, A., Benisty, M., Pinilla, P., et al. 2018, A&A, 620, A94, doi: 10.1051/0004-6361/201833872
  * Ginski et al. (2021) Ginski, C., Facchini, S., Huang, J., et al. 2021, ApJ, 908, L25, doi: 10.3847/2041-8213/abdf57
  * Gratton et al. (2019) Gratton, R., Ligi, R., Sissa, E., et al. 2019, A&A, 623, A140, doi: 10.1051/0004-6361/201834760
  * GRAVITY Collaboration et al. (2021) GRAVITY Collaboration, Sanchez-Bermudez, J., Caratti O Garatti, A., et al. 2021, A&A, 654, A97, doi: 10.1051/0004-6361/202039600
  * Huang et al. (2020) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2020, ApJ, 891, 48, doi: 10.3847/1538-4357/ab711e
  * Huang et al. (2018) Huang, P., Isella, A., Li, H., Li, S., & Ji, J. 2018, ApJ, 867, 3, doi: 10.3847/1538-4357/aae317
  * Isella et al. (2018) Isella, A., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L49, doi: 10.3847/2041-8213/aaf747
  * Kuffmeier et al. (2021) Kuffmeier, M., Dullemond, C. P., Reissl, S., & Goicovic, F. G. 2021, A&A, 656, A161, doi: 10.1051/0004-6361/202039614
  * Lin & Shu (1964) Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646, doi: 10.1086/147955
  * Lin & Youdin (2015) Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17, doi: 10.1088/0004-637X/811/1/17
  * Liu et al. (2022) Liu, Y., Flock, M., & Fang, M. 2022, Science China Physics, Mechanics, and Astronomy, 65, 269511, doi: 10.1007/s11433-021-1891-8
  * Lyra & Lin (2013) Lyra, W., & Lin, M.-K. 2013, ApJ, 775, 17, doi: 10.1088/0004-637X/775/1/17
  * Monnier et al. (2017) Monnier, J. D., Harries, T. J., Aarnio, A., et al. 2017, ApJ, 838, 20, doi: 10.3847/1538-4357/aa6248
  * Montesinos & Cuello (2018) Montesinos, M., & Cuello, N. 2018, MNRAS, 475, L35, doi: 10.1093/mnrasl/sly001
  * Montesinos et al. (2016) Montesinos, M., Perez, S., Casassus, S., et al. 2016, ApJ, 823, L8, doi: 10.3847/2041-8205/823/1/L8
  * Muro-Arena et al. (2018) Muro-Arena, G. A., Dominik, C., Waters, L. B. F. M., et al. 2018, A&A, 614, A24, doi: 10.1051/0004-6361/201732299
  * Muro-Arena et al. (2020) Muro-Arena, G. A., Benisty, M., Ginski, C., et al. 2020, A&A, 635, A121, doi: 10.1051/0004-6361/201936509
  * Nealon et al. (2019) Nealon, R., Pinte, C., Alexander, R., Mentiplay, D., & Dipierro, G. 2019, MNRAS, 484, 4951, doi: 10.1093/mnras/stz346
  * Pérez et al. (2018) Pérez, L. M., Benisty, M., Andrews, S. M., et al. 2018, ApJ, 869, L50, doi: 10.3847/2041-8213/aaf745
  * Pfeil & Klahr (2019) Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150, doi: 10.3847/1538-4357/aaf962
  * Pinilla et al. (2015) Pinilla, P., Birnstiel, T., & Walsh, C. 2015, A&A, 580, A105, doi: 10.1051/0004-6361/201425539
  * Pinte et al. (2023) Pinte, C., Teague, R., Flaherty, K., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 645, doi: 10.48550/arXiv.2203.09528
  * Pohl et al. (2017) Pohl, A., Benisty, M., Pinilla, P., et al. 2017, ApJ, 850, 52, doi: 10.3847/1538-4357/aa94c2
  * Rich et al. (2020) Rich, E. A., Wisniewski, J. P., Sitko, M. L., et al. 2020, ApJ, 902, 4, doi: 10.3847/1538-4357/abb2a3
  * Rich et al. (2019) Rich, E. A., Wisniewski, J. P., Currie, T., et al. 2019, ApJ, 875, 38, doi: 10.3847/1538-4357/ab0f3b
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Speedie et al. (2022) Speedie, J., Booth, R. A., & Dong, R. 2022, ApJ, 930, 40, doi: 10.3847/1538-4357/ac5cc0
  * Stolker et al. (2016) Stolker, T., Dominik, C., Avenhaus, H., et al. 2016, A&A, 595, A113, doi: 10.1051/0004-6361/201528039
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4, doi: 10.3847/1538-4365/ab929b
  * Sturm et al. (2020) Sturm, J. A., Rosotti, G. P., & Dominik, C. 2020, A&A, 643, A92, doi: 10.1051/0004-6361/202038919
  * Suzuki et al. (2016) Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74, doi: 10.1051/0004-6361/201628955
  * Tabone et al. (2022) Tabone, B., Rosotti, G. P., Cridland, A. J., Armitage, P. J., & Lodato, G. 2022, MNRAS, 512, 2290, doi: 10.1093/mnras/stab3442
  * Uyama et al. (2018) Uyama, T., Hashimoto, J., Muto, T., et al. 2018, AJ, 156, 63, doi: 10.3847/1538-3881/aacbd1
  * van Boekel et al. (2017) van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132, doi: 10.3847/1538-4357/aa5d68
  * van der Marel et al. (2016) van der Marel, N., Cazzoletti, P., Pinilla, P., & Garufi, A. 2016, ApJ, 832, 178, doi: 10.3847/0004-637X/832/2/178
  * Varga et al. (2021) Varga, J., Hogerheijde, M., van Boekel, R., et al. 2021, A&A, 647, A56, doi: 10.1051/0004-6361/202039400
  * Villenave et al. (2023) Villenave, M., Stapelfeldt, K. R., Duchene, G., et al. 2023, arXiv e-prints, arXiv:2311.07668, doi: 10.48550/arXiv.2311.07668
  * Wolff et al. (2016) Wolff, S. G., Perrin, M., Millar-Blanchaer, M. A., et al. 2016, ApJ, 818, L15, doi: 10.3847/2041-8205/818/1/L15
  * Zhang et al. (2021) Zhang, S., Hu, X., Zhu, Z., & Bae, J. 2021, ApJ, 923, 70, doi: 10.3847/1538-4357/ac2c82
  * Zhu et al. (2014) Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-n. 2014, ApJ, 785, 122, doi: 10.1088/0004-637X/785/2/122