# Boundary Layers of Accretion Disks: Acoustic, Vortex-Driven and Other Modes

Matthew S. B. Coleman,1,2Roman R. Rafikov,1,3 Alexander A. Philippov4
1Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, USA
2Department of Astrophysical Sciences, 4 Ivy Lane, Princeton University,
Princeton, NJ 08540, USA
3Centre for Mathematical Sciences, Department of Applied Mathematics and
Theoretical Physics, University of Cambridge,
Wilberforce Road, Cambridge CB3 0WA, UK
4Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue,
New York, NY 10010, USA Corresponding author<EMAIL_ADDRESS>N.
Bahcall Fellow at the IAS

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Disk accretion onto weakly magnetized objects possessing a material surface
must proceed via the so-called boundary layer (BL) — a region at the inner
edge of the disk, in which the velocity of accreting material abruptly
decreases from its Keplerian value. Supersonic shear arising in the BL is
known to be conducive to excitation of acoustic waves that propagate into the
accretor and the disk, enabling angular momentum and mass transport across the
BL. We carry out a numerical exploration of different wave modes that operate
near the BL, focusing on morphological characteristics of the modes in the
innermost parts of accretion disk. Using a large suite of simulations covering
a broad range of Mach numbers (of the supersonic shear flow in the BL), we
provide accurate characterization of the different types of modes, verifying
their properties against analytical results, when available. We discover new
types of modes, in particular, global spiral density waves launched by
vortices forming in the disk near the BL as a result of the Rossby wave
instability; this instability is triggered by the vortensity production in
that region caused by the nonlinear damping of acoustic waves. Azimuthal
wavenumbers of the dominant modes that we observe appear to increase
monotonically with the Mach number of the runs, but a particular mix of modes
found in a simulation is mildly stochastic. Our results provide a basis for
better understanding of the angular momentum and mass transport across the BL
as well as the emission variability in accreting objects.

###### keywords:

accretion, accretion discs – hydrodynamics – instabilities

††pubyear: 2015††pagerange: Boundary Layers of Accretion Disks: Acoustic,
Vortex-Driven and Other Modes–B.3

##  1 Introduction

Accretion disks are ubiquitous in astrophysics, with objects ranging from
active galactic nuclei to protostars being fundamentally tied to them. In the
cases where the central object (i.e. accretor) is not a black hole, but is a
neutron star, a white dwarf, a protostar, or a protoplanet (henceforth we
refer to any of these objects as a “star"), the accretor has a material
surface, which the accreted material must connect to in some fashion. If the
accretion rate is high and the magnetic field of the star is sufficiently low,
then the accretion flow does not get disrupted by magnetic stresses (Ghosh et
al., 1977; Koenigl, 1991) and the disk can extend all the way to the surface
of the star. This particular situation inevitably requires accreting gas to
transition from rapid, supersonic rotation (at Keplerian velocity) in the disk
to a slow rotation in the star. The region of the disk-star system where this
transition takes place is known as the boundary layer (BL). Systems where the
BLs are expected to emerge include e.g. FU Ori type young stellar objects
(Popham et al., 1993) and cataclysmic variables (CVs, Kippenhahn & Thomas
1978; Narayan & Popham 1993). Weakly magnetized neutron stars in low-mass
X-ray binaries (LMXBs) are believed to accrete in a broadly similar fashion
through the so-called spreading layer (Inogamov & Sunyaev, 1999, 2010;
Gilfanov et al., 2003; Revnivtsev & Gilfanov, 2006; Philippov et al., 2016).
Objects accreting gas directly onto their surfaces through the BLs are the
subject of the present paper.

In order for the material arriving from the Keplerian disk to become a part of
a slowly rotating star it must somehow lose its angular momentum. While the
magnetorotational instability (MRI; Velikhov 1959; Chandrasekhar 1960; Balbus
& Hawley 1991) is traditionally invoked as the favored angular momentum
transport mechanism in ionized Keplerian accretion disks, it would not operate
in the BL. This is because the MRI requires that the angular frequency
$\Omega$ of the fluid flow decays with the distance,
$\mathrm{d}\Omega/\mathrm{d}r<0$, whereas the BL naturally has
$\mathrm{d}\Omega/\mathrm{d}r>0$, preventing the MRI from operating (Pessah &
Chan, 2012). This conclusion has been verified by MHD simulations of the BLs
(Belyaev et al., 2013b).

Instead the disk must utilize a different mechanism to remove angular momentum
from the accreting gas, which passes through the BL on its way to the surface
of the star. Belyaev & Rafikov (2012) identified a robust mechanism for doing
that — a linear instability operating in a supersonic shear flow, which
generates acoustic waves in the BL where the azimuthal velocity of the flow
exhibits sharp supersonic variation. These waves propagate both out in the
disk and into the star, allowing energy and angular momentum of accreting gas
to be transported over significant distances before being dissipated.
Numerical simulations later confirmed that this angular momentum transport
mechanism robustly operates within the BL, both in hydrodynamic (Belyaev et
al., 2012, 2013a; Hertfelder & Kley, 2015a) and magnetohydrodynamic (Belyaev
et al., 2013b; Belyaev & Quataert, 2017) settings.

This discovery marked a significant paradigm shift compared to the local
transport mechanisms invoked in previous studies of the BL problem (Kippenhahn
& Thomas, 1978; Popham et al., 1993; Narayan & Popham, 1993; Hertfelder et
al., 2013). The intrinsically non-local nature of this mechanism could
substantially impact disk thermodynamics and the shape of the accretion disk
spectrum (Belyaev et al., 2012, 2013a). Another important implication follows
from the fact that the modes excited in the BL are intrinsically non-
axisymmetric. This should lead to the variability of emission produced in the
near-BL part of the disk and may explain the various types of quasi-periodic
variability observed in objects accreting through the BLs (e.g. CVs, see e.g.
Warner 2003).

These ramifications, as well as the ubiquity of the BLs in astrophysics,
motivate further efforts to better understand the physics of the BLs through
numerical simulations, building on the previous work of Belyaev et al. (2012,
2013a, 2013b). These past studies, while significantly advancing our
understanding of the BL structure, were often limited in terms of their
numerical resolution, duration of the simulations, and the number of model
parameters that have been varied.

In this paper, first in a series, we present a new set of long-term, high-
resolution, hydro simulations focused on exploring the BL physics. We provide
extensive exploration of both the physical and numerical parameter space to
test the sensitivity of outcomes to both types of simulation inputs. After
describing our simulation suite, we carry out an in-depth analysis of the mode
structure of the perturbations that arise in the vicinity of the BL as a
result of ongoing acoustic instability. One of the highlights of this study is
the discovery of new types of modes, naturally emerging in this disk region,
and our attempts at understanding their origin. In the future we will use this
numerical data set to analyze angular momentum and mass transport driven by
the different modes operating in the vicinity of the BL (Coleman et al. in
prep).

Our paper is organized as follows. In §2 we discuss our physical setup and
typical values of the Mach number in different astrophysical objects to
motivate the parameter choices for our simulations. We remind the basics of
the acoustic mode phenomenology in §3, and cover the details of our numerical
setup in §4. We provide detailed morphological description of the various
modes that we find in our runs for different values of the Mach number in §5.
Description of the vortex-driven modes and explanation of their origin are
provided in §6 & 7, respectively. We discuss properties of other modes found
in our simulations, as well as some other aspects of our work, in §7, and
summarize our main findings in §9.

##  2 Physical Setup and Typical Mach Numbers

In this work we consider a system consisting of a central object with a
surface (a star) and an accretion disk extending all the way to the star, i.e.
having a physical contact with its surface. We study the evolution of this
system in two-dimensional (2D, vertically integrated), hydrodynamic (i.e. no
magnetic fields) setup. The disk is non-self-gravitating and orbits in a
Newtonian potential of a central point mass $M_{\star}$. Very importantly, the
disk has no intrinsic viscosity111Our simulations have no explicit viscosity
and numerical viscosity is negligible. so that any mass re-distribution
(accretion) in the system can take place only due to the action of the waves
propagating in the disk and the star.

Similar to a number of previous studies of the BL (Belyaev et al., 2012,
2013a, 2013b), we treat disk thermodynamics using the globally isothermal
equation of state (EoS), $P=\Sigma c_{s}^{2}$, where $P$ is the vertically
integrated pressure, $\Sigma$ is the surface density, and $c_{\rm s}$ is the
sound speed which is constant. The advantage of using this simple EoS comes
from the fact that it keeps the disk thermodynamics unchanged through its
evolution, while at the same time naturally allowing conservation of the
angular momentum flux of the waves propagating in the disk (Miranda & Rafikov,
2019b, 2020). The latter point is very important for us since we are focusing
on the wave-driven evolution of the system. Another commonly used EoS, locally
isothermal, has been demonstrated by Miranda & Rafikov (2019b, 2020) to not
conserve the angular momentum of the waves, a phenomenon that would have
greatly complicated understanding of the wave-driven transport (Coleman et al.
in prep). On the other hand, an adiabatic EoS with $P\propto\Sigma^{\gamma}$,
$\gamma\neq 1$ would lead to evolution of the thermal state of the disk as a
result of entropy production at the shocks inevitably arising in the system
(see below), again complicating the interpretation of the results. Implicit in
our choice of the EoS is the assumption of gas pressure dominating the total
pressure, i.e. radiation pressure is neglected. Therefore, results of our
study are not directly applicable to accreting neutron stars, as radiation
pressure plays an important role in such systems.

Another reason for using the globally isothermal EoS is that it allows us to
reduce the number of parameters needed to characterize the system. With this
setup all details of the disk thermodynamics are captured in a single variable
— constant gas sound speed $c_{\rm s}$. As a result, the key physical
parameter governing the behavior of the system in our runs is the
dimensionless Mach number ${\mathcal{M}}$ defined as the ratio of the
Keplerian speed at the surface of the star ($r=R_{\star}$) to the sound speed
$c_{s}$:

$\displaystyle\mathcal{M}$
$\displaystyle\equiv\frac{\Omega_{K}(R_{\star})R_{\star}}{c_{\rm s}},$ (1)

where $\Omega_{K}=(GM_{\star}/r^{3})^{1/2}$ is the Keplerian rotation rate.

To provide motivation for the values of ${\mathcal{M}}$ explored in this work,
we estimate characteristic Mach numbers in some astrophysical systems that may
naturally host BLs. We do this by relating the midplane temperature ($T_{\rm
m}$) of an optically thick accretion disk to its optical depth ($\tau>1$),
accretion rate ($\dot{M}$), and local orbital frequency ($\Omega$) in a
standard fashion (Shakura & Sunyaev, 1973):

$\displaystyle T_{\rm m}^{4}$
$\displaystyle=\dfrac{3}{8\pi}\dfrac{\dot{M}\Omega^{2}}{\sigma_{\rm SB}}\tau,$
(2)

where $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant. Using this value of
the temperature to compute the characteristic sound speed in the disk via
$c_{\rm s}^{2}=k_{\rm B}T_{\rm m}/\mu$ (where $k_{\rm B}$ is the Boltzmann
constant and $\mu$ is the mean molecular weight) and assuming for simplicity
that Eqn. (2) holds all the way to the surface of the star (i.e. down to
$R=R_{\star}$), we obtain the following general expression for the
characteristic Mach number:

$\displaystyle\mathcal{M}=\left(\dfrac{8\pi}{3}\dfrac{G^{3}\sigma}{k_{\rm
B}^{4}}\right)^{1/8}M_{\star}^{3/8}\mu^{1/2}\dot{M}^{-1/8}R_{\star}^{-1/8}\tau^{-1/8}.$
(3)

One obvious class of astrophysical objects that may feature the BLs are the
accreting white dwarfs, CVs and AM CVn systems. Because of thermal instability
in the disk, these systems can exist in two states characterized by high or
low values of $\dot{M}$. In the high-$\dot{M}$ state one finds for the typical
parameters of CVs

$\displaystyle\mathcal{M}=32\left(\dfrac{M_{\star}}{0.6M_{\odot}}\right)^{3/8}\left(\dfrac{\mu}{0.6}\right)^{1/2}\left(\dot{M}_{-9}\,\dfrac{R_{\star}}{1.4R_{\earth}}\,\dfrac{\tau}{10^{4}}\right)^{-1/8},$
(4)

where $\dot{M_{i}}=\dot{M}/10^{i}M_{\odot}$ yr-1, while for the AM CVn
systems222AM CVns typically have relatively high accretor masses (see e.g.
Roelofs et al., 2007). This causes them to have small radii due to the
degenerate nature of the mass-radius relation for white dwarfs.

$\displaystyle\mathcal{M}=50\left(\dfrac{M_{\star}}{0.9M_{\odot}}\right)^{3/8}\left(\dfrac{\mu}{1.4}\right)^{1/2}\left(\dot{M}_{-9}\,\dfrac{R_{\star}}{0.8R_{\earth}}\,\dfrac{\tau}{5\\!\times\\!10^{4}}\right)^{-1/8}.$
(5)

The optical depth $\tau$ is estimated from stratified shearing-box simulations
in the respective regime333Note there is a factor of 2 difference here as
these papers define $\tau_{\rm tot}$ as twice the midplane $\tau$.: Hirose et
al. (2014) and Coleman et al. (2016) for CVs and Coleman et al. (2018) for AM
CVns. In the low-$\dot{M}$ states of these systems (when $\dot{M}$ drops by
2-3 orders of magnitude) these estimates would go up to ${\mathcal{M}}\sim
300$, however, it is not even clear if the disk would not be disrupted by even
the weak stellar magnetic field when $\dot{M}$ is so low.

For FU Ori stars episodically accreting at high $\dot{M}$ from the
protoplanetary disk we find

$\displaystyle\mathcal{M}=3.9\left(\dfrac{M}{M_{\odot}}\right)^{3/8}\left(\dfrac{\mu}{0.6}\right)^{1/2}\left(\dot{M}_{-5}\,\dfrac{R_{\star}}{2R_{\odot}}\,\dfrac{\tau}{6\times
10^{5}}\right)^{-1/8},$ (6)

where $\tau$ is taken from the simulations of Hirose (2015).

Motivated by these estimates, and taking into account the numerical
constraints that would permit an efficient parameter space exploration, in
this work we focus on exploring the BL physics for the values of
${\mathcal{M}}$ in the range $5\leq{\mathcal{M}}\leq 15$. We run at least one
simulation for each integer value of ${\mathcal{M}}$ in this range, although
for several characteristic values of ${\mathcal{M}}=6,9,12$ we run multiple
simulations to explore the sensitivity to the initial conditions, resolution,
etc.

##  3 Acoustic mode phenomenology

We now remind the reader some facts about the acoustic modes excited in the
BLs, that have been previously discussed in Belyaev & Rafikov (2012) and
Belyaev et al. (2012, 2013a). Here we simply summarize the main points made in
these studies, which will facilitate the description of the modes found in
this work.

Highly supersonic shear present in the BL efficiently drives the non-
axisymmetric acoustic (sonic) waves propagating on both sides of the shear
layer, in which the azimuthal velocity drop takes place (Belyaev & Rafikov,
2012). These waves are global and propagate both inside the star and outside
of it, i.e. in the disk. In general, there are three different types of modes
that can get excited in the system, but the simulations typically exhibit only
two of them, termed lower and upper modes in Belyaev et al. (2013a). These
modes have quite distinct appearance both in the disk and inside the star, and
obey very different dispersion relations. They are described in more detail
next and are illustrated in Figure 1 showing the 2D snapshots of
$rv_{r}\sqrt{\Sigma}$ — a quantity that should be conserved for a sound wave
propagating through the disk. We will routinely display the spatial maps of
$rv_{r}\sqrt{\Sigma}$ in the $r-\phi$ coordinate plane (with $\phi$ as the
vertical axis), to highlight the details of the morphological features of
acoustic modes near $r=R_{\star}$ where they are excited.

###  3.1 Upper modes

The defining feature of the upper modes is that inside the star they have
$k_{r}\neq 0$, i.e. the wave crests of the perturbation pattern associated
with this mode are substantially inclined with respect to the radial
direction. In the disk this mode starts off with $k_{r}=0$ as $r\to R_{\star}$
from above, but further out in the disk the perturbation pattern gets wrapped
up by the differential rotation into multiple trailing spiral arms, see
Section 3.3. This pattern is illustrated in Figure 1a. Note the sign change of
the perturbation variable ($v_{r}$ in this case) across the BL.

Figure 1:  Example snapshots from two simulations showing the quantity
$rv_{r}\sqrt{\Sigma}$ (a proxy for acoustic wave action) as a function of $r$
and $\phi$ in Cartesian coordinates. The white dotted lines at $r=1$ separate
the star from the disk. (a) A snapshot from an ${\mathcal{M}}=6$ run at 25
inner orbits (during the linear growth phase of the sonic instability) where
an $m=7$, $\Omega_{\rm p}=0.758$ global upper mode is clearly visible in the
disk and in the star. (b) A snapshot of an ${\mathcal{M}}=9$ run at 175 inner
orbits. The dominant global lower mode with $m=19$ and $\Omega_{\rm p}=0.315$
is clearly visible inside the star and is trapped in the inner disk,
$r\lesssim 2$ (Fig. 3 shows $m=19$ as having the most power at $r=0.92,\,1.20$
at this time). The dot-dashed white lines near $r\approx 2.2$ are the inner
and outer Lindblad resonances for the $m=19$ mode, while the solid line
indicates the corotation radius. The criss-cross structure results from the
outward-propagating $m=19$ mode reflecting off the Inner Lindblad Resonance
and self-intersecting. In both panels the dashed white curves show the
expected shape of the wave pattern for each of the modes, computed using the
WKB equations (10)-(11) and the measured values of $m$ and $\Omega_{\rm p}$.

Belyaev et al. (2013a) came up with the following dispersion relation between
the azimuthal wavenumber $m$ and the pattern speed $\Omega_{\rm p}$ for the
upper modes, which should hold approximately for $m\gg 1$:

$\displaystyle\left[\Omega(R_{\star})-\Omega_{\rm
p}\right]^{2}=\frac{c_{s}^{2}}{R_{\star}^{2}}+\frac{\kappa^{2}(R_{\star})}{m^{2}},$
(7)

where $\Omega(R_{\star})$ and $\kappa(R_{\star})$ are the values of the disk
angular and epicyclic frequencies as $r\to R_{\star}$. Note that in our runs
we often find deviations of $\Omega(r)$ from the Keplerian frequency
$\Omega_{K}$ in the disk near the star (but outside the BL). These deviations
are caused by the non-trivial contribution of pressure support to the radial
momentum balance

$\displaystyle\Omega^{2}(r)=\Omega_{K}^{2}(r)+\frac{1}{\Sigma r}\frac{\partial
P}{\partial r},$ (8)

enabled by the restructuring of the disk surface density near the star (see
§5.1). For that reason, epicyclic frequency $\kappa(r)$ in the disk near the
star is generally not equal to $\Omega_{K}(r)$, as it would be in a purely
Keplerian flow.

In this work we analyze our simulation outputs using a relation between $m$
and $\Omega_{\rm p}$, which is more accurate than the equation (7). The
derivation of this refined dispersion relation (23) can be found in Appendix A
and it is plotted for many values of ${\mathcal{M}}$ in Figure 14, clearly
showing that $\Omega_{\rm p}$ of the upper modes increases with increasing
$m$.

###  3.2 Lower modes

The lower modes exhibit the $k_{r}$ behavior, which is opposite to that of the
upper modes: inside the star they propagate in the azimuthal direction only,
i.e. their radial wavenumber inside the star vanishes, $k_{r}=0$. At the same
time, just outside the star the lower modes have $k_{r}\neq 0$, and further in
the disk $k_{r}$ gets modified by the differential rotation. This is
illustrated in Figure 1b, where inside the star, for $r<R_{\star}$, one can
see the perturbation pattern perfectly aligned with the radial direction,
whereas just outside the star the wave crests are inclined with respect to the
radial direction already at $r=R_{\star}$.

The dispersion relation for the lower modes was derived in Belyaev et al.
(2013a) who showed that the pattern speed $\Omega_{\rm p}$ of the lower mode
is related to its azimuthal wavenumber $m$ through the following relation

$\displaystyle\dfrac{\Omega_{\rm
p}}{\Omega\left(R_{\star}\right)}=\dfrac{R_{\star}}{r_{0}}\sqrt{\mathcal{M}^{-2}+\left(\dfrac{\mathcal{M}R_{\star}}{2mr_{0}}\right)^{2}},$
(9)

where $r_{0}$ is a parameter, which is close to unity, see Belyaev et al.
(2013a). Figure 14 illustrates this dispersion relation for a number of values
of ${\mathcal{M}}$, demonstrating that $\Omega_{\rm p}$ of the lower modes
decreases with increasing $m$.

###  3.3 Propagation of the BL-excited acoustic modes in the disk

While the appearance of the upper and lower modes just outside the BL is very
different (in terms of their $k_{r}$), their subsequent propagation in the
disk follows the classical behavior of the density waves in differentially
rotating disks (Binney & Tremaine, 2008). In particular, linear density wave
theory for a Keplerian disk predicts that an $m$-th azimuthal harmonic of a
perturbed fluid variable $f$ behaving as $f\sim\exp(i\Phi_{m})$, where
$\Phi_{m}=\int^{r}k_{r}(r^{\prime})dr^{\prime}+m\varphi$ and $k_{r}$ is the
radial wavenumber, obeys the standard WKB dispersion relation (Goldreich &
Tremaine, 1980):

$\displaystyle m^{2}\left[\Omega_{K}(r)-\Omega_{\rm
p}\right]^{2}=\Omega_{K}^{2}(r)+c_{s}^{2}k_{r}^{2}.$ (10)

Wave crests trace the trajectory along which perturbation maintains a constant
phase $\Phi_{m}$, so that $d\Phi_{m}=0$. Using the expression for $\Phi_{m}$,
one finds

$\displaystyle\frac{dr}{d\varphi}=-\frac{m}{k_{r}(r)}.$ (11)

Integrating this equation with the $k_{r}(r)$ determined by the equation (10)
one obtains the shape (i.e. the $\varphi(r)$ dependence) of the wave crests of
the $m$-th harmonic of the fluid perturbation with the pattern speed
$\Omega_{\rm p}$.

In Figure 1a the white dashed curve shows the analytical prediction for the
wave crest location computed using the equations (10)-(11) and the values of
$\Omega_{\rm p}$ and $m$ derived for the upper mode dominating in this
snapshot (see captions). One can see the analytical calculation agreeing with
the actual wave crest shape extremely well and predicting the mode to form a
pattern of spiral arms sheared by the differential rotation and propagating
out to large distances (this mode has a narrow evanescent region near the star
where $\Omega(r)>\Omega_{\rm p}$, in which $k_{r}=0$).

The same calculation done for the values of $\Omega_{\rm p}$ and $m$
characterizing the lower mode dominating in Figure 1b, are shown by the white
dashed curve in that panel and are seen to also agree very well with the shape
of the outgoing and the incoming sonic waves. One can see that the lower modes
are trapped in a resonant cavity extending between the stellar surface and the
Inner Lindblad Resonance (ILR), at which $k_{r}=0$ for that mode. For the
modes discussed in this work the radial location $r_{\mathrm{ILR}}$ of the ILR
typically ends up far enough in the disk where $\Omega(r)$ can be well
approximated by $\Omega_{K}$. Then $r_{\mathrm{ILR}}$ is determined by the
condition $\Omega_{K}(r_{\mathrm{ILR}})(1-m^{-1})=\Omega_{\rm p}$, given the
mode’s $m$ and $\Omega_{\rm p}$, so that

$\displaystyle
r_{\mathrm{ILR}}=R_{\star}\left[\frac{\Omega_{K}(R_{\star})}{\Omega_{\rm
p}}\frac{m-1}{m}\right]^{2/3}.$ (12)

For $m\gg 1$ the location of ILR is close to the corotation radius
$r_{\mathrm{c}}=R_{\star}\left[\Omega_{K}(R_{\star})/\Omega_{\rm
p}\right]^{2/3}$, see Figure 1b.

##  4 Numerical setup

We simulate the BL and its vicinity — outer layers of a star and inner regions
of an accretion disk — in (vertically integrated) cylindrical geometry, using
Athena++ (Stone et al., 2020) to solve the hydrodynamic equations

$\displaystyle\dfrac{\partial\Sigma}{\partial
t}+\nabla\cdot\left(\Sigma\mathbf{v}\right)$ $\displaystyle=0,$ (13)
$\displaystyle\dfrac{\partial\Sigma\mathbf{v}}{\partial
t}+\nabla\cdot\left[\Sigma\mathbf{vv}+\Sigma c_{s}^{2}\mathbf{I}\right]$
$\displaystyle=-\Sigma\mathbf{\nabla}\Phi,$ (14)

where $\mathbf{v}$ is the fluid velocity, $\mathbf{I}$ is the identity tensor,
and $\Phi$ is the gravitational point mass potential for the star. For all
runs we used the HLLE Riemann solver, second-order van Leer time integrator,
and second-order piecewise-linear primitive reconstruction.

Our simulation domain extends from $r_{\min}$ to $r_{\rm max}$ in the radial
direction and covers full $2\pi$ in azimuthal direction ($\phi$). Our grid is
uniformly spaced in $\phi$ and logarithmically spaced in $r$ (i.e. $\delta
r\propto r$). We choose $r_{\rm max}=4R_{\star}$, far enough in the disk to
ensure that the structures emerging in our runs can fit within the simulation
domain (e.g. see §6.4). The inner boundary inside the star is placed at
$r_{\min}$ such that the density contrast $\Sigma\left(r_{\rm
min}\right)/\Sigma\left(R_{\star}\right)=10^{7}$, where we assume a isothermal
hydrostatic atmospheric profile inside the star. This choice of $r_{\rm min}$
was made to minimize the simulations dependence of $r_{\rm min}$; see Section
4.3 for more details.

We carry out a detailed analysis of our simulation outputs, both on-the-fly
and in post-processing. In particular, we analyze the behavior of the fluid
variables in Fourier domain and develop a fully automated procedure for
detecting and measuring properties of the various wave-like perturbations
present in our simulations. These and other analysis modules are described in
more detail in §B. The ability to not only infer the existence of multiple
modes and derive mode wavenumber $m$ and pattern speed $\Omega_{\rm p}$, but
to also follow their evolution throughout the full duration of a simulation is
what makes our analysis extremely powerful. It enables us to see important
trends and patterns across dozens of the BL simulations performed for
different values of ${\mathcal{M}}$. Some other details of our numerical setup
are described below.

###  4.1 Units

To define simulation units we took the surface of the star as our fiducial
location. We chose $r=1$ and $\Sigma=1$ to correspond to this location (at
$t=0$), and for the Keplerian velocity at the surface of the star
$v_{K}(R_{\star})$ to be unity, making $c_{s}=\mathcal{M}^{-1}$ and
$GM_{\star}=1$. This choice makes the Keplerian period at the surface of the
star $\tau_{\star}=2\pi$; this is why we often state times in the form of
$t/2\pi$.

###  4.2 Initial Conditions

To create the initial conditions of our simulations we partitioned the
simulation domain into three regions: star ($r\leq 1-\delta/2$), transition
($1-\delta/2<r<1+\delta/2$), disk ($r\geq 1+\delta/2)$, with $\delta=0.05$. We
initialized the star in hydrostatic equilibrium (HSE),
$\Sigma=\Sigma_{0}\exp\left(\mathcal{M}^{2}/r\right)$, $\mathbf{v}=0$. The
disk is initialized with $\Sigma=r^{-3/2}$,
$\mathbf{v}=r^{-1/2}\hat{\mathbf{\phi}}$, i.e. a pure Keplerian disk
neglecting pressure support (in the beginning of a simulation the disk quickly
adjusts to an equilibrium state accounting for the radial pressure support).
The inner and outer (radial) boundary conditions maintain these initial
conditions, for the star and disk respectively. The transition region smoothly
transition from the star and disk such that $\Sigma$ and $\mathbf{v}$ as well
as their first and second derivatives (with respect to $r$) are continuous.

In order to seed the acoustic instability we add random seed perturbations in
the velocity field (in the disk region only). Let $R\in[0,1)$ be a random
number picked from a uniform probability distribution function (PDF), and
$A=10^{-2}$ be the amplitude of the initial perturbation. To examine any
possible dependence associated with this choice we tried four different
implementations of random seeds:

  1. 1.

Block-random: each block has the same series of random numbers $R$ (one per
cell) and the initial perturbation is $v_{r}=AR$.

  2. 2.

Block-phase-random: similar to Block-random but

$\displaystyle v_{r}(r,\phi)=\dfrac{A}{4}\sum_{n=2,3,5,7}\sin(n\phi+2\pi R).$
(15)

  3. 3.

Globally random: each block receives a different series of random numbers with
$v_{r}=AR$.

  4. 4.

Prime modes: only four random numbers are chosen ($R_{n}$) and

$\displaystyle v_{r}(r,\phi)=\dfrac{A}{4}\sum_{n=2,3,5,7}\sin(n\phi+2\pi
R_{n}).$ (16)

The seed type of each of our simulations can be found in Table 1. The final
letter in each simulation name indicates the seed used for the random number
generator.

###  4.3 Numerical Robustness and Convergence

We experimented with varying several numerical parameters to ensure that our
simulations are insensitive to these choices. To test convergence of our
results (see §8.4), for $\mathcal{M}=6,9,12$ we tried doubling and halving the
resolution. For a given value of Mach number ${\mathcal{M}}$ our fiducial
resolution is typically chosen so as to keep the radial grid scale relative to
the disk scale height $\delta r/h$ roughly constant, with $3\%\lesssim\delta
r/h\lesssim 6\%$. For $\mathcal{M}=9$ we also tried varying the aspect ratio
of the grid cells $N_{r}/N_{\phi}=7/8$ and $2$ and found no noticeable
differences with the fiducial choice of $N_{r}/N_{\phi}=1$.

The numerical parameter that impacts our results most is the inner radial
extent of the simulation, $r_{\rm min}$. As long as $r_{\rm min}$ is such that
$\Sigma\left(r_{\rm min}\right)/\Sigma\left(R_{\star}\right)\approx 10^{7\pm
1}$, we found that there are no substantial differences in the results or runs
with different $r_{\rm min}$. However, at density contrasts $\gtrsim 10^{8}$
tiny numeric fluctuations near $r_{\rm min}$ (likely caused by disagreement
between analytic and numerical hydrostatic equilibrium) would get amplified by
the density contrast and become the dominant source of waves in the
simulation, resulting in large radial oscillations of the star. On the other
hand, below density contrasts of $\sim 10^{5}$ (at $r_{\rm min}$) the acoustic
waves generated in the BL and traveling into the star get prematurely
truncated by the edge of the simulation domain.

###  4.4 Simulation Improvements

While the Athena++ simulations we run are similar to the Athena simulations
presented in Belyaev et al. (2012, 2013a), ours differ in a few key ways.
First, all of our runs extend the full $2\pi$ in the azimuthal direction while
only a handful of the previous simulations cover this angular extent. This is
important for properly capturing all non-axisymmetric structures emerging in
our simulations. Second, our simulations cover a larger radial extent, going
out to a maximum radius of $4R_{\star}$ compared to $2.5R_{\star}$, as before.
We found that this allows us to observe structures that have not been reported
previously, see e.g. §6.4. Third, while the old simulations used a uniform
radial grid, we use a logarithmic radial grid giving us higher effective
resolution at the surface of the star, where the acoustic instability occurs,
for the same number of cells444For our fiducial resolution $\mathcal{M}=9$
simulations this gives us $\sim 2$ times the radial resolution at $R_{\star}$
compared to a uniform radial grid with the same number of cells.. Forth, we
perform Fourier analysis of the outputs by running fast-Fourier transforms
(FFTs) on the fly (see Section B.1) giving us new, previously unattainable,
diagnostic capabilities. This allows us to identify a statistically
significant sample of modes emerging in our simulations, and to discover and
quantify new types of modes, as will be discussed in the following sections.

Name | $\mathcal{M}$ | $R_{\rm min}$ | $N_{R}$ | $N_{\phi}$ | Seed type
---|---|---|---|---|---
M05.FR.mix.a | 5 | 0.608 | 1024 | 1024 | block-phased-mixed
M05.FR.r.a | 5 | 0.608 | 1024 | 1024 | block-random
M06.FR.mix.a | 6 | 0.691 | 1024 | 1024 | block-phased-mixed
M06.FR.prime.a | 6 | 0.691 | 1024 | 1024 | prime modes
M06.FR.r.a | 6 | 0.691 | 1024 | 1024 | block-random
M06.FR.random.a | 6 | 0.691 | 1024 | 1024 | globally random
M06.HR.mix.a | 6 | 0.691 | 2048 | 2048 | block-phased-mixed
M06.HR.mix.b | 6 | 0.691 | 2048 | 2048 | block-phased-mixed
M06.HR.prime.a | 6 | 0.691 | 2048 | 2048 | prime modes
M06.HR.r.a | 6 | 0.691 | 2048 | 2048 | block-random
M06.HR.r.b | 6 | 0.691 | 2048 | 2048 | block-random
M06.HR.r.c | 6 | 0.691 | 2048 | 2048 | block-random
M07.FR.r.a | 7 | 0.752 | 2048 | 2048 | block-random
M08.FR.r.a | 8 | 0.799 | 2048 | 2048 | block-random
M09.LR.mix.a | 9 | 0.834 | 2048 | 2048 | block-phased-mixed
M09.LR.mix.b | 9 | 0.834 | 2048 | 2048 | block-phased-mixed
M09.LR.prime.a | 9 | 0.834 | 2048 | 2048 | prime modes
M09.LR.r.a | 9 | 0.834 | 2048 | 2048 | block-random
M09.LR.r.b | 9 | 0.834 | 2048 | 2048 | block-random
M09.LR.r.c | 9 | 0.834 | 2048 | 2048 | block-random
M09.LR.random.a | 9 | 0.834 | 2048 | 2048 | globally random
M09.FR.mix.a | 9 | 0.834 | 4096 | 4096 | block-phased-mixed
M09.FR.mix.b | 9 | 0.834 | 4096 | 4096 | block-phased-mixed
M09.FR.prime.a | 9 | 0.834 | 4096 | 4096 | prime modes
M09.FR.r.a | 9 | 0.834 | 4096 | 4096 | block-random
M09.FR.r.b | 9 | 0.834 | 4096 | 4096 | block-random
M09.FR.r.c | 9 | 0.834 | 4096 | 4096 | block-random
M09.FR.random.a | 9 | 0.834 | 4096 | 4096 | globally random
M09.HR.r.a | 9 | 0.834 | 8192 | 8192 | block-random
M10.FR.r.a | 10 | 0.861 | 4096 | 4096 | block-random
M11.FR.r.a | 11 | 0.882 | 4096 | 4096 | block-random
M12.FR.mix.a | 12 | 0.899 | 4096 | 4096 | block-phased-mixed
M12.FR.prime.a | 12 | 0.899 | 4096 | 4096 | prime modes
M12.FR.r.a | 12 | 0.899 | 4096 | 4096 | block-random
M12.FR.r.b | 12 | 0.899 | 4096 | 4096 | block-random
M12.FR.random.a | 12 | 0.899 | 4096 | 4096 | globally random
M13.FR.r.a | 13 | 0.913 | 8192 | 8192 | block-random
M14.FR.r.a | 14 | 0.924 | 8192 | 8192 | block-random
M15.FR.r.a | 15 | 0.933 | 8192 | 8192 | block-random
Table 1: See Section 4.2 for descriptions of seed type. Figure 2: Maps of
(left) the vortensity perturbation relative to its initial value, (center) the
perturbation measure $rv_{r}\sqrt{\Sigma}$ in the $r-\phi$ coordinate plane
for a range of $r$, and (right) $rv_{r}\sqrt{\Sigma}$ in physical space
($x(r,\phi),y(r,\phi)$) in the full simulation domain for ${\mathcal{M}}=9$
simulation M09.FR.r.a. Different panels correspond to different moments of
time labeled in each row. Color bars on top show the scale of the vortensity
perturbation and $rv_{r}\sqrt{\Sigma}$ for this run. See text in §5.1 for
detailed description of the evolution shown in this figure. Figure 3:  Time
evolution of the Fourier amplitudes of $v_{r}\sqrt{\Sigma}$ (top panels) and
pattern speeds $\Omega_{\rm p}$ (bottom panels) of a subset of the most
prominent modes in the disk-star system for the ${\mathcal{M}}=9$ run
M.09.FR.r.a. Different columns illustrate the mode amplitude and $\Omega_{\rm
p}$ at different radii: (left) $r=0.92$, inside the star, and (right) at
$r=1.2$, in the inner disk. Different curves are color coded according to the
azimuthal wave number of the mode that they represent, labeled in the inset.
The large ticks on the lower most horizontal axis indicate the temporal
midpoint of an automatically detected global mode (see Section B.1). If any of
these ticks corresponds to one of the plotted modes, then they are drawn in
the same color (e.g. the orange tick at $t/2\pi\approx 225$). Figure 4: Time
evolution of the azimuthally-averaged (a) angular frequency $\Omega(r)$ and
(b) surface density $\Sigma(r)$ for the ${\mathcal{M}}=9$ run M09.FR-r-a
described in Section 5.1. The colors indicated by the color bar denote the
time in inner orbital periods (i.e. $t/2\pi$). See text for details.

##  5 Morphological characterization of various modes for different values of
${\mathcal{M}}$.

We now provide a systematic description of the modes that emerge in
simulations with different ${\mathcal{M}}$ and their evolution during the
runs. Readers not interested in such details can skip this section.

We use several types of diagnostics to illustrate our observations. To
highlight the development, saturation and subsequent evolution of the wave
modes we use the 2D snapshots of $rv_{r}\sqrt{\Sigma}$ in physical coordinates
($x(r,\phi),y(r,\phi)$) at different moments of time starting from the linear
development of the instability until the end of the simulation, see the right
columns of Figures 2, 5, 7, 9. Also, to highlight the perturbation details
near $r=1$, we supplement these maps of $rv_{r}\sqrt{\Sigma}$ with their
projections onto $r-\phi$ coordinate plane, see the central columns of the
same figures. Left columns of these figures illustrate the evolution of the
flow vortensity and are discussed in §6.

Harmonic content of the perturbation pattern is illustrated in Figures 3, 6,
8, 10 for ${\mathcal{M}}=9,6,12$ and $15$, correspondingly. There we show the
amplitudes $A_{m}$ and pattern speeds $\Omega_{\rm p}$ of the dominant modes
identified by our automated mode detection procedure (see §B.2 & B.3) as a
function of time at two different radii, inside (at $r$ shown in each figure)
and outside (at $r=1.2$) the star. This plot allows us to see the transitions
between the different types of modes during the simulation. To interpret the
nature of the observed modes we will later (see §8.3) use the Figure 14, which
displays the different branches of the dispersion relation for upper and lower
modes.

Table 1 lists the details for the runs performed in this study, including
their Mach number ${\mathcal{M}}$, the value of $r_{\rm min}$, resolution and
the type of initial noise pattern used to trigger the acoustic instability in
the BL (see §4.2). To streamline the comparison of runs with different
${\mathcal{M}}$, in the following we focus on simulations that start with the
same initial setup — a particular realization of the block-random noise
pattern, see §4.2; for that reason all these simulations have ’r.a’ in their
label. At the same time, we also provide comparison with simulations using
other types of initial conditions (while keeping ${\mathcal{M}}$ fixed).

###  5.1 A typical ${\mathcal{M}}=9$ run.

We start off by examining mode patterns found in an ${\mathcal{M}}=9$
simulation M09.FR.r.a, which was run at resolution $4096\times 4096$ see Table
1. Figure 2 illustrates the development and operation of the different modes
in this run. The value of ${\mathcal{M}}$ was chosen to demonstrate the
differences with ${\mathcal{M}}=6$ case often studied in the past hydrodynamic
studies of the BL (Belyaev et al., 2012, 2013a); our own exploration of the
${\mathcal{M}}=6$ case is covered in Section 5.2.

In the beginning of the simulation, sonic instability starts off in the form
of an upper $m=27$ mode. At $t/2\pi=25$ (not shown in this figure) the
instability is still in the linear phase, with the perturbed radial velocity
at the level of $\sim 10^{-6}c_{s}$. However, already at $t/2\pi=50$ the
radial velocity reaches saturation at an amplitude of several per cent of
$c_{s}$. At this time, shown in Figure 2Ab-Ac, the pattern of the upper mode
inside the star (with $k_{r}\neq 0$ for $r<1$) is already significantly
affected by the growing lower $m=9$ mode (which has $k_{r}=0$ inside the
star).

What we observe outside the star — large scale spiral arms which freely extend
into the outer disk — may seem consistent with the upper mode still dominating
in the disk. However, the number of spiral arms in the disk is not equal to
azimuthal wavenumber $m=27$ of the upper mode, it is closer to $6$ or $7$. We
will discuss the origin of this pattern in §6.

By the time $t/2\pi=100$ shown in Figures 2Bb-Bc, the picture changes
appreciably. Upper mode weakens considerably (see Fig. 3) and the perturbation
pattern is dominated by a superposition of several lower modes with $m=19,20$
— note their radially elongated pattern inside the star with $k_{r}\approx 0$
as typical for lower modes. These modes have pattern speed $\Omega_{\rm
p}\approx 0.32$(see Fig. 3) for which the corotation radius lies at
$r_{c}=\Omega_{\rm p}^{-2/3}\approx 2.1$ (assuming a Keplerian rotation curve,
which is a good assumption at these radii). However, it looks like these lower
modes launched from the BL have not yet fully filled the resonant cavity
associated with this $\Omega_{\rm p}$ since there are no visible reflected,
inward-propagating waves.

At this time, one can also see the hints of the emergence of an $m=2$ pattern
in the disk, manifesting itself at $t/2\pi=100$ as two broad leading arms for
$r\lesssim 2.5$. The shape of these arms is broadly consistent with what one
would expect from an $\Omega_{\rm p}\approx 0.32$ lower mode, suggesting that
this low-$m$ pattern is somehow related to the $m=19,20$ lower modes present
in the system. This new pattern, as well as the superposition of a weak upper
and a growing lower modes, make the appearance of perturbation in the disk
rather chaotic at this time.

These transitions are accompanied by the evolution of disk surface density
$\Sigma$ and angular frequency $\Omega$ near the stellar surface, as
illustrated in Figure 4 at different moments of time. At around $t/2\pi=100$ a
number of interesting features start to develop in the $\Omega(r)$ profile in
the inner disk, see Figure 4a: an inflection point-like transition at
$\Omega\approx 0.4$ inside the BL, a plateau for $1\lesssim r\lesssim 1.2$,
and a slightly super-Keplerian rotation for $r\gtrsim 1.2$. All these features
are caused by accretion of gas from the disk onto the star, which is driven by
the dissipation of acoustic modes in the disk. Accretion is revealed by the
reduction of $\Sigma(r)$ compared to its initial profile for $r\lesssim 1.7$,
see Figure 4b. This depletion, or gap, is quite substantial near $r=1$
($\Sigma$ drops to $20-30\%$ of its initial value at $r=1.2$) and severely
modifies the radial pressure support in this part of the disk. In agreement
with the equation (8), this has a direct impact on the $\Omega(r)$ behavior:
$\Omega(r)$ develops a sub-Keplerian plateau in the part of the gap where
$\Sigma(r)$ decreases with $r$ (i.e. $r\lesssim 1.2$), and becomes slightly
super-Keplerian outside of this region, since $\Sigma(r)$ increases over a
range of $r$ there. An inflection in $\Omega(r)$ profile inside the BL appears
because the accreted mass accumulates there, modifying the radial pressure
support in the BL. These features will be discussed in more details and across
different values of ${\mathcal{M}}$ in Coleman et al. (in prep.).

Beyond $\sim 100$ orbits the system settles into a less chaotic state
(amplitude of $v_{r}$ variations decreases by $\sim 2-4$ to
$v_{r}/c_{s}<10^{-2}$), which persists until about $t/2\pi=300$, see 2C-F.
During this time the prominent $m=19$ lower mode becomes more and more
coherent both in the disk and inside the star (see Figure 3a,b). Its
relatively low pattern speed $\Omega_{\rm p}=0.315$ corresponds to the
corotation radius $r_{c}\approx 2.16$. This is in good agreement with the
value of the radius at which the outwardly propagating wave crests reach
$k_{r}=0$, turn around, and start propagating inwards, towards the star. Thus,
$m=19$ mode becomes trapped inside the radially extended region — the
_resonant cavity_ — between the stellar surface $r=1$ and the inner Lindblad
resonance which is close to $r_{c}$. The interference of these outward/inward
propagating waves locked inside the resonant cavity at $1<r\lesssim 2$ gives
rise to a regular pattern of criss-crossing leading and trailing spiral arms
confined to the resonant cavity and rotating with a fixed angular frequency
$\Omega_{\rm p}$ on top of the (largely) Keplerian flow in the disk. Such
long-lived regular patterns of lower modes have been identified previously in
Belyaev et al. (2012, 2013a).

Also, during the same period, the aforementioned $m=2$ mode grows in intensity
and very noticeably changes its morphology: it turns into a set of radially
elongated, azimuthally broad perturbation structures that undergo a phase
shift by $\pi$ at around $r=1.4$. This mode has very low pattern speed of
$\Omega_{\rm p}\approx 0.15$ putting its corotation radius around $r\approx
3.5$, still inside our simulation domain. We discuss this mode in more detail
in §8.2.2, but will note here that it persists until about $400$ orbits, co-
existing with the other modes produced at the BL.

At around $300$ orbits an upper $m=23$ mode starts emerging inside the
star555Despite being able to visually identify the $m=23$ mode, it only
becomes the dominant mode inside the star very briefly. Additionally, it’s
amplitude decays very rapidly with increasing radius, preventing our code from
choosing this as a mode to be shown in Fig. 3 (see Section B.3 for more
details).; it has non-zero $k_{r}$ for $r<1$ since its wave crests are at an
angle with respect to the radial direction. At the same time the significance
of the previously dominant lower $m=19$ mode goes down both inside and outside
the star; by $t/2\pi=350$ the associated regular pattern of crossing leading
and trailing waves trapped inside the resonant cavity essentially disappears.

At the same time, outside the star perturbation pattern evolves in a rather
peculiar fashion. Our data do not show an upper $m=23$ mode in the disk — such
a mode would have $k_{r}=0$ right next to the stellar surface. We do see
strong spiral arms with $k_{r}=0$ near the star and extending all the way into
the disk, but there are far fewer of them, only 5 or 6, instead of $m=23$ as
would be appropriate for the regular upper mode (which certainly exists inside
the star). On the other hand, no corresponding $m=5$ or $6$ modes are visible
inside the star; this situation is very similar to that shown in Figure 2Ab-
Ac. The origin of these low-$m$ spiral arms will be covered in §6,7.

Also, starting at around $t/2\pi\approx 400$, a strong $m=6$ perturbation
pattern develops in the disk, confined within $1<r\lesssim 1.25$. It is most
coherent around $t/2\pi\approx 450$, but can be easily traced until the end of
the run (using our automated mode detection algorithm), interfering with the
other modes operating in the system. The nature of this perturbation will be
discussed in §8.2.1.

As time goes by, the number of spiral arms in the disk at $r\gtrsim 1.3$
decreases, with spirals also broadening in azimuthal direction, as if they
were merging together. As a result, after $500$ orbits only 2 or 3 global
trailing arms, somewhat chaotic in appearance, remain in the disk. Such
global, low-$m$ spirals propagating far into the disk are seen in a number of
our runs and represent a novel feature of the BL simulations that will be
discussed in more details in Section 6.4. The disk-star system remains in this
state (with upper $m=23$ mode still visible inside the star and $m=6$ mode
trapped close to the stellar surface) until the end of this simulation at
$t/2\pi=600$.

####  5.1.1 Other ${\mathcal{M}}=9$ and similar runs

Apart from the run M09.FR.r.a described in §5.1, we have also carried out 14
more ${\mathcal{M}}=9$ runs to test the sensitivity of outcomes to our adopted
resolution and initial conditions. With very few exceptions, all runs start
with a brief period (around 50 orbits) of an upper mode operation, switching
to a long period of the lower $m\sim 19-21$ mode dominance, almost always in
conjunction with the $m=2$, $k_{r}=0$ pattern in the disk. In the last
$30-50\%$ of the simulation run time we quite robustly find the $m=4-6$ large
scale spiral arms in the disk (§6) co-existing with the high $m\sim 20-28$
upper mode inside the star (in a few cases the system manages to revert back
to the lower mode with $m\sim 10$ inside the star by the end of the run).
Also, in more than $50\%$ of the runs, after $t/2\pi=300-400$, we see an
$m=4-7$ disk mode confined to the vicinity of the stellar surface (§8.2.1) and
interacting with the other modes. The sensitivity of ${\mathcal{M}}=9$
simulations to numerical inputs is further discussed in §8.4.

We have also carried out ${\mathcal{M}}=10$ and $11$ runs with the same
initial conditions as the ${\mathcal{M}}=9$ run described in §5.1. They show
an evolutionary sequence very similar to that in ${\mathcal{M}}=9$ runs. The
only notable feature worth mentioning is the appearance of a small number of
narrow spiral arms in the beginning of the ${\mathcal{M}}=11$ run. Such arms
are illustrated in Figure 11 and their origin is discussed in §6.

####  5.1.2 Long-term evolution

To provide uniform comparison, all runs in our simulation suite have identical
duration of $600$ inner orbits. However, we did run one ${\mathcal{M}}=9$
simulation for much longer, up to $t/2\pi=3780$, to explore the long-term
evolution of the disk+star system.

Up to $t/2\pi=600$ the evolution of the system is quite similar to all other
${\mathcal{M}}=9$ runs: initial development of an $m=21$ global lower mode
that later becomes chaotic and is replaced by a combination of a number of
global spirals and some low-$m$ modes trapped near the stellar surface. Beyond
that point the system settles into a state in which a high-$m$ upper mode
dominates inside the star, whereas an $m=3-4$ mode trapped near the stellar
surface (to be discussed in §8.2.1) dominates in the disk. The system persists
in this state for very long time, until $t/2\pi\sim 3000$. Then another re-
arrangement happens, with a sequence of several periodic global spirals
merging into a single narrow one-armed spiral, followed by a strong $m=13$
lower mode developing in the system. This is how this simulation ends, giving
us a hint of the complicated long-term evolution mediated by the instabilities
in the BL.

###  5.2 ${\mathcal{M}}=6$ and other low-${\mathcal{M}}$ runs.

Some features of the behavior of the disk-star system for ${\mathcal{M}}=6$
have been discussed in previous studies (Belyaev et al., 2012, 2013a),
although those simulations were typically not advanced for as long as our
current runs. Here we describe the details of a M06.HR.r.a run performed at
the resolution of $2048\times 2048$. The differences with ${\mathcal{M}}=6$
simulations run at other resolutions and/or with different initial noise
patterns are briefly mentioned in §8.4.

In the beginning, at $t/2\pi=25$, this run features a strong upper $m=7$ mode,
see Figure 5a. This mode is rather short lived and goes away already by
$t/2\pi=50$, yielding to the lower mode with $m=10-11$ and $\Omega_{\rm
p}\approx 0.47$. Multiple shocks associated with this mode rapidly drive a
substantial re-arrangement of the disk surface density and $\Omega$ profile in
the inner disk, similar to what is seen in Figure 4. Inside the star this
lower mode (with $m$ slowly evolving from 11 to 9) dominates until the end of
the simulation, see Figures 5b & 6. In the disk it can be traced until
$t/2\pi\lesssim 200$ in the form of criss-crossing leading and trailing
spirals confined to a resonant cavity, $r\lesssim 1.7-2$. This evolution is in
agreement with the past studies (Belyaev et al., 2012, 2013a).

Figure 5: Same as Figure 2 but for the ${\mathcal{M}}=6$ run M06.HR.r.a. Time
corresponding to each snapshot is indicated at the top. See §5.2 for details.
Figure 6: Same as Fig. 3, but for the ${\mathcal{M}}=6$ simulation M06.HR.r.a.
Figure 7: Same as Fig. 2 but for the ${\mathcal{M}}=12$ simulation M12.FR.r.a.
Note a very different radial extent ($r<1.7$) of the Cartesian maps chosen to
illustrate the most important structures. See §5.3 for details.  Figure 8:
Same as Fig. 3, but for the ${\mathcal{M}}=12$ simulation M12.FR.r.a.  Figure
9: Same as Fig. 2 but for the ${\mathcal{M}}=15$ simulation M15.FR.r.a. The
radial range ($r<1.4$) is considerably smaller in Cartesian maps. See §5.4 for
details.  Figure 10: Same as Fig. 3, but for the ${\mathcal{M}}=15$ simulation
M15.FR.r.a.

At around $t/2\pi=225$ (Figure 5c) a new feature emerges in the disk: an $m=2$
mode that replaces the lower $m=9-10$ mode (which still persists in the star).
The low $\Omega_{\rm p}\approx 0.15$ of this mode (Figure 6d) places its
corotation radius at $r_{\mathrm{c}}\approx 3.5$, while its ILR is at
$r_{\mathrm{ILR}}\approx 2.3$, see equation (12). Most of the mode power is
concentrated at $r\lesssim 2.5$, which is compatible with it being trapped
interior to the $r_{\mathrm{ILR}}$.

Despite its apparent radial confinement close to the star, this mode is not a
lower mode; in particular, its location in the dispersion relation Figure 14b
(an orange cross in the lower left corner) is far from the branch
corresponding to the classical lower modes. In fact, it lies close to a branch
of the upper mode dispersion relation, tempting one to identify it as such.
However, it has a morphology that is incompatible with the classical upper
mode behavior — its $k_{r}(r)$ is non-zero for $r\to 1$, it does not produce
spiral arms in the outer disk, etc. This mode is discussed in more detail in
§8.2.1, where we show that it has the same nature as the $m=6$ mode confined
to the inner disk in the ${\mathcal{M}}=9$ run described in §5. Similar
persistent low $m=2-3$ patterns in the disk have been previously reported in
Belyaev et al. (2012). In this particular simulation, the $m=2$ mode dominates
fluid perturbation in the disk almost until the end of the simulation, when it
becomes more chaotic and a strong $m=1$ perturbation (discussed in §6)
eventually develops on top of it.

Other runs carried out for similar values of ${\mathcal{M}}$ ranging from $5$
to $8$ paint a picture similar to the ${\mathcal{M}}=6$ case: lower modes
dominating inside the star, initially with $m=11-13$, then dropping to
$m=5-9$. At the same time, the trapped low-$m$ modes mentioned above (e.g.
$m=2$ for ${\mathcal{M}}=5$ and $m=4$ for ${\mathcal{M}}=8$) dominate in the
disk, typically with very low $\Omega_{\rm p}$ allowing them to extend far
out.

An interesting exception is the ${\mathcal{M}}=7$ case, in which the disk is
dominated for the majority of the run by a small number, $2-5$, of global
spiral arms propagating out, see Figure 11b. This picture is very similar to
what is observed in the end of ${\mathcal{M}}=9$ run described in §5.1. These
modes may look like upper modes, but their $\Omega_{\rm p}\approx 0.75$ is too
high for an upper mode (which can be deduced from Figure 14c); the nature of
these modes will be addressed in §6. This run also features a superposed
prominent $m=16$ lower mode (twice the $m$ of the lower mode dominating inside
the star) locked in the resonant cavity $1<r\lesssim 1.8$, which makes the
${\mathcal{M}}=7$ run look even more different from ${\mathcal{M}}=6$ or $8$
simulations. These observations demonstrate that a particular mix of modes
emerging in simulations is not a strictly monotonic function of
${\mathcal{M}}$.

###  5.3 ${\mathcal{M}}=12$ runs.

Runs at Mach number ${\mathcal{M}}=12$ show surprisingly little variation in
the outcomes, regardless of the initial conditions. All of them have been run
at the resolution of $4096\times 4096$. Below we will briefly describe the
evolution of the system in the M12.FR.r.a simulation, but all of
${\mathcal{M}}=12$ runs exhibit a very consistent pattern dominated by the
lower modes.

In the beginning (at $t/2\pi=50$) the system features a dominant lower $m=9$
mode inside the star with $k_{r}=0$, see Figure 7. At the same time, there are
multiple spiral arms in the disk freely propagating towards the edge of the
simulation domain, which are very similar in appearance to the disk modes
present in the beginning of the ${\mathcal{M}}=9$ run, see Figure 2Ab,Ac and
§6. The absence of the disk counterpart of the strong lower $m=9$ mode visible
inside the star may be caused by the fact that it has a high pattern frequency
$\Omega_{\rm p}\approx 0.8$ (see Figure 8): already by that time $\Omega(r)$
profile develops a sub-Keplerian plateau just outside $r=1$ (similar to that
seen in Figure 4a) with the maximum $\Omega$ close to this $\Omega_{\rm p}$.
As a result, the resonant cavity in which this mode would be confined in the
disk is too narrow, suppressing its signal.

By $t/2\pi=125-150$ perturbation pattern changes and the dominance switches to
the global lower $m=16$ mode. This mode has a low $\Omega_{\rm p}\approx
0.45$, which guarantees a substantial radial extent of the resonant cavity in
which the mode is trapped: its corotation radius is $r_{\mathrm{c}}\approx
1.7$. This mode persists for a long time, roughly until $t/2\pi=450$, see
Figure 7.

Beyond that point, by 550 orbits, the lower $m=16$ mode fades away, and the
dominant mode inside the star becomes the $m=11$ lower mode with $\Omega_{\rm
p}\approx 0.65$ and ILR at $r_{\mathrm{ILR}}\approx 1.25$. However, in the
disk we do no find 11 criss-crossing waves confined between the star and
$r_{\mathrm{ILR}}$, as one would expect for a global $m=11$ mode. Instead, we
see an $m=7$ perturbation pattern confined to $1<r\lesssim 1.2$. This is quite
reminiscent of the situation described in §5.2, where we found the disk to
feature an unusually low-$m$ mode ($m=2$ in that case) trapped near the
stellar surface, simultaneously with a higher-$m$ lower mode inside the star.

In addition to the $m=7$ mode trapped near the star, Figure 7Fb,Fc also
reveals presence of the large scale global spiral arms extending to the outer
boundary of our domain. Although somewhat less coherent, these spirals are
similar to the spiral arms found closer to the end of ${\mathcal{M}}=9$ run,
see Figure 2. They will be discussed further in §6. The pattern seen at 550
orbits stays unchanged until the end of the simulation.

Remarkably, all other ${\mathcal{M}}=12$ runs show very similar regular
behavior, down to minor details: a dominant $m=9$ lower mode with high
$\Omega_{\rm p}\approx 0.8$ in the beginning, changing around $t/2\pi=125-175$
to lower $m=16$ mode confined between $r=1$ and $r\approx 1.5-1.6$. Only the
time at which $m=16$ lower mode fades away shows some variation, roughly
between $t/2\pi=400$ and $550$.

The regular and stable behavior exhibited by ${\mathcal{M}}=12$ runs is
somewhat unique. The ${\mathcal{M}}=10,11$ runs exhibit the behavior which is
much closer to ${\mathcal{M}}=9$ simulations. At the same time
${\mathcal{M}}=13,14$ runs are also very different from ${\mathcal{M}}=12$
ones, and resemble ${\mathcal{M}}=15$ run, see next.

###  5.4 ${\mathcal{M}}=13-15$ runs.

Starting at ${\mathcal{M}}=13$ the simulation outcomes change quite
dramatically, with a general pattern that persists at least until
${\mathcal{M}}=15$. Simulations for ${\mathcal{M}}=13-15$ have resolution
$8192\times 8192$ and initial conditions (’r.a’) identical to those used in
the ${\mathcal{M}}=6,9,12$ runs discussed earlier. Their general behavior is
illustrated here using the ${\mathcal{M}}=15$ run M15.FR.r.a in Figure 9.

Pretty much all the time since the start of the simulation, after the sonic
instability reaches saturation, and until the very end, the star supports a
strong upper mode with $m$ varying between 12 to 25 at different moments of
time (e.g. $m=19$ for $t/2\pi>200$), see Figure 10.

Outside the star situation is more complicated. In the beginning (see Fig.
9Ab,Ac), the disk features around 8 global spiral arms that extend from the
stellar surface all the way to the outer domain of the simulation. These arms
are quite narrow in azimuthal direction. Given that their number does not
coincide with the azimuthal wavenumber of the upper mode present in the disk,
they cannot be the manifestation of a global upper mode in the star-disk
system. Later on, at $175$ orbits (panels Bb,Bc), the morphology of the spiral
arms changes: their number goes down to about 5 while their azimuthal width
increases. At later times the arms become less coherent and their number
steadily decreases, with just one or two reasonably strong arm-like structures
present in the end of the simulation (accompanied by numerous weaker arms).
The origin of these spiral structures will again be covered in §6.

Starting at around 275 orbits and until the end of the simulation, the disk
also exhibits a concentration of perturbation power in the radial range
$1<r\lesssim 1.17$. It comes in the form of crossing wakes with positive and
negative $k_{r}$, locked inside the resonant cavity. Despite the appearance
similar to the lower mode in the disk, we interpret these waves to be the
high-${\mathcal{M}}$ analogues of the low-$m$ modes trapped near the stellar
surface that we saw in other runs at lower $M$, see Figs. 5 ($m=2$ after 300
orbits for ${\mathcal{M}}=6$) and 5 ($m=7$ at 550 orbits for
${\mathcal{M}}=12$). A notable difference with the other runs is that for
${\mathcal{M}}=15$ the azimuthal wavenumber of these wakes is around 20, much
higher than at lower ${\mathcal{M}}$.

A very similar evolutionary sequence is observed in the ${\mathcal{M}}=13,14$
runs — persistent dominance of the upper modes inside the star, evolving
pattern of the global spirals in the disk, multiple crossing wakes at
$r\lesssim 1.2$. The only notable trend that we see is the increase of the
azimuthal wavenumber of the observed features with growing ${\mathcal{M}}$.

##  6 Vortex-driven modes in the near-BL region

In addition to spatial distributions of $rv_{r}\sqrt{\Sigma}$, which
illustrate the amplitude of the wave-like perturbation, we also examined the
maps of vortensity (or potential vorticity, related to the vorticity
$\mathbf{\omega}\equiv\mathbf{\nabla}\times\mathbf{v}$)

$\displaystyle\zeta\equiv\dfrac{\mathbf{\omega}}{\Sigma}=\dfrac{\mathbf{\nabla}\times\mathbf{v}}{\Sigma},$
(17)

which are shown in the $r-\phi$ coordinate plane in Figures 2, 5, 7, 9 (left
columns).

These maps reveal that many of the morphological structures observed in our
simulations and mentioned in §5 are, in fact, caused by the localized
structures in the spatial distribution of vortensity that emerge in the near-
BL region. Quite generally, we find two types of vortensity structures that
give rise to global waves in disks. Their typical appearance is illustrated in
Figure 11, where we plot both vortensity and $rv_{r}\sqrt{\Sigma}$ in the
$r-\phi$ coordinate plane for a couple of representative runs.

Figure 11:  Two main types of the localized vortensity structures (vortensity
change relative to its initial value) emerging in our simulations (left) and
the fluid perturbation $rv_{r}\sqrt{\Sigma}$ that they drive in the disk
(right). (a,b) Compact vortices present close to the BL, around $r\approx 1.1$
in the ${\mathcal{M}}=11$ run M11.FR.r.a at $t/2\pi=173$. (c.d) Regular
“rolls" forming in the inner disk around $r\approx 1.5$ in the M07.FR.r.a run
at $t/2\pi=450$. White dashed curves (on the right) represent the locations of
the wave crests according to the WKB dispersion relation (10), and their
association with the corresponding vortex structures (on the left). See §6 for
further discussion.

The first type of vortensity structures reveals itself in Figure 11a (showing
a snapshot of ${\mathcal{M}}=11$ run at 175 orbits) as sharp, elliptical,
anticyclonic features in $\zeta$ map located very close to the stellar
surface. We call these structures simply vortices. They appear rather narrow
in azimuthal direction but this is simply a result of the aspect ratio chosen
in this figure — in reality they are rather elongated azimuthally.
Nevertheless, these vortices are typically well-isolated in azimuthal
direction while sharing the same radial range $1<r\lesssim 1.2$ like beads on
a wire. Looking now at Figure 11b, one immediately notices that azimuthal
positions of these vortices coincide extremely well with the starting
azimuthal locations (at $r\approx 1$) of a number of sharp, narrow spiral arms
that propagate out into the inner disk. The strength of the arm (amplitude of
its $rv_{r}\sqrt{\Sigma}$) appears to scale with the size of the vortex to
which the arm in connected. The number of arms — about 7 — is the same as the
number of noticeable vortices in panel (a) of the figure. One can also see
that the global spirals in the disk co-exist with the global lower $m=13$
acoustic mode (easily visible inside the star) — a very different type of the
wave-like perturbation.

Second type of vortensity structure is illustrated in Figure 11c, which shows
a snapshot of the ${\mathcal{M}}=7$ run at 450 orbits. This vortensity map
reveals a set of four azimuthally elongated “rolls", as we call666We often
collectively refer to isolated vortices and rolls as just “vortices”. these
structures, which are centered at $r\approx 1.5$ and have approximately equal
azimuthal extent. Unlike vortices, the rolls are not isolated and touch each
other, collectively covering the full circumference of the disk. Another
difference with respect to vortices is that the rolls are always found at some
separation from the stellar surface; in Figure 11c they occupy a radial range
$1.4\lesssim r\lesssim 1.6$.

Comparing panels (c) and (d) of Figure 11 one can see that each roll is
associated with a broad spiral arm easily visible in $rv_{r}\sqrt{\Sigma}$.
Just as the rolls, the spiral arms are azimuthally broad, which distinguishes
them from the narrow arms launched by the vortices. This results in a distinct
$m=4$ pattern of global spirals in the outer disk. The leading half of each
roll is connected to the $v_{r}>0$ part of the corresponding spiral arm, while
the opposite is true for the trailing half of the roll, indicating their
anticyclonic nature (same as vortices). Also, Figure 11b shows that the roll-
driven spirals arms can naturally co-exist with the acoustic modes, in that
case a lower mode which is rather strong in the disk out to $r\approx 1.75$.

Since the starting points of the global spirals (their azimuthal locations at
$r\to 1$) always coincide with the positions of their associated
vortices/rolls, the pattern speed of the spiral arms in our runs is the same
as the orbital frequency of these vortensity structures. As both vortices and
rolls are passively advected with the fluid, angular frequency $\Omega(r)$ of
the disk fluid at their orbital radii sets $\Omega_{\rm p}$ of their global
spirals. Given that rolls are more distant from the stellar surface than the
vortices, $\Omega_{\rm p}$ of the spirals associated with vortices is higher
than $\Omega_{\rm p}$ of the spirals related to rolls.

The two kinds of vortensity structures described above emerge at different
times in many (but not all) of our simulations, and can even co-exist for
brief periods of time. Moreover they tend to evolve and exhibit transitions
from one type of structure to another. We now briefly describe the typical
evolutionary patterns of vortex-driven modes in runs with different
${\mathcal{M}}$.

###  6.1 Vortensity structures in ${\mathcal{M}}=9$ runs

Our fiducial ${\mathcal{M}}=9$ run illustrated in Figure 2 shows a set of
isolated vortices emerging near the stellar surface ($r\lesssim 1.25$) by
$t/2\pi=50$. These vortices are the true reason behind a set of strong global
spirals that are visible in panel (Ab) of this figure (and not the upper mode,
as mentioned in §5.1). They persist at 100 orbits, and their associated spiral
arms are discernible in the disk even in the face of a strong lower mode that
develops in the system. However, beyond that point vortices get washed out and
merge with each other — their number goes down. Correspondingly, the
characteristic narrow spiral arms in the disk disappear leaving only the lower
mode.

Beyond 400 orbits a new transition takes place in the system — a set of rolls
starts to emerge at $r\approx 1.3$. At $t/2\pi=450$ one can see 5 regular,
roughly equally spaced rolls connected to a set of 5 strong global spirals in
the disk. These rolls evolve by merging with each other: only 3 of them remain
at 500 orbits (still at roughly equal azimuthal separation from each other),
connected to an $m=3$ set of global spirals in the disk, see panels (Ja)-(Jc).
Only 2 rolls (and spirals) remain at 550 orbits, separated by roughly
$180^{\circ}$. However, by 600 orbits they drift azimuthally towards one
another (while remaining at the same radial distance) and would merge into a
single roll if we ran this ${\mathcal{M}}=9$ simulation for longer.

###  6.2 Vortensity structures in the high-${\mathcal{M}}$ runs

At higher values of ${\mathcal{M}}>9$ we typically find rolls to emerge quite
early. For example, in Figure 7Aa illustrating an ${\mathcal{M}}=12$ run
described in §5.3, a number (7 or 8) of rolls become apparent at $r\approx
1.17$ already at 50 orbits, when a number (9 or 10) of strong, isolated
vortices is still present closer to the star. Careful examination of the panel
(Ab) of that figure reveals two complexes of global spirals — one due to the
vortices next to the BL and another one associated with the rolls, forming
further out in the disk. They can be distinguished by their different pitch
angles: roll-driven spiral have lower $\Omega_{\rm p}$ and are less tightly
wound than the vortex-driven spiral arms, which have higher $\Omega_{\rm p}$.
Because of the difference of their $\Omega_{\rm p}$, the two sets of spirals
drift azimuthally relative to each other.

Co-existence of rolls and vortices persists in this ${\mathcal{M}}=12$ run for
quite a while, with both types of structures (and their associated spirals)
visible up to $250$ orbits. However, the number of both vortensity
substructures goes down as they merge, while maintaining roughly the same
radial distance. Vortices near the stellar surface stop being visible only
after $\approx 300$ orbits, see panel (Da).

In this particular simulation vortensity distribution also tends to develop a
banded structure after about 150 orbits. Radially narrow and almost azmuthally
symmetric bands in $\zeta$ maps appear to give rise to weaker rolls at larger
separation from the star. This complicated radial distribution of $\zeta$ goes
away only at the end of the simulation, although the rolls at $r\approx 1.1$
still persist in some form.

A rather similar evolution of vortensity structures is also found in
${\mathcal{M}}=15$ run described in §5.4. Left row of Figure 9 shows strong
vortices early on (panel Aa), which co-exist with a number of rolls later on
(panels Ba and Ca), with rolls dominating after $\approx 300$ orbits. These
vortensity structures explain the global spirals visible in the maps of
$rv_{r}\sqrt{\Sigma}$ at various degree of coherence throughout the
simulation.

###  6.3 Vortensity structures in the low-${\mathcal{M}}$ runs

Situation is quite different in our runs with low values of ${\mathcal{M}}<9$.
We find that only the ${\mathcal{M}}=7$ run shows the development of strong
vortices and, subsequently, rolls, reminiscent of the ${\mathcal{M}}=9$ run;
similarity of the perturbation morphology between the ${\mathcal{M}}=7$ and 9
runs has been previously noted in §5.2. On the other hand, ${\mathcal{M}}=8$
run does not show any strong or long-lasting azimuthal vortensity structures —
the distribution of $\zeta$ in this run looks quite axisymmetric throughout
its duration. And the simulations with ${\mathcal{M}}=5$ and 6 do develop
rather peculiar vortensity structure, illustrated in the left column of Figure
5, which is very distinct from the higher ${\mathcal{M}}$ runs.

This ${\mathcal{M}}=6$ run (M06.HR.r.lc.a) shows near-stellar surface vortices
only for a very brief interval of time around 75 orbits (not shown). And soon
after a strong $m=2$, low-$\Omega_{\rm p}$ mode (described in §5.2) appears in
the disk, the distribution of $\zeta$ develops a characteristic wavy $m=2$
pattern, in which contours of constant $\zeta$ oscillate in $\phi$ with large
radial amplitude ($1.1\lesssim r\lesssim 1.4$). These oscillations result from
passive advection of vortensity by the periodic large amplitude perturbations
of $v_{r}$ associated with the $m=2$, low-$\Omega_{\rm p}$ mode.

Later on, at 425 orbits, one notices two localized vortices (blue dots in
Figure 5Ea near $\phi/\pi\approx 0.2$ and 1.6) appearing quite far from the
star, around $r=2.1$. These vortices drift radially inwards and eventually
merge, resulting in a single vortex visible at 525 orbits at $r\approx
1.9,\phi/\pi\approx 0.2$, which is responsible for the strong $m=1$
perturbation in $rv_{r}\sqrt{\Sigma}$ that develops in the outer disk for
$r\gtrsim 1.8$. However, careful examination of the vortensity patterns at
larger radii reveals that these vortices form early on near the outer boundary
of our simulation domain, as a result of a numerical artefact related to our
outer boundary condition. Their subsequent inward drift is a natural outcome
of the vortex dynamics in the disk, see Paardekooper et al. (2010).

This sequence of vortensity evolution is very typical for our
${\mathcal{M}}=5$ and 6 runs: we see essentially no vortensity structures
produced near the stellar surface (except for the wavy advective patterns),
but at late time vortices resulting from numerical artefacts at the outer
boundary migrate in and disturb the global vortensity distribution. However,
starting at ${\mathcal{M}}=7$ and higher we never see these numerical
artefacts appear in our runs.

Figure 12:  Example of a single-armed spiral emerging in one of our
${\mathcal{M}}=12$ simulations M12.FR.random.a White dashed curve shows the
analytical fit given by equation (19), which matches the shape of the spiral
density wave very well. Note a prominent global lower mode active in the inner
disk, at $r\lesssim 1.5$. See §6.4 for details.

###  6.4 Emergence of narrow, single-armed spirals

In roughly one third of our simulations we observe vortices or rolls to
gradually merge into a single strong, coherent vortex, which launches a
narrow, single-armed spiral density wave in the disk. A typical example is
shown in Fig. 12 illustrating one of our ${\mathcal{M}}=12$ runs
(M12.FR.random) at 400 orbits. These narrow spiral features form almost
exclusively in runs with ${\mathcal{M}}\geq 7$. This is because, as discussed
earlier in §6.1-6.3, single isolated vortices tend to form only in simulations
with sufficiently high values of ${\mathcal{M}}$. These arms are rather long-
lived and can last for $\sim 100$ orbits. They are important because they can
lead to interesting observational manifestations in the time domain.

Such single-armed features have much smaller azimuthal width than the $m=1$
patterns emerging in some of the low-${\mathcal{M}}$ runs, e.g. the one shown
in Figure 5(Fa)-(Fc). They closely resemble the spiral arms that appear in
simulations of protoplanetary disks with embedded, moderately massive planets.
Because of the narrow azimuthal width, such arms must be superpositions of a
number of high-$m$ acoustic modes (as in the case of planet-driven spirals),
with pattern speed $\Omega_{\rm p}$ set by the angular frequency $\Omega$ of
their parent vortex (or roll).

This allows us to better understand the shape of these arms. Indeed, for $m\gg
1$, the first term in the right-hand side of the WKB dispersion relation (10)
can be neglected (at large radii $\Omega(r)$ also becomes small compared to
$\Omega_{\rm p}$), allowing us to express

$\displaystyle k_{r}(r)\approx-m\frac{\Omega(r)-\Omega_{\rm p}}{c_{s}(r)},$
(18)

where we chose sign so that $k_{r}>0$ in the outer disk, far from the BL.
Integrating the relation (11) with this expression for $k_{r}$ gives the
equation for the shape of the wave crest in the form (Rafikov, 2002)

$\displaystyle\varphi(r)=\varphi_{\rm ref}+\int\limits_{r_{\rm
ref}}^{r}\frac{\Omega(r^{\prime})-\Omega_{\rm
p}}{c_{s}(r^{\prime})}dr^{\prime},$ (19)

where $\varphi_{\rm ref}$ is the azimuthal coordinate of the wave crest at
some reference radius $r_{\rm ref}$.

The fact that the modes with $m\gg 1$ have $\varphi(r)$ essentially
independent of $m$ means that these modes can constructively interfere,
maintaining the one-armed profile in a narrow azimuthal range over large
radial intervals. For disk-planet interaction this observation was made
previously by Ogilvie & Lubow (2002) and Rafikov (2002), whereas Bae & Zhu
(2018) and Miranda & Rafikov (2019a) pointed out that this coherence works
particularly well in the outer disk (whereas in the inner disk it gets
gradually lost). This is relevant for our case since the narrow global arms
that we observe are exterior to the vortices that launch them.

Single-armed spirals that we see in our runs were not observed in previous
simulations of the BLs (e.g. Belyaev et al., 2012, 2013a). Many of these
earlier simulations did not extend over the full $2\pi$ in the azimuthal
direction, which would both not support single-armed features and affect the
emergence and evolution of vortices driving the single-armed spiral. Other
simulations that did cover the full $2\pi$ in azimuth had limited radial
extent ($r_{\rm max}=2.5R_{\star}$), which likely prevented them from
revealing single-armed spirals. To verify this hypothesis we preformed test
runs in which we varied $r_{\rm max}$ and, indeed, did not find any single-
armed spirals in simulations with $r_{\rm max}<3R_{\star}$. This suggests that
a large radial extent is necessary for capturing the development of such wave
phenomena in simulations.

##  7 Origin of the vortex-driven modes

In §6 we uncovered a clear connection between the multiple spiral arms and the
vortex-like structures in the near-BL part of the disk. In particular,
azimuthal locations of vortices in the near-BL region coincide with the
launching sites of the major global spiral arms in the disk. The multiplicity
and pattern speeds of these spiral arms are controlled by the number of the
corresponding vortices and their radial location. This naturally raises a
question of the underlying cause behind this relationship.

Local perturbations of vortensity, confined both in radius and azimuthal
angle, which we call vortices for simplicity, are known to trigger density
waves in accretion disks through the velocity perturbations that they induce
in the underlying disk. This has been demonstrated both numerically (Li et
al., 2001; Johnson & Gammie, 2005) and through detailed analytical exploration
(Heinemann & Papaloizou, 2009; Paardekooper et al., 2010). In many ways the
action of vortices is similar to that of planets (or other massive orbiting
perturbers), that launch density waves via their gravitational coupling to the
disk at the Lindblad resonances (Goldreich & Tremaine, 1980). Thus, as long as
vortices are present in the inner part of the disk, the excitation of global
spiral arms propagating over large distances in the disk is quite natural.

However, this brings up the next obvious question: what causes the emergence
of vortices in the near-BL region in the first place? We address this question
next.

###  7.1 Origin of vortices in the near-BL region

Examination of vortensity evolution in Figs. 2, 5, 7, 9 demonstrates that in
the beginning of the run $\zeta$ grows in the near-BL region above its initial
value, which is equal to $\Omega_{K}/(2\Sigma)$ and is radially constant in
the disk region for our initial conditions, see §4.2. Evolution of $\zeta(r)$
is illustrated in Fig. 13, where we plot the azimuthally-averaged profiles of
the vortensity at different moments of time for the four simulations described
in §5. One can see that in all four runs vortensity experiences substantial
evolution in the near-BL region. This raises a possibility of a Rossby wave
instability (RWI), which operates in presence of radially structured
vortensity, being triggered in this part of the disk.

The importance of RWI in astrophysical disks has been pointed out by Lovelace
et al. (1999) who demonstrated, in particular, that infinitesimal
perturbations can grow exponentially provided that the underlying radial
profile of $\zeta$ has an extremum. It was subsequently studied by a number of
authors both analytically (Li et al., 2000; Ono et al., 2016) and numerically
(Li et al., 2001; Johnson & Gammie, 2005; Ono et al., 2018). A natural outcome
of the nonlinear stage of the RWI is the formation of multiple vortices (each
of them launching their individual spiral arms) with their subsequent merger
into a single major vortex (Ono et al., 2018). This sequence of events is
precisely what we observe in our simulations. Previous studies typically
considered the RWI triggered by features in $\zeta(r)$ profile arising due to
localized steps, bumps, or gaps in the surface density. The latter — a drop in
$\Sigma$ — is always found in our simulations, see Figure 4.

Note that, according to Papaloizou & Lin (1989), in barotropic disks, such as
the globally isothermal disk considered in this work, exponentially growing
modes of the RWI require the minima of the vortensity profile to exist in the
disk. However, a smooth drop in $\Sigma$ in a Keplerian disk would give rise
to a maximum of $\zeta(r)$, which should be stable according to Papaloizou &
Lin (1989). Nevertheless, in many of our simulations we also find the minima
of $\zeta(r)$ to emerge quite naturally. Fig. 13 shows that at different
moments of time ${\mathcal{M}}=9,12,15$ runs exhibit $\zeta(r)$ profiles with
(multiple) local minima, which would give rise to RWI. Interestingly, the
profile of $\zeta(r)$ in ${\mathcal{M}}=6$ run tends to show only a single
broad maximum and no minima, see Fig. 13d. This is consistent with the lack of
the near-BL vortices in the low-${\mathcal{M}}$ runs, see §6 and Fig. 5. At
the same time, it should be remembered that derivation of the RWI excitation
criterion in Papaloizou & Lin (1989) was based on certain simplifying
assumptions: static, axisymmetric background vortensity profile, infinitesimal
perturbations, etc. In real near-BL region, we see that $\zeta(r)$ is
generally non-axisymmetric, rapidly changes in time, and is being constantly
perturbed by the acoustic waves, which are at least mildly nonlinear. For
these reasons the RWI criterion formulated in Papaloizou & Lin (1989) may not
be directly applicable for interpreting the results of our simulations, even
if it seems to work qualitatively.

###  7.2 Vortensity evolution in the near-BL region

A final step in closing the logical loop of understanding vortex-driven modes
is to explain the apparent growth of vortensity in the vicinity of the stellar
surface that we observe in Fig. 13, which is necessary for the RWI to get
triggered. In barotropic disks $\zeta$ is strictly conserved,
$\mathrm{d}\zeta/\mathrm{d}t=0$. However, conservation of $\zeta$ gets broken
in presence of shocks. In our BL simulations mildly nonlinear lower modes
evolve into shocks very naturally, driving the growth of vortensity within the
resonant cavity where they are trapped. Upper modes do not seem to be very
efficient at driving the growth of $\zeta$.

Figure 13: Evolution of the azimuthally-averaged vortensity profile
$\zeta(r)=\langle\omega/\Sigma\rangle$ in simulations with different
${\mathcal{M}}$ (labeled in each panel). Profiles of $\zeta(r)$ at different
times are color-coded according to the legend in the upper panel. See §7.1 for
the discussion.

The local rate at which vortensity evolves due to shocks depends on a variety
of factors: multiplicity of the waves (i.e. azimuthal wavenumber $m$ of the
lower modes), their amplitude, their pitch angle (which itself is determined
by the pattern speed $\Omega_{\rm p}$ of the underlying modes), see Kevlahan
(1997), Lin & Papaloizou (2010) and Dong et al. (2011). In addition to $\zeta$
production at shocks, vortensity is also passively advected into the star as a
result of mass accretion. Intricate interplay between these processes
determines the radial distribution of the vortensity growth rate in such a way
that $\zeta(r)$ exhibits a complicated structure in the near-BL zone, allowing
the RWI to operate.

To summarize, vortex-driven modes emerge as a result of multi-stage process
driven by the sonic modes. First, sonic instability in the supersonic shear
layer produces (lower) acoustic modes. Second, these modes, being mildly
nonlinear, evolve into shocks and drive vortensity production within the
resonant cavity near the stellar surface. Third, accumulation of vortensity
creates the conditons for excitation of the RWI that gives rise to multiple
vortices in its nonlinear phase. Finally, each vortex launches a density wave
that propagates out from the BL region as a vortex-driven spiral arm. A very
similar sequence of steps occurs in tidal coupling of protoplanetary disks
with massive embedded planets (Koller et al., 2003; Li et al., 2005; de Val-
Borro et al., 2007; Lin & Papaloizou, 2010): planet-driven density waves shock
near the planet, modifying the vortensity profile and triggering RWI, which
produces vortices at the edges of the forming gap, with secondary spiral arms
being driven by such vortices in the disk.

###  7.3 Compact vortices vs “rolls"

The two main types of vortensity structures that we identify in our
simulations — isolated vortices and rolls — are quite different from one
another in a number of ways.

First, isolated vortices have smaller azimuthal extent and are less regular
than rolls: rolls tend to appear as azimuthally periodic chains of regular
patterns of vortensity, whereas isolated vortices are more irregular in their
morphology. Second, isolated vortices are most prominent in the very beginning
of the simulation, whereas rolls appear only after sufficient time has passed
for the disk surface density and vortensity structure to be substantially
modified near the BL. Third, isolated vortices exist only very close to the
BL, at $r\to 1$, whereas rolls tend to form at some separation from the BL,
typically at $r\sim(1.1-1.5)$.

At least some of these observations can be interpreted by comparing $\zeta$
maps in Figs. 2, 5, 7, 9 with the radial profiles of $\zeta$ in Figure 13. In
particular, in the beginning of the simulations $\zeta$ profiles show strong
peak of vortensity at $r\to 1$, which is the natural consequence of the
initial sharp gradient of $\Omega(r)$ across the BL. These peaks are what
gives rise to vortices early on in the simulation. As the run progresses and
the BL broadens, radial gradients of $\Omega(r)$ diminish, lowering $\zeta$
peaks at $r\approx 1$ and reducing the significance of the strong, sharp
vortices over time.

Figure 13 also shows other vortensity peaks, appearing in the disk at some
distance from the BL at later stages. It is easy to see that the radial
locations of these peaks coincide with the locations of the chains of rolls
that emerge in our runs at roughly the same moments time. In other words,
rolls in the inner disk appear to be driven by vortensity generation at weak
shocks, caused by the nonlinear evolution of the near-BL acoustic modes.

Given this difference in origin, one may wonder if isolated vortices are
purely an artefact of our initial conditions in the form of a sharp
$\Omega(r)$ gradient. This is only partly true, since such gradient has to
persist through our runs because of the $\Omega$ drop in the BL. The amplitude
of this gradient (which directly translates into the height of the vortensity
peak) is a strong function of ${\mathcal{M}}$ since the BL width is a
sensitive function of the Mach number and scales roughly as
${\mathcal{M}}^{-2}$, see Belyaev et al. (2012) and Coleman et al. (in prep.).
For that reason vortices at $r\to 1$ are never strong in our
low-${\mathcal{M}}$ runs. Whereas in the high-${\mathcal{M}}$ runs, in which
the BL is narrow, $\zeta(r)$ maintains a peak near the star and vortices at
$r\approx 1$ tend to be long lived; see e.g. $\zeta$ distribution in the left
panels of Figure 9, where some vortex-like structures are present right next
to the BL throughout the full duration of this ${\mathcal{M}}=15$ run.

###  7.4 Robustness of the vortex-driven modes

Formation of a depression in $\Sigma$ near the stellar surface and the
associated peak of $\zeta$ appears essential for providing the conditions for
vortex/roll excitation in the near-BL region. Our simulations are inviscid,
and such forming gap does not get replenished by the material arriving from
larger radii in the disk. However, in real accretion disks there is mass
inflow (e.g. due to the MRI), which would tend to refill the gap with gas
brought in from larger radii, and might prevent vortex-driven modes from
appearing. This possibility may be difficult to realize because of the
efficiency with which sonic modes transport mass near the stellar surface. It
is plausible that even with the continuous mass inflow from larger radii,
sonic modes will still be able to modify $\Sigma(r)$ near the star, sufficient
to keep RWI going. The gap does not need to be very deep for vortex-driven
modes to emerge; for example, ${\mathcal{M}}=12,15$ runs exhibit rather
shallow (only $\sim 30\%$ deep) gaps but still support vortex-driven modes.

Another potential issue with the vortex-driven waves is the fact that our
simulations are 2D. While vortices can certainly form in 3D simulations, there
is an ongoing debate about their longevity in realistic protoplanetary disks
with vertical structure (Barranco & Marcus, 2005; Lithwick, 2009; Lesur &
Papaloizou, 2009; Meheut et al., 2012b; Meheut et al., 2012a; Lin, 2012; Lin &
Pierens, 2018). In this regard we note that our own 3D simulations (to be
analyzed in the future) do show the emergence of the vortex-driven modes and
their survival over long time intervals.

##  8 Discussion

The main goal of this work is a systematic exploration of the acoustic mode
activity in the vicinity of the BL. We do this in a rather simple but easy to
control setup, with the Mach number ${\mathcal{M}}$ being the only free
parameter distinguishing our runs. The initial distribution of the disk
surface density is chosen to ensure a flat profile of vortensity $\zeta$, to
avoid possible biases related to the initial conditions.

The equation of state that we use in this work — globally isothermal — greatly
simplifies our analysis of the angular momentum and mass transport in the
near-BL region that will be presented in Coleman et al. (in prep). This is
crucial since recent studies (Lin, 2015; Miranda & Rafikov, 2019b, 2020) have
shown that the often used non-barotropic, locally isothermal equation of state
leads to non-conservation of the angular momentum flux carried by the waves
even in the absence of explicit dissipation. Use of the barotropic globally
isothermal equation of state also allows us to not worry about the long-term
heating/cooling effects on the disk temperature, reducing the number of
physical parameters that we need to follow.

While carrying out this exploration we discovered new types of hydrodynamic
wave-like phenomena that emerge in the inner disk near the stellar surface.
Probably the most interesting are the vortex-driven waves (§6), and we already
covered their origin and properties at length in §7. In the following we
provide a discussion of several other notable results of our simulations,
among them the analysis of the regular acoustic (§8.1) and other (§8.2) modes,
as well as the dependence of their harmonic content on ${\mathcal{M}}$.

Figure 14: Dispersion relation, i.e. a plot of a pattern speed $\Omega_{\rm
p}$ vs. the azimuthal wavenumber $m$, for a number of representative
simulations with different values of the Mach number ${\mathcal{M}}$ (labeled
in the panels). Runs for all integer values of ${\mathcal{M}}$ between 5 and
15 are shown, with two runs for ${\mathcal{M}}=9$ to illustrate the effects of
resolution. We show the $(m,\Omega_{\rm p})$ pairs for modes automatically
detected in the disk only (at $r=1.2$, yellow pluses), in the star only
($r=(r_{\rm min}+R_{\star})/2$, blue circles), and global modes present both
in the disk and the star with the same $m$ and $\Omega_{\rm p}$ at a given
time (green stars). Red dotted and dot-dashed curves display the dispersion
relations for the upper (§3.1) and lower (§3.2) modes, correspondingly; black
dotted and dot-dashed curves are their respective higher-order azimuthal
harmonics. Horizontal dashed lines show the maximum value of $\Omega(r)$ in
the disk at late times. See §8.1.1 and other text for details.

###  8.1 Acoustic modes

Acoustic modes excited by supersonic shear in the BL are interesting not only
on their own but also because they are the ultimate drivers of accretion onto
the central object (Coleman et al., in prep.) and are intimately involved in
generation of other types of modes, see §7. Both lower and upper modes
described in §3 are observed in our simulations. Only very rarely do we see
the third — middle —- mode described in Belyaev et al. (2013a) temporarily
appear early on in some of our runs.

We generally find the upper mode to be prominent in the beginning of all our
runs. Later on its significance tends to go down in simulations with
${\mathcal{M}}\lesssim 9$, whereas in simulations with higher ${\mathcal{M}}$
it may reappear later on, but not always: the upper mode is absent in our
${\mathcal{M}}=12$ runs but persists through the whole duration of the
simulation in ${\mathcal{M}}=15$ case, see Figures 7 & 9.

The lower mode is seen in most of our runs, often through the full simulation
duration, like in ${\mathcal{M}}=12$ case (but we remind that
${\mathcal{M}}=12$ runs are quite unique in maintaining extremely stable lower
mode, see §5.3). They are far less prominent in ${\mathcal{M}}\geq 13$ runs,
but are still present there at some level, see below.

The general expectation following from the theory of acoustic mode excitation
outlined in Belyaev & Rafikov (2012) and Belyaev et al. (2013a) is that the
modes should have comparable strength (in $rv_{r}\sqrt{\Sigma}$) immediately
inside and outside the star. However, very often it is much easier to detect a
particular mode inside the star than outside. For example, ${\mathcal{M}}=6$
run (Figure 5) at $t/2\pi>300$ shows a telltale $k_{r}=0$ (i.e. radially
elongated perturbation pattern) signature of the lower $m=8,9,10,12$ modes
inside the star, which do not have a counterpart with matching pattern
frequencies outside the star, see Figure 6. We speculate that this departure
from the theoretical expectation can be at least partly caused by the non-
uniform surface density distribution in the inner disk.

In other cases the apparent lack of the disk counterpart for a mode may be
caused by its overlap with some other modes, complicating its identification.
This is likely the case for the upper $m=19$ mode with $\Omega_{\rm p}\approx
0.65$ in the ${\mathcal{M}}=15$ run shown in Figure 9: this mode is obvious
inside the star (note its non-zero $k_{r}$ there), whereas it is hardly
visible in the disk. However, Figure 10d shows that this mode is in fact also
present in the disk (at $r=1.2$, albeit with a substantially reduced
amplitude) with the same $\Omega_{\rm p}$; it is hard to detect by eye in
simulation snapshots because of its interference with other modes in the disk.
This comparison demonstrates the benefits of automatic mode detection
procedure that we employ in analyzing our simulations.

We now examine how the mix of modes detected by our automated analysis
compares with the dispersion relations derived in Belyaev et al. (2013a) and
this work, see §3.1,3.2.

####  8.1.1 Dispersion relation for acoustic modes

Figure 14 displays the $(m,\Omega_{\rm p})$ pairs for all the modes found by
our automatic mode detection procedure in simulations with different values of
${\mathcal{M}}$. Some of these modes are truly global (green stars), i.e. they
are detected as a wave pattern with the same $m$ and $\Omega_{\rm p}$ at the
same interval of time both inside the star and in the disk (at $r=1.2$). In
most cases modes are found only in the disk (yellow pluses) or only in the
star (blue circles), a possibility that we mentioned earlier.

We also display in red the dispersion relations (9) for the lower modes (dot-
dashed) and (23) for the upper modes (dotted). Note that equation (9) depends
on a parameter $r_{0}$ (specific for each ${\mathcal{M}}$), which we fix by
aligning the lower mode dispersion relation curves with the clusters of
$(m,\Omega_{\rm p})$ points in Figure 14; this procedure is not very
straightforward for ${\mathcal{M}}=14,15$, see panels (k) and (l). We also
note that the dispersion relation (23) assumes that a plateau in $\Omega(r)$
has already developed near the stellar surface (see Figure 4a), so that the
epicyclic frequency is $\kappa\approx 2\Omega$; this may not be true early on
in the simulation. The “height" of the plateau at late times $\Omega_{\rm
max}$, i.e. the maximum value of $\Omega$, is shown by the horizontal dashed
curves.

In general we see good correspondence between the dispersion relations and the
detected modes, as typically a significant number of $(m,\Omega_{\rm p})$
points fall on top of the red curves. These modes often cover a significant
range of azimuthal wavenumbers (e.g. the lower modes in panels (f) or (h)),
and sometimes come in clusters, i.e. are grouped in $m$ (e.g. $m=10-16$ in
panel (e), or $m=21-26$ in panel (g)). Such groupings may result from the
nonlinear evolution of a single dominant mode: nonlinear distortion of a
perturbation profile, natural for even moderately strong acoustic waves that
we observe in our runs, transfers power to other modes, primarily to those
with similar $m$. This likely explains the slow but persistent changes in the
modes that we see in our runs: as the acoustic waves in the disk are
dispersive, the new modes produced as a result of the non-linear evolution of
a parent mode eventually lose coherence with it, smearing out the original
wave packet. Thus, the finite lifetime of the dominant modes that we see in
our runs should not come as a surprise.

At the same time, there are some modes that do not line up with the red
curves. Some of these modes are simply not the usual upper and lower acoustic
modes, see e.g. §8.1.1. But many of the outliers are the higher-order
azimuthal harmonics of the primary modes. To illustrate that we show the
harmonics of the main dispersion relation branches with twice and three times
higher azimuthal wavenumber $m$ and the same $\Omega_{\rm p}$ as black dotted
and dot-dashed curves in Figure 14. One can see that in some cases the modes
lying on the main branch of the dispersion relation have counterparts with the
same $\Omega_{\rm p}$ at or close to one of the higher-order branches of that
dispersion relation. Clear examples of this can be seen in panel (b) for the
lower modes with $\Omega_{\rm p}\in(0.45,0.55)$, and in panel (g) for the
upper modes with $\Omega_{\rm p}\in(0.5,0.6)$. Such higher-order azimuthal
counterparts of the modes naturally arise from the non-sinusoidal shape of the
wave packets with certain azimuthal periodicity.

For almost all values of ${\mathcal{M}}$ we also see some modes that stay
close to the $\Omega_{\rm p}=\Omega_{\rm max}$ horizontal dashed line. These
modes are likely trapped in the innermost part of the disk where at late times
$\Omega(r)$ features a plateau with $\Omega(r)\approx\Omega_{\rm max}$. These
modes are likely related to the vortensity structures forming in this part of
the disk — vortices or rolls — as the latter are passively advected with the
fluid at almost constant orbital frequency $\Omega_{\rm max}$. Stability of
$\Omega_{\rm max}$ (see Figure 4) should help these modes maintain their
coherence over long intervals of time, which may have important implications
for the variability associated with the BL (such as dwarf nova oscillations);
this issue will be explored in future work. Also, for ${\mathcal{M}}<8$ some
lower modes feature $\Omega_{\rm p}$ exceeding $\Omega_{\rm max}$; these modes
must have been present in the disk early on, when the $\Omega(r)$ profile was
still close to Keplerian.

Finally, panels (e) and (f) of Figure 14 compare two simulations with the same
Mach number ${\mathcal{M}}=9$ but different resolutions. The higher resolution
run M09.HR.r.a ($8192\times 8192$) appears to show no disk modes, in contrast
with the run M09.FR.r.a ($4096\times 4096$), which reveals a number of both
global and disk-only modes. However, this outcome is caused simply by the
difficulty of mode detection by our automated mode-finding algorithm in the
high-resolution simulation: by examining its outputs by eye we do find a
number of disk modes. But their lifetime is likely too short for them to be
registered as waves with a well-defined $\Omega_{\rm p}$ by our software.

###  8.2 Other disk-only modes

Next we briefly discuss a couple of other wave structures that occasionally
emerge in our simulations and cannot be classified as upper or lower modes (or
their harmonics). Both of these modes are present only in the disk, with no
counterpart inside the star.

####  8.2.1 Resonant modes

One type of such modes are the relatively low-$m$ modes in the disk trapped in
the resonant cavity near the star; for this reason these modes may be confused
with the usual lower modes. However, unlike the lower modes they (1) do not
have a strong counterpart with the same $m$ inside the star, (2) usually do
not exhibit densely criss-crossed pattern, and (3) obey a very different
dispersion relation, as we show next. The difference in appearance between the
two types of modes can be easily spotted in Figure 7 for ${\mathcal{M}}=12$,
where a strong lower mode is present in panels with $t/2\pi=150-450$, whereas
at $t/2\pi=550$ there is a strong $m=7$ disk-only mode with no crossings of
the incoming and outgoing wakes, trapped at $r<1.25$.

Figure 15: Dispersion relation (21) for the resonant modes for $q=2$, plotted
together with the $(m{\mathcal{M}},\Omega_{\rm p})$ values for the resonant
modes identified in our simulations, shown with points colored by their value
of ${\mathcal{M}}$. The bright magenta points correspond to the
${\mathcal{M}}=8$ simulation, which show the largest deviation from the trend.
See §8.2.1 for details.

This type of mode manifests itself also at other values of ${\mathcal{M}}$: as
a strong $m=2$ mode for ${\mathcal{M}}=6$ (at $t/2\pi=325-425$), confined to
$r<2.5$; as a $m=6$ mode for ${\mathcal{M}}=9$ (at $t/2\pi=400-450$), confined
to $r<1.3$ (although strongly disturbed by the $m=5$ vortex-driven mode); and
a $m\approx 20$ mode at ${\mathcal{M}}=15$ (at $t/2\pi=275-5756$), confined to
$r<1.17$. We also find this mode type to persist in our very long
${\mathcal{M}}=9$ run, where it strongly manifests itself as either $m=3$ or
$m=4$ pattern for more than 2000 orbits.

Such modes were previously described in Belyaev et al. (2012), who traced
their origin to a geometric resonance for a trapped acoustic wave: if, after
multiple reflections off the stellar surface and the Inner Lindblad Resonance,
the density wave closes on itself (after its azimuthal phase wraps around the
star $q$ times, where $q$ is a small integer), then this reinforces its
strength and gives rise to a stable mode. Mathematically, Belyaev et al.
(2012) have shown that this leads to a relationship between $m$ and
$\Omega_{\rm p}$ for these modes, which can be cast as

$\displaystyle
q\frac{\pi}{{\mathcal{M}}}=\frac{1}{\Omega(R_{\star})R_{\star}}\int_{R_{\star}}^{r_{\rm
ILR}}dr\sqrt{m^{2}\left[\Omega(r)-\Omega_{\rm p}\right]^{2}-\kappa^{2}(r)},$
(20)

where $r_{\rm ILR}$ is given by equation (12). It was also shown in that work
that such resonant modes indeed obey the dispersion relation (20), see their
Figure 18.

One could turn the integral relationship between $m$ and $\Omega_{\rm p}$ in
equation (20) into an approximate algebraic one by dropping the $\kappa^{2}$
term. This is equivalent to approximating $r_{\rm ILR}\approx r_{c}$ and is
accurate for $m\gg 1$, see equation (12) and Figure 1b. Then one finds

$\displaystyle m{\mathcal{M}}\approx q\pi\left[2+\frac{\Omega_{\rm
p}}{\Omega_{K}(R_{\star})}-3\left(\frac{\Omega_{\rm
p}}{\Omega_{K}(R_{\star})}\right)^{1/3}\right]^{-1}.$ (21)

As ${\mathcal{M}}$ in our runs increases, we find the resonant mode wavenumber
$m$ to increase as well. As a result, equation (21) predicts that $\Omega_{\rm
p}$ for this mode should also increase with both ${\mathcal{M}}$ and its $m$,
see Figures 17 and 18 of Belyaev et al. (2012). This leads to narrowing of the
resonant cavity for this mode as ${\mathcal{M}}$ goes up.

One can see that $m$ and ${\mathcal{M}}$ enter equation (21) only in
combination $m{\mathcal{M}}$. This allows us to plot the dispersion relation
(21) as a single curve for the different values of ${\mathcal{M}}$ in
$(m{\mathcal{M}},\Omega_{\rm p})$ coordinates. We do this in Figure 15, where
we also plot $(m{\mathcal{M}},\Omega_{\rm p})$ points for all resonant modes
that we were able to reliably identify in our simulations. One can see that
with $q=2$ the dispersion relation (21) fits the simulation results quite
well. The only exception are the two occurrences of the resonant mode in our
${\mathcal{M}}=8$ run, for which a different value of $q$ might have worked
better as we see multiple crossings of resonant modes in this run (usually we
see only a single crossing). Note that $q=2$ that we find in this work is
different from $q=1$ found in Belyaev et al. (2012), not clear why.

The dispersion relation shown in Figure 15 is clearly different from that of
the lower modes, for which $\Omega_{\rm p}$ always decreases with $m$, see
Figure 14. This is despite the fact that the two types of modes have similar
morphological appearance and are confined to a resonant cavity in the disk;
they also have a similar effect on the angular momentum and mass transport in
the disk (Coleman et al., in prep.).

At the same time, the dispersion relation (7) for the upper acoustic mode (as
well as its more refined version (23)) leads to $\Omega_{\rm p}$ increasing
with $m$, similar to the behavior predicted by the equation (21). For that
reason, in Figure 14 we often find the $(m,\Omega_{\rm p})$ pairs for the
resonant modes to lie close to the main branch of the upper mode dispersion
relation (e.g. see $m=2$, $\Omega_{\rm p}\approx 0.15$ resonant mode for
${\mathcal{M}}=6$ in Figure 14b, or $m=6$, $\Omega_{\rm p}\approx 0.5$
resonant mode for ${\mathcal{M}}=9$ in Figure 14e).

####  8.2.2 Low-$m$, $k_{r}=0$ modes

As mentioned in §5.1, our fiducial ${\mathcal{M}}=9$ run shows yet another
disk-only mode with $m=2$, $\Omega_{\rm p}\approx 0.15$, between roughly 150
and 400 orbits in Figure 2. It has a very unusual appearance, with $k_{r}=0$
and azimuthally extended perturbation pattern (i.e. not a narrow feature),
confined to $r\lesssim 2.2$, which is close to the $r_{\rm ILR}$ for this
mode. Its perturbation also undergoes a flip by $\pi$ in azimuthal phase at
$r\approx 1.5$. Vortensity maps in Figure 2 show no structures in $\zeta$ at
this radius or beyond it.

The emergence of this mode is not unique to the run M09.FR.r.a displayed in
Figure 2, as we observe it in several other ${\mathcal{M}}=9$ runs with
different kinds of initial perturbation. The low $m$ and $\Omega_{\rm p}$ of
this mode places it very close to the main branch of the upper mode dispersion
relation (23), see Figure 14e. This is not surprising since that dispersion
relation was derived assuming $k_{r}=0$ (see §A), which is true at all $r$ for
the $m=2$ mode that we see. At the moment we do not have an explanation for
the origin or properties of this unusual disk-only mode.

###  8.3 Dominant modes as a function of ${\mathcal{M}}$

Given that we have BL simulations for every integer value of ${\mathcal{M}}$
between 5 and 15, we can explore how the harmonic content of the modes that we
detect changes with ${\mathcal{M}}$.

In general, we find that a particular mix of modes that exist at different
times during a given simulation is pretty stochastic. This means that a
different realization of the same simulation, especially with the different
model of initial noise introduced to trigger the instability (§4.2), would
result in a somewhat different outcome in terms of the mode types and
azimuthal wavenumbers. Resolution of the simulations also plays a role, see
Figure 14e,f, but the differences there often depend on the performance of our
mode-finding analysis software. The only notable exception are our
${\mathcal{M}}=12$ simulations, in all of which we robustly see the $m=16$
lower mode with a pattern speed $\Omega_{\rm p}=0.45$ dominating both inside
and outside of the star over hundreds of orbits.

Figure 16: Histogram of the azimuthal wave number $m$ for the dominant modes
identified in simulations with different values of Mach number $\mathcal{M}$.
For each $\mathcal{M}$ there are three histograms with arbitrary horizontal
displacement. Left (red): dominant star-only modes. Middle (blue): dominant
disk-only modes. Right (grey): global modes. For each ${\mathcal{M}}$ the
dominant modes are defined as the three modes with the highest time-integrated
power, volume averaged over a specific region. For some values of
$\mathcal{M}$ we have multiple simulations, so the numeric value of the
corresponding bin reflects the number of runs, in which the mode meets the
aforementioned criteria. The dotted line corresponds to $\mathcal{M}=m$ which
qualitatively follows the trend seen in the data. See §8.3 for details.

On the other hand, we do observe certain trends with ${\mathcal{M}}$. In
particular, by carefully examining Figure 14 one can notice that the numbers
of $(m,\Omega_{\rm p})$ points aligned with the branches for either the upper
or lower acoustic modes change more or less monotonically as ${\mathcal{M}}$
is varied. More specifically, lower modes are more common for
${\mathcal{M}}\leq 12$, whereas the upper modes start showing up in noticeable
clusters along the upper dispersion relation branches for ${\mathcal{M}}\geq
9$.

To examine possible trends with ${\mathcal{M}}$ at a more quantitative level,
we came up with the following procedure. First, we compute the power in all
modes with $m<32$ for the fluid variable $v_{r}\sqrt{\Sigma}$ and then
integrate it over time for $t/2\pi>100$ and over radius in three distinct
domains: “star" defined as $r<1$, “disk" defined as $1<r<2.2$, and “global"
defined as $r<2.2$, i.e. “disk+star". In a given domain, the three modes with
the highest time- and radius- integrated power are considered to be the
dominant mode. This data is summarized in Fig. 16, where the histograms of
different color characterize the distribution of $m$ for the dominant modes in
three respective regions for all ${\mathcal{M}}$.

By examining this figure we find that at each $\mathcal{M}$ there is a
substantial spread in the values of $m$, even for a given domain. This spread
is caused by a number of factors: stochasticity of the mix of modes, different
types of modes involved (e.g. upper, lower, disk-only), resolution of
different simulations, etc. Also, we have reasonably representative statistics
on the distribution of $m$ only for ${\mathcal{M}}=6,9,12$, for which there
are multiple runs with different initial conditions and resolutions; for most
other values of ${\mathcal{M}}$ only a single run is available.

Despite these issues, we can still make some observations. In particular,
inside the star the dominant mode’s $m$ typically obeys $3/4\lesssim
m/\mathcal{M}\lesssim 3/2$, although we have seen $m$ as high as $\sim
2.5\mathcal{M}$ for the fiducial and high resolution simulations. Low
resolution simulations typically do not show mode numbers this high. Dominant
disk modes tend to have lower $m$, most likely because of the vortex-driven
modes at intermediate values of ${\mathcal{M}}$ and the low-$m$ resonant
modes. Dominant global modes tend to have $m$ somewhat higher than the modes
inside the star.

Qualitatively, there is an overall trend of increasing the dominant mode wave
number $m$ with growing Mach number $\mathcal{M}$. Just as a guide, dotted
line in Figure 16 shows a linear relation $m={\mathcal{M}}$. This line does
not represent a fit of any kind and is merely shown to guide the eye. One
conclusion that we can draw from this exercise is that a complete
characterization of the mix of the dominant modes operating in the vicinity of
the BL may require a substantially larger number of simulations than we have
presented in this work.

###  8.4 Sensitivity to numerical parameters

Our simulation suite allows us to probe the sensitivity of the results to
certain numerical inputs for some values of ${\mathcal{M}}$, namely the
initial noise pattern used to trigger the sonic instability in the BL (§4.2)
and resolution (§4.3), see Table 1.

When comparing the simulations with the different forms of the initial noise
(run at the same resolution and ${\mathcal{M}}$), we generally do not find
strong differences or trends with the noise pattern. For ${\mathcal{M}}=6,9$
the qualitative behavior of the simulations remains the same, although as we
alluded to in §5.1.1,8.3, the detailed outcomes of individual simulations are
stochastic. And all ${\mathcal{M}}=12$ runs are similar to one another even at
the quantitative level, see §5.3.

Resolution has a more substantial effect on our results. For ${\mathcal{M}}=9$
it affects dominant azimuthal wavenumber $m$ of the modes that we observe,
with $m$ increasing with resolution. For example, in simulations with the same
block random initial condition “r.a", we find that the dominant lower mode has
$m=7$ at lowest resolution ($2048\times 2048$), increasing to $m=19$ in the
fiducial resolution ($4096\times 4096$) case, and reaches $m=26$ at the
highest resolution ($8192\times 8192$). The transitions between the different
types of modes described in §5.1 are captured quite reliably between the high
and fiducial resolution cases, suggesting that their results are converged at
least at the qualitative level, but less so at the lowest resolution.

The dependence on resolution is stronger in the ${\mathcal{M}}=6$ simulations.
In particular, high resolution ($2048\times 2048$) ${\mathcal{M}}=6$ runs
demonstrate the early development of the low-$m$ ($m=2,3$) resonant modes,
typically around $t/2\pi=200$. At the same time, the low resolution
($1024\times 1024$) runs either do not show these modes at all for the
duration of the simulation, or exhibit them very late. Thus, high resolution
is clearly necessary for revealing important features of the BLs with low
${\mathcal{M}}$.

###  8.5 Comparison with the existing studies

There is a number of past studies on the BL structure, both (semi-)analytical
(Kippenhahn & Thomas, 1978; Narayan & Popham, 1993; Popham et al., 1993) and
numerical (Kley & Lin, 1996; Armitage, 2002; Steinacker & Papaloizou, 2002;
Balsara et al., 2009; Hertfelder & Kley, 2017), that postulate some form of
local shear stress to enable angular momentum transport inside the BL. Since
in practice transport in the BL is mediated by the global acoustic modes
(Belyaev et al., 2012, 2013a), these studies cannot be directly compared to
our work.

Our study goes beyond (in ways already discussed in §4.4) the similar past
work of Belyaev et al. (2012, 2013a, 2013b) and Hertfelder & Kley (2015b), who
also simulated BLs mediated by the acoustic waves. We explore a larger set of
Mach number values, use higher resolution and longer run times, and carry out
an extensive exploration of the sensitivity of our results to resolution and
initial conditions. We also provide a very careful analysis of our results and
extensively study the harmonic content of our simulations. All this led to new
important findings such as the vortex-driven modes (§6), one-armed spirals
(6.4), and so on.

Belyaev (2017) considered a different way of exciting acoustic modes in the
disk, namely by coupling them to the incompressible inertial waves inside the
star. Our use of the globally isothermal equation of state precludes us from
exploring this possibility, which should be addressed by future simulations
with more sophisticated treatment of gas thermodynamics. Finally, in their 3D,
unstratified MHD simulations Belyaev & Quataert (2018) found that acoustic
waves do not efficiently remove angular momentum from the accreting material
in the BL, causing a dense belt of rapidly spinning material to form near the
stellar equator. While this is an important issue, which should be addressed
in the future using stratified MHD simulations with realistic thermodynamics,
Belyaev & Quataert (2018) do find acoustic waves to be active in the disk,
which is what our study has focused on.

###  8.6 Observational implications

Observational implications of the wave-driven angular momentum transport in
the BL have been previously discussed in Belyaev et al. (2012, 2013a). One of
them is the modification of the spectral signature associated with the energy
dissipation in the BL. While the energy conservation implies that the total
amount of energy released by the accreted matter must be large, a particular
band in which the associated emission is released should be dramatically
affected by the global nature of the angular momentum and energy transport by
the acoustic modes. This is likely to have important ramifications for the so-
called “missing boundary layer" problem (Ferland et al., 1982).

Long-lived mode patterns that we see in our simulations should also lead to
characteristic variability associated with them. Our current work opens up new
possibilities in this regard, by revealing the existence of the vortex-driven
modes and one-armed spiral arms in the disk. Both of them may affect the light
curves of objects accreting through the BL in characteristic ways.

These (and other) implications of the wave-mediated accretion through the BL
will be further investigated in the upcoming work.

##  9 Summary

In this work we explored the physics of the BLs (forming at the surfaces of
weakly magnetized accreting objects), which are mediated by the instabilities
excited in a supersonic shear layer at the inner edge of the disk. In this
paper, first in a series, we presented a suite of global, 2D, hydrodynamic
simulations of the BLs using simple thermodynamics (globally isothermal) and
encompassing both the outermost layers of the accreting object and a
substantial region of the inner disk. Using this data set we carried out a
systematic exploration of the different waves that emerge in the disk in the
vicinity of the BL (and regulate its properties), as a function of Mach number
of the system ${\mathcal{M}}$ (for which we explored a range of values, from
${\mathcal{M}}=5$ to ${\mathcal{M}}=15$). Our results on the morphological
characteristics of these modes can be summarized as follows.

  * •

Our simulations reveal a plethora of different modes operating in the vicinity
of the BL, some of which have been known previously, e.g. upper and lower
acoustic modes (Belyaev et al., 2013a), whereas others are identified for the
first time. Different types of modes can easily co-exist in the inner disk.
They have finite life times, although some of them can operate for hundreds of
orbits.

  * •

We discover a new class of modes that are triggered by the emergence of non-
axisymmetric, localized vortensity structures in the vicinity of the BL. These
vortex-driven modes are quite prominent in many of our simulations. We argue
that their parent vortices result from the Rossby wave instibility, triggered
by the vortensity production near the BL, which is caused by the nonlinear
damping of the acoustic modes.

  * •

In many simulations we observe multiple near-BL vortices to merge into a
single one, giving rise to a prominent global, one-armed spiral density wave
in the disk. Structures like this may cause periodicity of emission associated
with the BL accretion.

  * •

While some of the modes that we find are global, i.e. operate both in the disk
and the star, others are star-only or disk-only modes. We analyzed the
properties of several disk-only modes trapped in the vicinity of the BL.

  * •

We compared characteristics of many different modes identified in our
simulations with their analytical dispersion relations and found good
agreement.

  * •

As the Mach number of our simulations changes we find the mix of near-BL modes
to evolve, although not always in a monotonic fashion. Nevertheless, we find
indications that the azimuthal wavenumber of the dominant modes emerging in
our runs grows with ${\mathcal{M}}$ roughly linearly.

  * •

We find moderate sensitivity of our simulation outcomes to the initial
conditions, i.e. the form of the seed perturbation triggering the instability
in the BL. With rare exceptions, a particular mix of modes that we find in our
runs is somewhat stochastic. High resolution is important for the accurate
numerical description of the BL physics.

Our results pave the way for future efforts to explore angular momentum and
mass transport, as well as the associated evolution of the disk in the
vicinity of the BL, driven by the modes identified in this study.

## Acknowledgements

We thank Changgoo Kim for his assistance in modifying the FFT module in
Athena++ and Jim Stone for making the code publicly available. We gratefully
acknowledge financial support from NSF via grant AST-1515763, NASA via grant
14-ATP14-0059, and Institute for Advanced Study via the John N. Bahcall
Fellowship to R.R.R. Research at the Flatiron Institute is supported by the
Simons Foundation. Resources supporting this work were provided by the NASA
High-End Computing (HEC) Program through the NASA Advanced Supercomputing
(NAS) Division at Ames Research Center. Through allocation AST160008, this
work used the Extreme Science and Engineering Discovery Environment (XSEDE),
which is supported by National Science Foundation grant number ACI-1548562
(Towns et al., 2014).

## References

  * Armitage (2002) Armitage P. J., 2002, MNRAS, 330, 895
  * Bae & Zhu (2018) Bae J., Zhu Z., 2018, ApJ, 859, 118
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Balsara et al. (2009) Balsara D. S., Fisker J. L., Godon P., Sion E. M., 2009, ApJ, 702, 1536
  * Barranco & Marcus (2005) Barranco J. A., Marcus P. S., 2005, ApJ, 623, 1157
  * Belyaev (2017) Belyaev M. A., 2017, ApJ, 835, 238
  * Belyaev & Quataert (2017) Belyaev M. A., Quataert E., 2017, preprint,  (arXiv:1709.01197)
  * Belyaev & Quataert (2018) Belyaev M. A., Quataert E., 2018, MNRAS, 479, 1528
  * Belyaev & Rafikov (2012) Belyaev M. A., Rafikov R. R., 2012, ApJ, 752, 115
  * Belyaev et al. (2012) Belyaev M. A., Rafikov R. R., Stone J. M., 2012, ApJ, 760, 22 (BRS)
  * Belyaev et al. (2013a) Belyaev M. A., Rafikov R. R., Stone J. M., 2013a, ApJ, 770, 67 (BRS1)
  * Belyaev et al. (2013b) Belyaev M. A., Rafikov R. R., Stone J. M., 2013b, ApJ, 770, 68 (BRS2)
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition
  * Chandrasekhar (1960) Chandrasekhar S., 1960, Proceedings of the National Academy of Science, 46, 253
  * Coleman et al. (2016) Coleman M. S. B., Kotko I., Blaes O., Lasota J.-P., Hirose S., 2016, MNRAS, 462, 3710
  * Coleman et al. (2018) Coleman M. S. B., Blaes O., Hirose S., Hauschildt P. H., 2018, ApJ, 857
  * Dong et al. (2011) Dong R., Rafikov R. R., Stone J. M., 2011, ApJ, 741, 57
  * Ferland et al. (1982) Ferland G. J., Langer S. H., MacDonald J., Pepper G. H., Shaviv G., Truran J. W., 1982, ApJ, 262, L53
  * Ghosh et al. (1977) Ghosh P., Lamb F. K., Pethick C. J., 1977, ApJ, 217, 578
  * Gilfanov et al. (2003) Gilfanov M., Revnivtsev M., Molkov S., 2003, A&A, 410, 217
  * Goldreich & Tremaine (1980) Goldreich P., Tremaine S., 1980, ApJ, 241, 425
  * Heinemann & Papaloizou (2009) Heinemann T., Papaloizou J. C. B., 2009, MNRAS, 397, 52
  * Hertfelder & Kley (2015a) Hertfelder M., Kley W., 2015a, A&A, 579, A54
  * Hertfelder & Kley (2015b) Hertfelder M., Kley W., 2015b, A&A, 579, A54
  * Hertfelder & Kley (2017) Hertfelder M., Kley W., 2017, A&A, 605, A24
  * Hertfelder et al. (2013) Hertfelder M., Kley W., Suleimanov V., Werner K., 2013, A&A, 560, A56
  * Hirose (2015) Hirose S., 2015, MNRAS, 448, 3105
  * Hirose et al. (2014) Hirose S., Blaes O., Krolik J. H., Coleman M. S. B., Sano T., 2014, ApJ, 787, 1
  * Inogamov & Sunyaev (1999) Inogamov N. A., Sunyaev R. A., 1999, Astronomy Letters, 25, 269
  * Inogamov & Sunyaev (2010) Inogamov N. A., Sunyaev R. A., 2010, Astronomy Letters, 36, 848
  * Johnson & Gammie (2005) Johnson B. M., Gammie C. F., 2005, ApJ, 635, 149
  * Kevlahan (1997) Kevlahan N. K. R., 1997, Journal of Fluid Mechanics, 341, 371
  * Kippenhahn & Thomas (1978) Kippenhahn R., Thomas H. C., 1978, A&A, 63, 265
  * Kley & Lin (1996) Kley W., Lin D. N. C., 1996, ApJ, 461, 933
  * Koenigl (1991) Koenigl A., 1991, ApJ, 370, L39
  * Koller et al. (2003) Koller J., Li H., Lin D. N. C., 2003, ApJ, 596, L91
  * Lesur & Papaloizou (2009) Lesur G., Papaloizou J. C. B., 2009, A&A, 498, 1
  * Li et al. (2000) Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ, 533, 1023
  * Li et al. (2001) Li H., Colgate S. A., Wendroff B., Liska R., 2001, ApJ, 551, 874
  * Li et al. (2005) Li H., Li S., Koller J., Wendroff B. B., Liska R., Orban C. M., Liang E. P. T., Lin D. N. C., 2005, ApJ, 624, 1003
  * Lin (2012) Lin M.-K., 2012, MNRAS, 426, 3211
  * Lin (2015) Lin M.-K., 2015, MNRAS, 448, 3806
  * Lin & Papaloizou (2010) Lin M.-K., Papaloizou J. C. B., 2010, MNRAS, 405, 1473
  * Lin & Pierens (2018) Lin M.-K., Pierens A., 2018, MNRAS, 478, 575
  * Lithwick (2009) Lithwick Y., 2009, ApJ, 693, 85
  * Lovelace et al. (1999) Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805
  * Meheut et al. (2012a) Meheut H., Yu C., Lai D., 2012a, MNRAS, 422, 2399
  * Meheut et al. (2012b) Meheut H., Keppens R., Casse F., Benz W., 2012b, A&A, 542, A9
  * Miranda & Rafikov (2019a) Miranda R., Rafikov R. R., 2019a, ApJ, 875, 37
  * Miranda & Rafikov (2019b) Miranda R., Rafikov R. R., 2019b, The Astrophysical Journal, 878, L9
  * Miranda & Rafikov (2020) Miranda R., Rafikov R. R., 2020, ApJ, 892, 65
  * Narayan & Popham (1993) Narayan R., Popham R., 1993, Nature, 362, 820
  * Ogilvie & Lubow (2002) Ogilvie G. I., Lubow S. H., 2002, MNRAS, 330, 950
  * Ono et al. (2016) Ono T., Muto T., Takeuchi T., Nomura H., 2016, ApJ, 823, 84
  * Ono et al. (2018) Ono T., Muto T., Tomida K., Zhu Z., 2018, ApJ, 864, 70
  * Paardekooper et al. (2010) Paardekooper S.-J., Lesur G., Papaloizou J. C. B., 2010, ApJ, 725, 146
  * Papaloizou & Lin (1989) Papaloizou J. C. B., Lin D. N. C., 1989, ApJ, 344, 645
  * Pessah & Chan (2012) Pessah M. E., Chan C.-k., 2012, ApJ, 751, 48
  * Philippov et al. (2016) Philippov A. A., Rafikov R. R., Stone J. M., 2016, ApJ, 817, 62
  * Popham et al. (1993) Popham R., Narayan R., Hartmann L., Kenyon S., 1993, ApJ, 415, L127
  * Rafikov (2002) Rafikov R. R., 2002, ApJ, 569, 997