# Boundary Layers of Accretion Disks: Wave-Driven Transport and Disk Evolution

Matthew S. B. Coleman,1,2 Roman 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 E-mail<EMAIL_ADDRESS>N. Bahcall Fellow
at the IAS

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

###### Abstract

Astrophysical objects possessing a material surface (white dwarfs, neutron
stars, etc.) accrete gas from the disk through the so-called surface boundary
layer (BL), in which the angular velocity of the accreting gas experiences a
sharp drop. Acoustic waves excited by the supersonic shear in the BL play an
important role in mediating the angular momentum and mass transport through
that region. Here we examine the characteristics of the angular momentum
transport produced by the different types of wave modes emerging in the inner
disk, using the results of a large suite of hydrodynamic simulations of the
BLs. We provide a comparative analysis of the transport properties of
different modes across the range of relevant disk parameters. In particular,
we identify the types of modes which are responsible for the mass accretion
onto the central object. We find the correlated perturbations of surface
density and radial velocity to provide an important contribution to the mass
accretion rate. Although the wave-driven transport is intrinsically non-local,
we do observe a clear correlation between the angular momentum flux injected
into the disk by the waves and the mass accretion rate through the BL. We find
the efficiency of angular momentum transport (normalized by thermal pressure)
to be a weak function of the flow Mach number. We also quantify the wave-
driven evolution of the inner disk, in particular the modification of the
angular frequency profile in the disk. Our results pave the way for
understanding wave-mediated transport in future three-dimensional,
magnetohydrodynamic studies of the BLs.

###### keywords:

accretion, accretion discs – hydrodynamics – instabilities

††pubyear: 2015††pagerange: Boundary Layers of Accretion Disks: Wave-Driven
Transport and Disk Evolution–A

##  1 Introduction

Disk accretion onto an object possessing a material surface (i.e. not a black
hole) is an important problem emerging in a variety of astrophysical settings.
When the magnetic field of an accretor is weak enough (Ghosh & Lamb, 1978;
Koenigl, 1991), the accretion flow can extend all the way to its surface. In
this case the transfer of incoming matter onto the central object (hereafter
"a star") must proceed through the so-called boundary layer (hereafter BL) — a
narrow region between the disk and the accretor, in which the angular velocity
of the accreting matter adjusts to the rotation of the central star. In order
for that to happen, the gas arriving from the disk must lose its angular
momentum inside the BL.

The magneto-rotational instability (MRI; Velikhov 1959; Balbus & Hawley 1991),
traditionally invoked as the favored angular momentum transport mechanism in
sufficiently ionized accretion disks, cannot enable this process in the BL: it
does not operate in this region (Pessah & Chan, 2012) since the angular
velocity of the fluid there necessarily increases with the distance. Instead,
Belyaev & Rafikov (2012) proposed that a sonic instability, driving the
excitation of acoustic waves by the supersonic shear flow inside the BL,
mediates the angular momentum transport in that region in a non-local fashion,
very different from the local $\alpha$-models Shakura & Sunyaev (1973), that
have been previously invoked in the BL context (Popham et al., 1993; Balsara
et al., 2009; Hertfelder & Kley, 2017). This instability has been subsequently
verified to robustly operate within the BL and drive angular momentum
transport using hydrodynamic and magneto-hydrodynamic (MHD) simulations
(Belyaev et al., 2012, 2013a, 2013b).

Recently Coleman et al. (2021) (hereafter Paper I) presented a new, extensive
suite of two-dimensional (2D) hydrodynamic simulations of the BLs run for
multiple values of the Mach number ${\mathcal{M}}=v_{K}(R_{\star})/c_{s}$ —
the ratio of the Keplerian velocity $v_{K}(R_{\star})$ at the inner edge of
the disk (i.e. at the stellar radius $R_{\star}$) to the sound speed $c_{s}$
(which was assumed constant in that work). They carefully analyzed the
different waves emerging in the inner disk as a result of the sonic
instability operating in the BL, and discovered some new types of modes
emerging in the vicinity of the BL, e.g. the vortex-driven modes. That study
focused predominantly on the morphological characteristics of the waves and
their correspondence to the known analytical results.

In this work we will use the simulation suite presented in Paper I to
characterize the BL from a different angle, namely to explore the angular
momentum and mass transport in its vicinity driven by the wave activity. A
number of past studies explored the connection between the two transport
processes in accretion disks. For example, Balbus & Papaloizou (1999) argued
that the MRI-driven transport can be characterized as a local
$\alpha$-viscosity (Shakura & Sunyaev, 1973), whereas the transport driven by
disk self-gravity is global and cannot be represented using the
$\alpha$-anzatz (see also Larson 1984; Gnedin et al. 1995). A number of
studies have also looked at the angular momentum transport by the global
spiral waves in disks (Larson, 1990; Spruit, 1987; Rafikov, 2016; Arzamasskiy
& Rafikov, 2018), in particular those driven by embedded planets Goodman &
Rafikov (2001); Rafikov (2002). In the BL context, Belyaev et al. (2012,
2013a) explored some global transport characteristics of the acoustic waves,
whereas Belyaev et al. (2013b) did the same in the MHD case111In their global,
unstratified MHD simulations of the near-BL region Belyaev & Quataert (2018)
found accumulation of accreted material in a belt-like structure in the BL,
spinning at sub-Keplerian velocity. To allow the accreted material to join the
star this belt must dissipate somehow, but it is not yet clear how that
happens. Nevertheless, even with the belt Belyaev & Quataert (2018) still
found a substantial wave activity in the inner disk, which is what matters for
us in this study.. Also, Dittmann (2021) studied how the efficiency of
transport varies as the spin of the accreting object changes.

The aim of our present work is to extend the latter studies by looking at the
various characteristics of the wave-driven transport in the vicinity of the BL
— angular momentum and mass fluxes — in a systematic fashion across the range
of ${\mathcal{M}}$ values. A distinct feature of this work compared to
previous studies is that we examine the transport properties of individual
modes emerging in our simulations and determine their variation across the
different types of modes. Another goal is to explore the wave-driven evolution
of the inner parts of the accretion disk adjacent to the BL.

This work is organized as follows. After briefly describing in §2 the
simulation suite on which this study is based, we remind the reader the basics
of the mass and angular momentum transport in accretion disks in §3. Our
simulation results on the two kinds of transport are described in §4 and 5,
respectively. We examine the different contributions to the angular momentum
budget in §6 and describe the correlation between the angular momentum and
mass fluxes found in our simulations in §7. The wave-driven evolution of the
inner disk is characterized in §8. Finally, we discuss and summarize our
results in §9 and 10, respectively.

##  2 Description of the numerical data set

Our present study is based on the data produced by a set of 2D simulations (39
runs in total) presented in Paper I. These runs used Athena++ (Stone et al.,
2020) to solve hydrodynamic equations with the globally isothermal equation of
state (EOS), i.e. $P=\Sigma c_{s}^{2}$ for a constant sound speed $c_{s}$,
where $P$ is the vertically integrated pressure and $\Sigma$ is the surface
density. Magnetic fields and disk self-gravity have been ignored, and the
accretor was initially non-spinning.

These simulations are set in polar $r-\phi$ coordinates, with the azimuthal
coordinate covering full $2\pi$. In the radial direction the simulation domain
starts below the stellar radius $R_{\star}$ and extends in the disk out to
$4R_{\star}$. The radial grid is uniformly spaced in $\log r$ to increase
resolution in the vicinity of the BL. Simulations were run for every integer
value of ${\mathcal{M}}$ in the interval $5\leq{\mathcal{M}}\leq 15$. For
three values of ${\mathcal{M}}=6,9,12$ we run multiple simulations with the
different forms of the initial noise spectrum (triggering sonic instability in
the BL) and numerical resolution, to test their impact on the outcomes.

A distinctive feature of our simulations is their detailed, purpose-built,
real-time analysis. They employed a high-cadence sampling of the outputs,
allowing a highly informative analysis of the outputs to be performed. In
particular, we ran on-the-fly fast-Fourier transforms (FFTs) of various fluid
variables, giving us new, previously unobtainable diagnostic capabilities.
This allowed us to detect and characterize a number of different modes present
in the system, with their distinct azimuthal wavenumbers $m$ and pattern
speeds $\Omega_{P}$, at every moment of time. The ability to automatically
detect and analyze different modes present in the system is an important
improvement of the study presented in Paper I compared to the existing works.

In the following, when presenting our results, we will be using units in which
$R_{\star}=1$ and $\Sigma(R_{\star})=1$ at $t=0$; in these units the initial
profile of the surface density is $\Sigma=r^{-3/2}$. We also set the Keplerian
velocity at the surface of the star $v_{K}(R_{\star})$ to be unity, which
implies that $c_{s}=\mathcal{M}^{-1}$ and $GM_{\star}=1$. This choice makes
the Keplerian period at $R_{\star}$ to be $\tau_{\star}=2\pi$, and we will
often express time in the form of $t/2\pi$, or in the units of $\tau_{\star}$.

###  2.1 Main findings of Paper I

We now recount the main conclusions of the wave morphology study presented in
Paper I, that will allow us to better interpret the results of the current
work.

Our suite of simulations reveals a complicated pattern of wave activity in the
vicinity of the BL, with a number of different modes operating in the inner
disk. In addition to the upper and lower acoustic modes previously described
in Belyaev et al. (2013a), we have also discovered a new type of modes, so
called vortex-driven modes. They appear as global spiral arms extending from
the vicinity of the BL into the upper disk, and owe their existence to the
localized vortex-like structures that form in the disk next to the BL and
launch these density waves.

Another type of modes that we routinely detect in our runs are the so called
resonant modes (owing their existence to a particular geometric resonance
condition, Belyaev et al. 2012), which are present only in the disk. Similar
to the lower acoustic modes, the resonant modes are trapped between the
stellar surface and the inner Lindblad resonance, whose location
$r_{\mathrm{ILR}}$ in a Keplerian disk with
$\Omega=\Omega_{K}=(GM_{\star}/r^{3})^{1/2}$ is given by

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

In Paper I we also found a general tendency of the azimuthal wavenumber $m$ of
the most prominent modes to increase with ${\mathcal{M}}$. These results
appear robust with respect to variations of the initial conditions and
numerical resolution.

##  3 Wave-driven mass and angular momentum transport: basics

The main goal of the present work is to better understand how matter and
angular momentum are transported through the inner disk and the BL, resulting
in accretion onto the star. To that effect, we examine how the characteristics
of the different modes that we see in our runs are linked to the global
evolution of the star-disk system. We do this by exploring the behavior of
certain variables derived from our simulations — mass and angular momentum
fluxes — for different values of ${\mathcal{M}}$ and connecting them to the
changes of the disk properties.

In all our simulations transport of mass is tracked by measuring the mass
accretion rate,

$\displaystyle\dot{M}\left(r\right)\equiv-r\int_{0}^{2\pi}v_{r}\left(r,\phi\right)\Sigma\left(r,\phi\right)\,{\rm
d}\phi=-2\pi r\langle\Sigma v_{r}\rangle,$ (2)

(see Figure 1) where we introduced a shorthand notation

$\displaystyle\langle f(r)\rangle=(2\pi)^{-1}\int_{0}^{2\pi}f(r,\phi)d\phi$
(3)

for azimuthal averaging of any variable $f(r,\phi)$. Accretion rate is defined
such that $\dot{M}>0$ for inflow of mass towards the star.

To characterize angular momentum transport we adopt the standard procedure of
measuring the total angular momentum flux (AMF) $C_{\rm L}$, defined as

$\displaystyle C_{\rm L}=2\pi r^{2}\left<\Sigma v_{r}v_{\phi}\right>,$ (4)

where $v_{r}$ and $v_{\phi}$ are the radial and azimuthal velocity components.
By introducing a reference azimuthal velocity of the fluid $v_{\phi,0}(r)$, to
be specified later, and the azimuthal velocity perturbation $\delta
v_{\phi}=v_{\phi}-v_{\phi,0}$, we can further decompose $C_{\rm L}$ into the
advective angular momentum flux $C_{\rm A}$, and the wave angular momentum
flux222Despite the typographic similarities $c_{\rm s}$ and $C_{\rm S}$ are
used to denote two distinct quantities: sound speed and wave-induced angular
momentum flux, respectively. (or Reynolds stress) $C_{\rm S}$, defined as
follows:

$\displaystyle C_{\rm A}$ $\displaystyle=2\pi r^{2}\left<\Sigma
v_{r}\right>v_{\phi,0}=-\dot{M}~{}rv_{\phi,0},$ (5) $\displaystyle C_{\rm S}$
$\displaystyle=2\pi r^{2}\left<\Sigma
v_{r}\left(v_{\phi}-v_{\phi,0}\right)\right>=2\pi r^{2}\left<\Sigma
v_{r}\delta v_{\phi}\right>,$ (6)

so that $C_{\rm L}=C_{\rm A}+C_{\rm S}$. In our case the stress contribution
$C_{\rm S}$ arises primarily because of the oscillatory wave-like motions and
is driven by the wave modes rather than some turbulence, as would be the case
for e.g. the MRI.

In this work, when analyzing simulation outputs, we will most often be
choosing $v_{\phi,0}$ to be the mean azimuthal velocity of the fluid,

$\displaystyle v_{\phi,0}=\left<v_{\phi}\right>,$ (7)

following Fromang & Nelson (2006) and Flock et al. (2011), among others. Note
that alternative definitions of $C_{\rm S}$ using other choices of
$v_{\phi,0}(r)$ are also possible, see e.g. Ju et al. (2016), Arzamasskiy &
Rafikov (2018), as well as §3.1.

Because of oscillatory, intrinsically time-dependent nature of the acoustic
mode-driven transport, we find $\dot{M}$ and the different angular momentum
flux contributions to be highly time-dependent. For that reason, we resort to
averaging these physical quantities over sufficiently long time intervals to
provide meaningful comparison between them. We describe mathematical details
of such averaging procedures in Appendix A.

###  3.1 Relation between $\dot{M}$ and $C_{\rm S}$

The mass and angular momentum fluxes through the disk are closely related to
each other, which can be demonstrated quite generally starting from the
hydrodynamic equations of motion and continuity. In particular, it was shown
in Belyaev et al. (2013a) that by choosing $v_{\phi,0}(r)$ in the mass-
weighted form (Balbus & Hawley, 1998; Balbus & Papaloizou, 1999)

$\displaystyle v_{\phi,0}=v_{\Sigma}\equiv\frac{\left<\Sigma
v_{\phi}\right>}{\left<\Sigma\right>},$ (8)

different from that given by the equation (7), and defining the wave AMF (6),
or Reynolds stress, $C_{\rm S}$ accordingly (i.e. with
$v_{\phi,0}=v_{\Sigma}$), the equation representing conservation of the
angular momentum can be cast in the following form:

$\displaystyle\dot{M}\frac{\partial l}{\partial r}=\frac{\partial C_{\rm
S}}{\partial r}+2\pi
r^{3}\langle\Sigma\rangle\frac{\partial\Omega_{0}}{\partial t}.$ (9)

Here $\Omega_{0}=v_{\Sigma}/r$ and $l=rv_{\Sigma}=\Omega_{0}r^{2}$ are the
angular frequency and the specific angular momentum corresponding to the
reference azimuthal velocity $v_{\Sigma}$. If instead of (8) we adopted
$v_{\phi,0}$ in the form (7), then the last term in the equation (9) would
look differently.

Equation (9) directly relates the transport of mass — the term proportional to
$\dot{M}$ in the left hand side — to the transport of the angular momentum in
the disk — the divergence of the wave angular momentum flux $C_{\rm S}$ in the
right-hand side. The second term in the right hand side, usually neglected in
studies of accretion processes, arises in disks which evolve sufficiently
rapidly for their mean angular frequency $\Omega_{0}$ to change in time. If we
neglect this contribution to the angular momentum balance for a moment, then
we find

$\displaystyle\dot{M}=\left(\frac{\partial l}{\partial
r}\right)^{-1}\frac{\partial C_{\rm S}}{\partial r}=\frac{\partial C_{\rm
S}}{\partial l}.$ (10)

Note that this expression is fully general and applies for any choice of
$v_{\phi,0}$; a particular form of $v_{\phi,0}$ affects the definition of both
$C_{\rm S}$ and $l$.

In truly viscous disks with kinematic viscosity $\nu$ one should use viscous
angular momentum flux $F_{J}=-2\pi\nu\Sigma r^{3}d\Omega/dr$ instead of
$C_{\rm S}$ (Lynden-Bell & Pringle, 1974). In particular, in disks with
radially constant $\dot{M}$ one finds $F_{J}=\dot{M}l$.

However, as we will demonstrate in §6, the second term in the right hand side
of equation (9) often plays an important role near the BL. This contribution
to the angular momentum balance has been previously considered in Belyaev et
al. (2013a), and its relation to the evolution of the disk properties was
explored in Arzamasskiy & Rafikov (2018).

##  4 Transport of mass

We start by describing our results for the mass flux $\dot{M}$. Given that
$\dot{M}$ is closely related to $C_{\rm A}$, see equation (5), and in the disk
$v_{\phi,0}$ is typically pretty close to the Keplerian velocity (see §8.1).
One can then approximate $C_{\rm A}\approx-\dot{M}l_{K}$ outside the BL, where
$l_{K}(r)=\sqrt{GM_{\star}r}$ is the specific angular momentum for a Keplerian
disk. The profile of $C_{\rm A}(r)$ is plotted in Figures 2-3 (blue curve in
the middle row of each individual subpanel) for simulations at varied
${\mathcal{M}}$ and at different moments of time. It is discussed in more
details in §5.1, but for now we will mention some key features pertinent for
the $\dot{M}$ behavior.

First, in all plotted cases $C_{\rm A}(r)$ (and $\dot{M}$) exhibits a deep
minimum (maximum) close to the BL. This is naturally explained by the enhanced
wave activity near the BL, and results in important consequences for the
surface density evolution in the inner disk discussed in §8. Second, in most
cases both $C_{\rm A}(r)$ and $\dot{M}$ substantially diminish in amplitude
far from the accretor. We will relate this behavior to the angular momentum
fluxes carried by the different types of waves in §5.

###  4.1 $\dot{M}$ decomposition

Our simulations reveal an interesting feature of the $\dot{M}$ behavior near
the BL, that we discuss next. Let us write down
$\Sigma(r,\phi)=\langle\Sigma\rangle+\delta\Sigma(r,\phi)$, where
$\delta\Sigma(r,\phi)$ is the perturbation of the surface density relative to
its azimuthally-averaged value $\langle\Sigma\rangle$. Then the definition (2)
can be cast as

$\displaystyle\dot{M}\left(r\right)=-2\pi r\langle\Sigma\rangle\langle
v_{r}\rangle-2\pi r\langle\delta\Sigma v_{r}\rangle,$ (11)

illustrating the separation of $\dot{M}$ into two distinct contributions — the
one due to the advection of the mean surface density $\langle\Sigma\rangle$
(first term) and the contribution due to the correlation between the non-
axisymmetric fluctuation of $\Sigma$ and $v_{r}$. We illustrate the behavior
of these contributions in Figure 1 for three simulations with different values
of ${\mathcal{M}}$. One can see that in all cases both contributions to
$\dot{M}$ are, in general, equally important: they are comparable in magnitude
and tend to offset each other so that the total $\dot{M}$ is considerably
lower in magnitude than each of these terms over large radial ranges (which is
especially noticeable near the BL). For that reason, one cannot directly
connect $\dot{M}(r)$ to the radial profiles of $\langle v_{r}\rangle$ and
$\langle\Sigma\rangle$.

Importance of the $\dot{M}$ contribution due to the correlation between
$\delta\Sigma$ and $v_{r}$ is one of the distinctive features of the global,
acoustic wave-driven transport. In the classical picture of laminar viscous
disk accretion (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974) this
contribution is identically zero since the disk is axisymmetric and
$\delta\Sigma=0$. In our case, this contribution is non-zero in the vicinity
of the BL because of highly correlated $\delta\Sigma$ and $v_{r}$ for wave-
like fluid motions and the nonlinear dissipation of these waves, which leads
to damping of the angular momentum flux carried by the waves. Interestingly,
$\langle\delta\Sigma v_{r}\rangle$ is often more regular than
$\langle\Sigma\rangle\langle v_{r}\rangle$, since the latter sometimes
exhibits small-scale spatial variability even upon the long-term time
averaging, see Figure 1b.

We are not aware of any studies mentioning the role of the mass transport term
proportional to $\langle\delta\Sigma v_{r}\rangle$ in simulations featuring
turbulence, self-consistently driven by a local hydrodynamic (e.g. the
vertical shear instability) or MHD (e.g. MRI) mechanism. Examining this issue
may be interesting for figuring out whether in our case the non-zero second
term in equation (11) arises due to the intrinsically global nature of the
wave-driven transport near the BL (see §9.2), or it is also present when
accretion is mediated by local processes.

Figure 1: Different contributions to the total mass accretion rate $\dot{M}$
(black curve) in equation (11), plotted as a function of $r$ for several
simulations (see Paper I for their naming notation). The advective term
$2\pi\langle\Sigma v_{r}\rangle\langle v_{r}\rangle$ is shown in blue, while
the diffusive term $-2\pi\langle\delta\Sigma v_{r}\rangle$ is the orange
curve. Averaging over the interval $t/2\pi=100-200$ is performed.

##  5 Angular momentum transport

We now turn to the details of the angular momentum transport in the vicinity
of the BL. Figures 2,3 illustrate the behavior of various representative
variables during three different periods of time for the four representative
simulations with ${\mathcal{M}}=7,9,12$ and $15$. The morphological
description of the modes emerging in two of the runs — for ${\mathcal{M}}=9$
(M09.FR.r.a) and ${\mathcal{M}}=15$ (M15.FR.r.a) — has been provided in Paper
I, whereas the ${\mathcal{M}}=7$ (M07.FR.r.a) and ${\mathcal{M}}=12$
(M12.FR.mix.a) runs have not been discussed in detail in Paper I.

For each ${\mathcal{M}}$ we select three time intervals of length $50$ or
$100$ orbits and perform averages (as described in Appendix A) of various
angular momentum fluxes — total $C_{\rm L}$, advective $C_{\rm A}$, and wave
$C_{\rm S}$, defined by the equations (4)-(6) with
$v_{\phi,0}=\left<v_{\phi}\right>$ given by the equation (7). Their radial
profiles are shown in the middle row of each panel. As mentioned earlier in
§4, the behavior of the advective flux $C_{\rm A}$ in the disk is closely
related to the radial profile of the mass accretion rate $\dot{M}$.

We also select a representative moment of time within each interval, which
illustrates a set of modes typical for that period, and show snapshots of
various fluid variables in the top row the panel corresponding to that time
interval. The format is the same as in Figure 2 of Paper I, namely, we show a
polar map of the wave amplitude variable $rv_{r}\sqrt{\Sigma}$ (right), a
Cartesian map of the same variable over the reduced radial range to highlight
the near-BL details (middle), and a Cartesian map of the vperturbation of
ortensity $\zeta=\omega/\Sigma$ (here $\omega=\nabla\times\mathbf{v}$ is the
vorticity) relative to its initial value near the BL (left).

Figure 2: Figure illustrating angular momentum flux data for simulations with
${\mathcal{M}}=7$ (run M07.FR.r.a, left panels) and ${\mathcal{M}}=9$ (run
M09.FR.r.a, right panels). Each panel (labeled with capital letters) consists
of 5 sub-panels, labeled by their subscript and described in the beginning of
§5. Every panel is associated with a time interval over which the AMF data are
computed ($t/2\pi$, shown above the panel) and a representative moment of time
$T/2\pi$ within this interval, indicated as follows.

Finally, the bottom row shows the decomposition of the wave angular momentum
flux $C_{\rm S}$ into the contributions provided by the individual Fourier
modes of the fluid perturbation, which are dominant through the plotted time
interval (in a time integrated sense). To second order in perturbed variables
the Fourier contribution $C_{{\rm S},m}$ to $C_{\rm S}$ from the $m$-th mode
is (Belyaev et al., 2013a)

$\displaystyle C_{{\rm S},m}=2\pi
r^{2}\left<\Sigma\right>\left(v_{r,m}v_{\phi,m}^{\ast}+v_{r,m}^{\ast}v_{\phi,m}\right),$
(12)

where the asterisks denote complex conjugates, so that $C_{\rm
S}\approx\sum_{m=1}^{\infty}C_{{\rm S},m}$. The particular $C_{{\rm S},m}$
profiles shown in the figure are for the dominant modes chosen according to
the procedure described in Appendix A. The gray curves, representing the sum
of only the dominant harmonics indicated in each panel, sometimes show
deviations from the full $C_{\rm S}$ (black), computed using equation (6).
These differences arise because a number of other, less significant modes not
shown in these plots also contribute to the full $C_{\rm S}$, and also because
equation (12) is only second order accurate in fluid perturbations.

The decomposition of $C_{\rm S}$ into contributions from the different
azimuthal harmonics allows us to explore the role played by the individual
modes (which, once their $m$ is known, can be easily associated with the
different types of waves) in transporting the angular momentum in the vicinity
of the BL. This is the key advantage of our present work compared to Belyaev
et al. (2013a) and Dittmann (2021), who did not perform such decomposition.
Another improvement is the larger range of Mach numbers explored in our study,
providing us with a better understanding of the dependence of transport
properties on ${\mathcal{M}}$.

###  5.1 Different AMF contributions

Middle sub-panels of Figures 2-3 display the advective angular momentum flux
$C_{\rm A}$, the wave flux $C_{\rm S}$, and their sum $C_{\rm L}$ on the same
scale. A notable feature of the $C_{\rm A}$ curves is the small-scale
oscillations that they often exhibit. These oscillations are likely caused by
the intrinsic time-variability, since the inner disk is pervaded by the
multiple wave modes. This variability ends up manifesting itself in the
spatial domain, even though we perform time averaging over a substantial
interval of time, see Appendix A.

One can see that $|C_{\rm A}|$ peaks not too far from the stellar surface and
its maximum amplitude is always larger than that of $C_{\rm S}$ by a factor of
several. This makes sense since, to zeroth order (certainly in steady state),
$C_{\rm A}\propto\dot{M}\propto\partial C_{\rm S}/\partial r$, see equations
(5)), (9). Since $C_{\rm S}$ typically varies over a narrow region near the
BL, its radial derivative can reach high values. Naturally, the narrower is
the near-BL region over which $C_{\rm S}$ varies, the larger should we expect
the difference in amplitudes between $C_{\rm A}$ and $C_{\rm S}$ to be. And
indeed, in Figure 2 (for ${\mathcal{M}}=7$) we find $C_{\rm S}$ variation in
the range $1<r\lesssim 1.5$ to result in the ratio of the maximum amplitudes
$C_{\rm A}$ and $C_{\rm S}$ being $\sim(2-2.5)$, while in Figure 3Fd (for
${\mathcal{M}}=15$) $C_{\rm S}$ variation in the range $1<r\lesssim 1.2$ leads
to the maximum ratio of fluxes $\sim 5-6$.

In most cases we find $C_{\rm A}$ to be negative, as expected for mass inflow
with $\dot{M}>0$. Notable exceptions include ${\mathcal{M}}=9$ run at
$t/2\pi=400-500$ (Figure 2Fd) and ${\mathcal{M}}=15$ run at $t/2\pi=100-200$
(Figure 3Dd), for which $C_{\rm A}>0$ (i.e. $\dot{M}<0$) at large distances.
This can again be understood on the basis of the equations (9) or (10), by
noticing that in these cases $C_{\rm S}$ is positive and decays with $r$ far
from the star as a result of dissipation (it will become clear in §5.2,5.3
that this behavior is caused by the nonlinear damping of the vortex-driven
modes). As a consequence, $\partial C_{\rm S}/\partial r<0$ there, leading to
$\dot{M}<0$ (outflow) and positive $C_{\rm A}$.

Figure 3: Same as Fig. 2 but for ${\mathcal{M}}=12$ (run M12.FR.mix.lc.a, left
panels) and ${\mathcal{M}}=15$ (run M15.FR.r.a, right panels). For
${\mathcal{M}}=12$: (A) $T/2\pi=150$, $t/2\pi=100-200$; (B) $T/2\pi=300$,
$t/2\pi=250-350$; (C) $T/2\pi=575$, $t/2\pi=550-600$. For ${\mathcal{M}}=15$:
(D) $T/2\pi=150$, $t/2\pi=100-200$; (E) $T/2\pi=350$, $t/2\pi=300-400$; (F)
$T/2\pi=550$, $t/2\pi=500-600$.

###  5.2 Wave AMF ${C_{\rm S}}$ as a function of ${\mathcal{M}}$ and $m$

We now describe the behavior of the wave AMF $C_{\rm S}$ for different values
of ${\mathcal{M}}$ shown in Figures 2,3 in conjunction with the morphological
characteristics of the different waves present in the system.

####  5.2.1 ${\mathcal{M}}=7$ run

We start our description of the ${C_{\rm S}}$ behavior with the
${\mathcal{M}}=7$ run, shown in Figures 2A-C. At early times333Note that by
that time a significant depression in surface density has already formed in
the inner disc, see Figure 4a. ($t/2\pi=125$) the inner disk is dominated by
the $m=11$ lower mode trapped within $r\lesssim 1.5$ (panel Ab), which roughly
corresponds to the $r_{\rm ILR}$ of this mode (its $\Omega_{P}\approx 0.46$).
This mode gives rise to the negative AMF $C_{\rm S,11}$, see the blue
curve444The orange curve in that panel corresponding to the $m=10$ lower mode
also provides a substantial contribution to ${C_{\rm S}}$ near the star. This
mode is not seen in panels (Aa)-(Ac) because the data shown in panels (Ad),
(Ae) are integrated over an extended time interval and $m=10$ signal gets
washed out. in panel (Ae). As a result, ${C_{\rm S}}$ starts negative at
$r=R_{\star}$ and crosses zero around $r_{\rm ILR}$.

Outside the resonant cavity there is a set of rather incoherent vortex-driven
spiral arms, excited by a number of chaotic vortices residing at $r\lesssim
1.5$, visible in panel (Aa). These weak arms drive low-amplitude but positive
$C_{\rm S}$ outside $r=1.5$, see panel (Ad).

By $t/2\pi=225$ the amplitude of the wave activity decreases; $m=11$ lower
mode gets replaced with the lower $m=10$ mode, which is prominent inside the
star but is rather weak in the disk, see panel (Bb). As a result, the
magnitude of the (negative) segment of the $C_{\rm S}$ curve in the inner disk
drops by about two orders of magnitude, and it gets confined close to the star
(panel (Bd)). Also by that time, vortensity evolution produces three prominent
vortices at $r\approx 1.1$ (panel (Ba)), which excite strong vortex-driven
spirals in the outer disk, increasing the amplitude of $C_{\rm S}>0$ in the
outer disk (panel (Bd)).

In the final stages of this ${\mathcal{M}}=7$ simulation, at $t/2\pi=575$, a
lower $m=16$ mode clearly shows in the disk at $r\lesssim 1.8$ (although
inside the star $m=8$ mode dominates), although with rather modest amplitude
(panel (Cb)). This, again, makes ${C_{\rm S}}$ negative close to the star, see
the red curve for $C_{{\rm S},16}$ in panel (Ce). Also, the vortensity map in
panel (Ca) reveals two strong, azimuthally elongated vortices centered at
$r\approx 1.5$. These "rolls", as we termed them in Paper I, drive a pair of
azimuthally extended vortex-driven spiral arms with $m=2$ in the outer disk.
Their $C_{{\rm S},2}>0$ makes ${C_{\rm S}}$ positive in the outer disk (see
panel (Ce)), with the amplitude $\approx 2.5$ times higher than at
$t/2\pi=225$.

Already this run alone demonstrates that depending on a type of the mode
dominating in a particular region of the disk, one can get rather different
behaviors of ${C_{\rm S}}$. Our subsequent description of the runs for other
values of ${\mathcal{M}}$ will reinforce this conclusion.

####  5.2.2 ${\mathcal{M}}=9$ run

The ${\mathcal{M}}=9$ run in its early phases $t/2\pi=150$ reveals a set of
lower $m=19-21$ modes confined within $r_{\rm ILR}\approx 2.2$, corresponding
to their $\Omega_{P}\approx 0.3$, see Figure 2Db,Dc. These modes give rise to
${C_{\rm S}}<0$ near the star. Similar picture persists also at $t/2\pi=275$,
see panels (Eb)-(Ee). At the same time, in the outer disk ${C_{\rm S}}$ stays
close to zero: unlike the late stages of the ${\mathcal{M}}=7$ run, coherent
vortices do not emerge near the star (see panels (Da), (Ea)), so vortex-driven
spirals are very weak.

An interesting feature of the ${\mathcal{M}}=9$ run pointed out in Paper I, is
the emergence of $m=2$, radially elongated mode with very low pattern speed
$\Omega_{P}\approx 0.15$. It is readily noticeable in panel (Eb), however its
signal is not present in panel (Ee), which means that its AMF contribution
$C_{{\rm S},2}$ is negligible. This is hardly surprising since this mode has
radial wave number $k_{r}=0$.

Finally, at $t/2\pi=150$ the picture changes dramatically, see panels
(Fb),(Fc). Near the star, for $r\lesssim 1.4$, the disk is dominated by the
$m=6$ resonant mode (see Paper I for details), whereas in the outer disk there
are global $m=5$ spiral arms driven by a set of five vortex rolls located at
$r\approx 1.3$, see panel (Fa). The resonant mode gives rise to a small (but
negative near the star) AMF contribution $C_{{\rm S},6}$, whereas the vortex-
driven modes produce $C_{{\rm S},5}>0$ in the outer disk, similar to
${\mathcal{M}}=7$ run.

####  5.2.3 ${\mathcal{M}}=12$ run

An ${\mathcal{M}}=12$ run that we use in this work is different from the run
employed in Paper I. This is done in part to illustrate the fact, pointed out
in Paper I, that all ${\mathcal{M}}=12$ runs look very similar to each other,
in a much more homogeneous way than for any other value of ${\mathcal{M}}$.

Indeed, comparing Figure 3A-C with the Fig. 7 of Paper I, one can see a very
familiar pattern: for most of the simulation, including $t/2\pi=150$ and $300$
shown in Figure 3A,B, the inner disk (out to $r\approx 1.6$) is dominated by
the lower $m=16$ mode, accompanied by a set of narrow vortex-driven spirals
clearly visible in the outer disk. As expected, this pattern results in
$C_{{\rm S},16}$ driving ${C_{\rm S}}<0$ near the star, with ${C_{\rm S}}$
becoming positive further out, at $r\gtrsim 1.3-1.5$, see panels (Ae), (Be).

Closer to the end of the run, at $t/2\pi=550$, one can see a combination of a
resonant and a lower $m=11$ modes driving ${C_{\rm S}}$ below zero near the
star, with the weak vortex-driven modes carrying ${C_{\rm S}}>0$ outside
$r\approx 1.3$.

####  5.2.4 ${\mathcal{M}}=15$ run

The ${\mathcal{M}}=15$ simulation is the highest ${\mathcal{M}}$ run carried
out as a part of our simulation suite described in Paper I. Unlike other
simulations described so far, this run shows quite substantial upper acoustic
mode activity inside the star (Belyaev et al., 2013a), see Figure 3Db,Eb,Fb,
although the manifestations of this mode type in the disk are not as clear.

In this run we witness the appearance of 7-8 vortex rolls at $r\approx 1.15$
rather early on, already by $t/2\pi=150$, see panel Da. These rolls drive a
set of global, relatively narrow spiral arms propagating in the outer disk and
making ${C_{\rm S}}$ positive beyond $r\approx 1.1$. For $1<r\lesssim 1.1$ we
find ${C_{\rm S}}<0$, which could be due to the resonant modes starting to
develop in this part of the disk.

These resonant modes do become rather prominent (and more radially extended)
at later time, e.g. at $t/2\pi=350$ and $550$, see panels (Eb),(Fb). They keep
${C_{\rm S}}$ negative within $r\approx 1.2$, see panels (Ee),(Fe). At the
same time, the vortex rolls appearing earlier in the simulation merge into two
(at $t/2\pi=350$) and then one vortex (at $t/2\pi=550$), see panels (Ea),
(Fa). They give rise to the vortex-driven mode activity in the outer disk,
which maintains positive ${C_{\rm S}}$ there, albeit at a reduced amplitude
compared to $t/2\pi=150$.

###  5.3 Transport by the different types of modes

Examination of the bottom sub-panels in all panels of Figures 2-3 reveals that
at any moment of time there is often only one or two modes in the disk that
dominate the wave AMF $C_{\rm S}$. For example, in the ${\mathcal{M}}=7$ run
at $t/2\pi=100-200$ we find two lower modes with $m=10,11$ to provide the
dominant contribution to $C_{\rm S}$ for $r<1.4$ (Figure 2Ae), while at
$t/2\pi=550-600$ the vortex-driven $m=2$ mode dominates $C_{\rm S}$ for
$r>1.4$ (Figure 2Ce); in the ${\mathcal{M}}=12$ run the lower $m=16$ mode
dominates $C_{\rm S}$ for $r<1.5$ at both $t/2\pi=100-200$ and $250-350$
(Figure 3Ae,Be). As these dominant modes have a particular sign of $C_{\rm
S}$, they affect disk evolution in a certain way, as described later in §8.

A wave propagating in the disk carries angular momentum $\Delta J_{\rm
w}\propto[\Omega_{P}-\Omega(r)]$ (see below), so that the flux of angular
momentum associated with it is

$\displaystyle C_{\rm S}\propto k_{r}\Delta J_{\rm w}\propto
k_{r}[\Omega_{P}-\Omega(r)],$ (13)

where $k_{r}$ is the radial wavenumber. In particular, an outgoing ($k_{r}>0$)
wave has $C_{\rm S}>0$ in the region where $\Omega_{P}>\Omega(r)$ (e.g. a
vortex-driven mode far from the star), whereas an incoming ($k_{r}<0$) wave
has $C_{\rm S}>0$ in the region where $\Omega_{P}<\Omega(r)$ (e.g. a lower
mode reflected off the ILR). This is what one also finds for the density waves
excited by planets, outside and inside of the planetary semi-major axis. On
the other hand, in the disk region where $\Omega_{P}<\Omega(r)$ an outgoing
($k_{r}>0$) wave has $C_{\rm S}<0$ (e.g. a lower mode propagating away from
the star).

Nonlinear wave damping provides an intuitive way to understand why $\Delta
J_{\rm w}\propto{\rm sgn}[\Omega_{P}-\Omega(r)]$. Let us consider a shock wave
with a pattern speed $\Omega_{P}$. Whenever $\Omega_{P}>\Omega(r)$ (which is
true far from the accretor) the shock overtakes the disk fluid, accelerating
it in the azimuthal direction upon crossing the shock as a result of the shock
jump conditions. As a result, disk fluid receives positive $\Delta J$ from the
wave. On the contrary, for $\Omega_{P}<\Omega(r)$ (fulfilled close to the BL,
in the inner disk) the disk fluid catches up with the shock front and
decelerates upon crossing it; this means that the wave has transferred
negative angular momentum to the disk, i.e. that $\Delta J_{\rm w}<0$.

####  5.3.1 Transport by the lower acoustic modes

Lower modes are trapped inside the resonant cavity at $1<r<r_{\rm ILR}$, and
in this part of the disk $\Omega_{P}<\Omega(r)$. The mode is excited in the BL
(Belyaev et al., 2013a) and first propagates out (i.e. $k_{r}>0$) towards the
ILR. At the ILR (where its $k_{r}\to 0$) it gets reflected and propagates back
to the stellar surface with $k_{r}<0$ (see Figure 3Ab for an illustration). As
a result, the outgoing lower mode has angular momentum flux $C_{\rm S}^{\rm
out}<0$, whereas the incoming (reflected) one carries $C_{\rm S}^{\rm in}>0$.
The full angular momentum flux is the sum of the two,

$\displaystyle C_{\rm S}(r)=C_{\rm S}^{\rm out}(r)+C_{\rm S}^{\rm in}(r).$
(14)

Because of the nonlinear wave damping, $|C_{\rm S}^{\rm out}(r)|$ gradually
decreases in amplitude as $r$ increases, whereas $|C_{\rm S}^{\rm in}(r)|$
decreases in amplitude with decreasing $r$. As a result, for any $1<r<r_{\rm
ILR}$ one finds that the negative $C_{\rm S}^{\rm out}(r)<C_{\rm S}^{\rm
out}(r_{\rm ILR})$, whereas the positive $C_{\rm S}^{\rm in}(r)<C_{\rm S}^{\rm
in}(r_{\rm ILR})$. Also, $C_{\rm S}^{\rm out}(r_{\rm ILR})=-C_{\rm S}^{\rm
in}(r_{\rm ILR})$ since $C_{\rm S}(r_{\rm ILR})=0$. With this in mind,
equation (14) can be cast as

$\displaystyle C_{\rm S}(r)=\left[C_{\rm S}^{\rm out}(r)-C_{\rm S}^{\rm
out}(r_{\rm ILR}\right]+\left[C_{\rm S}^{\rm in}(r)-C_{\rm S}^{\rm in}(r_{\rm
ILR})\right]<0,$ (15)

showing that ${C_{\rm S}}$ of the lower modes is always negative, in agreement
with Figures 2-3. As a result, these modes drive mass flow towards the
accretor, $\dot{M}>0$, see §8.

####  5.3.2 Transport by the resonant modes

Resonant modes are very similar to the lower modes, since they are also
trapped in the resonant cavity, where $\Omega_{P}<\Omega(r)$. They also
propagate both towards and away from the ILR. As a result, by applying the
same logic as in §5.3.1, we conclude that resonant modes have $C_{\rm S}(r)<0$
and their dissipation gives rise to $\dot{M}>0$.

Typically, the magnitude of $C_{\rm S}$ for the resonant modes is lower than
$|C_{\rm S}|$ of the lower modes. This is caused by the smaller azimuthal
wavenumber $m$ of the resonant modes: lower $m$ implies azimuthally wider wave
profile, slower nonlinear wave steepening and shocking, and weaker
dissipation. As a result, for resonant modes $C_{\rm S}^{\rm out}(r)$ and
$C_{\rm S}^{\rm in}(r)$ in equation (14) are closer in amplitude to each other
(but still different in sign) than for the lower modes, resulting in their
more substantial cancellation and lower $|C_{\rm S}|$.

####  5.3.3 Transport by the vortex-driven modes

Vortex-driven modes have $\Omega_{P}$ set by the angular frequency of their
parent vortices, which is typically high since the vortices reside not too far
from the BL. As a result, these modes propagate with $k_{r}>0$ and
$\Omega_{P}>\Omega(r)$ in the outer disk, outside the corotation region of the
mode. This means that these modes have $C_{\rm S}>0$ and (in steady state)
result in $\dot{M}<0$, i.e. mass outflow in the outer disk.

As shown in Paper I, vortex-driven waves typically fall into two classes:
those produced by the compact, isolated vortices (e.g. see Figure
2Ba,Bb,3Aa,Ab), and the ones driven by the more azimuthally elongated "rolls"
(e.g. see Figure 2Ca,Cb,Fa,Fb). The former have small azimuthal width and are
superpositions of a number of high-$m$ perturbation modes (see Paper I). As a
result, their constituent modes are often not singled out individually in our
plots, see e.g. the $C_{\rm S}>0$ segments in the outer disk in Figures
2Be,3Ae, which are not dominated by a mode with a single value of $m$. On the
contrary, the azimuthally extended rolls typically produce waves with a well
defined, low value of $m$. As a result, in such cases a particular $C_{{\rm
S},m}$ dominates the full $C_{\rm S}$ in the outer disk, see Figures 2Ce,Fe.

The way in which vortex-driven modes dissipate depends on their amplitude, as
well as ${\mathcal{M}}$ and $m$. Higher $m$ accelerates wave evolution into a
shock because of the smaller azimuthal scale. For that reason, the vortex-
driven $m=2$ mode in the ${\mathcal{M}}=7$ run shown Figure 2Ce experiences
little damping (its $C_{{\rm S},2}$ is essentially constant in the outer
disk), whereas the $m=5$ vortex-driven mode in the ${\mathcal{M}}=9$ run shown
Figure 2Fe starts damping appreciably outside $r\approx 2.5$. Note that these
modes have similar ${\mathcal{M}}$ and amplitude, so it is the difference in
their $m$ that is responsible for the difference in their decay. Also,
everything else being equal, higher amplitude and higher ${\mathcal{M}}$
(Goodman & Rafikov, 2001; Rafikov, 2002) of a wave promote its faster shocking
and dissipation.

####  5.3.4 Transport by the upper acoustic modes

Upper acoustic modes show up in the disk early on in low ${\mathcal{M}}$ runs.
At high ${\mathcal{M}}$ they also operate inside the star over long periods of
time, see Paper I. These modes are evanescent near the BL, where their
$\Omega_{P}<\Omega(r)$, but become propagating outside their Outer Lindblad
Resonance, where their $\Omega_{P}>\Omega(r)$. In that regard, they are
similar to the vortex-driven modes considered above, meaning that they promote
mass outflow in the disk, $\dot{M}<0$. However, the amplitude of the $\dot{M}$
driven by these modes in the disk is typically too low to make them a
significant agent in the disk evolution.

####  5.3.5 Transport by other modes

Another interesting mode type that we see in our simulations is the low $m=2$
mode in ${\mathcal{M}}=9$ runs clearly visible in Figure 2Eb. This mode has a
characteristic radially elongated perturbation pattern (with a clear shift of
phase around $r\approx 1.5$) with $k_{r}=0$. This implies that this mode
propagates purely azimuthally and thus does not transport angular momentum in
the radial direction. This conclusion is supported by the lack of $C_{{\rm
S},2}$ among the different $C_{{\rm S},m}$ appearing in Figure 2Ee.

###  5.4 Effective transport coefficients

Figure 4: Space-time diagrams of $(\Sigma(t)-\Sigma(0))/\Sigma(0)$ (top),
$\alpha_{\rm stress}$ (middle), and $\alpha_{\rm acc}$ (bottom), defined by
equations (16)-(17) for the four simulations described in §5. Note that the
two $\alpha$ parameters often have different signs, even though their
amplitudes vary in similar fashion (with ’bursts’ over the same time
intervals). Large values of $\alpha$ typically lead to significant depletion
of material from the disk. See text for details.

A standard way of quantifying transport processes in an accretion disc is
through the use of the dimensionless $\alpha$ parameter, which effectively
normalizes stress by the thermal pressure $\Sigma c_{s}^{2}$ (Shakura &
Sunyaev, 1973). If the stress is caused by some form of local shear viscosity
(e.g. due to MRI turbulence), then the same $\alpha$ also characterizes mass
accretion through the disk. However, the transport of mass and energy
associated with the acoustic modes is intrinsically non-local as sonic waves
are dissipated far from their launching site. Nevertheless, we can still
formally define the two dimensionless quantities based on Reynolds stress
associated with the wave angular momentum flux $C_{\rm S}$ and accretion mass
flux $\dot{M}$ as follows:

$\displaystyle\alpha_{\rm stress}$ $\displaystyle\equiv\dfrac{C_{\rm S}}{2\pi
r^{2}\Sigma c_{s}^{2}}=\dfrac{C_{\rm S}}{2\pi
r^{2}\Sigma}\left[\dfrac{\mathcal{M}}{v_{K}(R_{\star})}\right]^{2},$ (16)
$\displaystyle\alpha_{\rm acc}$ $\displaystyle\equiv\dfrac{\dot{M}}{2\pi\Sigma
c_{s}^{2}}\Omega=\dfrac{\dot{M}}{2\pi\Sigma}\left[\dfrac{\mathcal{M}}{v_{K}(R_{\star})}\right]^{2}\Omega.$
(17)

Using equation (10) we can also write $\alpha_{\rm acc}$ as

$\displaystyle\alpha_{\rm acc}$
$\displaystyle=\dfrac{\Omega}{2\pi\Sigma}\left[\dfrac{\mathcal{M}}{v_{K}(R_{\star})}\right]^{2}\left(\frac{\partial
l}{\partial r}\right)^{-1}\left[\frac{\partial C_{\rm S}}{\partial r}+2\pi
r^{3}\langle\Sigma\rangle\frac{\partial\Omega_{0}}{\partial t}\right]$ (18)
$\displaystyle=\dfrac{\Omega}{2\pi\Sigma}\left[\dfrac{\mathcal{M}}{v_{K}(R_{\star})}\right]^{2}\left(\frac{\partial
l}{\partial r}\right)^{-1}\frac{\partial C_{\rm S}}{\partial
r}~{}~{}~{}~{}~{}~{}\mbox{in steady state.}$ (19)

Any local transport mechanism acting as a source of shear viscosity would
naturally have $\alpha_{\rm stress}=\alpha_{\rm acc}$, but in our case of the
transport being global this is no longer true. We illustrate this difference
in Figure 4, where we show $\alpha_{\rm stress}$ (middle row) and $\alpha_{\rm
acc}$ (bottom row) as a function of time and radial location in the vicinity
of the stellar surface for several values of ${\mathcal{M}}$. This data has
also been box-car smoothed in the time dimension only, with a width of $\delta
t/2\pi=5$.

First, one can see that in many cases $\alpha_{\rm stress}$ and $\alpha_{\rm
acc}$ are different not only in magnitude but also in sign, in drastic
difference with the Shakura & Sunyaev (1973) model. In particular, for all
${\mathcal{M}}$ we find $\alpha_{\rm acc}$ to be positive over extended
intervals of time and space, while $\alpha_{\rm stress}$ is negative (clearly
visible for $t/2\pi\lesssim 200$, $r<1.4$ in ${\mathcal{M}}=7$ case, almost
always in the ${\mathcal{M}}=12$ case, and so on). Comparison with the Figures
2,3 shows that this situation is typical whenever the parts of the disk
adjacent to the star are dominated by the modes trapped in the resonant cavity
near the star — lower and resonant modes — which carry $C_{\rm S}<0$ and have
large $\partial{C_{\rm S}}/\partial r>0$. This naturally results in
$\alpha_{\rm stress}<0$ and $\alpha_{\rm acc}>0$, see equations (16)-(19).

Further from the star, the disk is typically dominated by the upper and
vortex-driven modes carrying low-amplitude $C_{\rm S}>0$, which results in
$\alpha_{\rm stress}>0$. As for $\alpha_{\rm acc}$, we find it to be negative
in some cases, e.g. at $t/2\pi\lesssim 100$ in ${\mathcal{M}}=12$ case, or at
$t/2\pi\sim 100$ in ${\mathcal{M}}=15$ case, when strong vortex-driven modes
are present for $r\gtrsim 1.5$. This situation is typical when $C_{\rm S}>0$
carried by the upper or vortex-driven modes decays sufficiently rapidly with
$r$ (as a result of nonlinear dissipation), leading to large $\partial{C_{\rm
S}}/\partial r<0$ and mass outflow, see (18)-(19). But more often we find
$\alpha_{\rm acc}>0$ even far from the star, despite ${C_{\rm S}}$ being
positive there (e.g. for $t/2\pi\lesssim 200$ in ${\mathcal{M}}=7$ case). This
may seem surprising, since $\partial C_{\rm S}/\partial r$ is typically small
but negative during such episodes, see Figure 3Be, which would result in mass
outflow according to the equation (19). However, this steady-state equation
does not apply in such situations, as the (second) time-dependent term in the
brackets in equation (18) is often more important than the (negative)
divergence of ${C_{\rm S}}$, driving mass inflow and $\alpha_{\rm acc}>0$ even
far from the star.

Second, both $\alpha_{\rm stress}$ and $\alpha_{\rm acc}$ are clearly highly
variable both in time and space. This has important implications for the
mediation of accretion flow by the modes, which will be discussed in §9.1. The
peak values of $\alpha_{\rm stress}$ and $\alpha_{\rm acc}$ are reached during
the relatively short-lived bursts of activity. Outside of these episodes both
$|\alpha_{\rm stress}|$ and $|\alpha_{\rm acc}|$ can maintain much lower
levels $\lesssim 10^{-3}$ for extended periods of time.

Third, with our data we can compare the values of $\alpha$ typical for
simulations with different ${\mathcal{M}}$. One can see that $\alpha_{\rm
stress}$ varies roughly from $-10^{-2}$ to $10^{-3}$, with negative (positive)
values marking the local dominance of either lower or resonant (upper or
vortex-driven) modes. At the same time, $\alpha_{\rm acc}$ only rarely attains
negative values (at the level of $-10^{-3}$) implying mass outflow. Its peak
amplitude can be as high as $\sim 0.1$ when it is positive. Results shown in
Figure 4 suggest that $\alpha_{\rm stress}$ and $\alpha_{\rm acc}$ do not show
any obvious dependence on ${\mathcal{M}}$. These transport parameters
certainly show little variation with ${\mathcal{M}}$ at higher $\dot{M}$, e.g.
during the bursts of activity confined to relatively short intervals in the
beginning of each simulation. This observation is important since it suggests
that similar values of $\alpha_{\rm stress}$ and $\alpha_{\rm acc}$ may be
expected also for significantly higher ${\mathcal{M}}$, typical for many
astrophysical objects, see Paper I.

High amplitudes of $\alpha_{\rm stress}$ and $\alpha_{\rm acc}$ during the
initial phases of our simulations lead to rapid re-adjustment of the disk
structure in response to the vigorous deposition of the angular momentum
carried by the waves, as described in §8. The top row of Figure 4 illustrates
this process by showing the space-time diagram of the local surface density
change in the disk. It is clear that injection of angular momentum by the
acoustic modes into the inner disk leads to the depletion of gas in the disk
out to $r\sim(1.2-1.6)$, which gets accreted onto the central object. The
depletion usually takes place early on, following the large bursts of
$\alpha_{\rm acc}$.

There is also an inverse effect: a significant depletion of mass in the inner
disk Figure 4 acts to suppress the mode amplitudes and to reduce the
efficiency of the angular momentum and mass transport. This is quite clear in
the ${\mathcal{M}}=7,9$ simulations, where the amplitudes of $\alpha_{\rm
stress}$ and $\alpha_{\rm acc}$ sharply drop right after a substantial gap has
been carved up in the inner disk during the initial burst of accretion. In our
purely hydrodynamic simulations, there is no mechanism (e.g. MRI) to replenish
this accreted material. As a result, we do not achieve a true steady state.
However, situation would be different when other transport mechanisms
operating in the bulk of the disk are taken into account, see §9.1.

##  6 Detailed angular momentum balance

Next we examine the role played by the different contributions in the angular
momentum balance equation (9). In Figure 5 we plot the time averages (see
Appendix A for details of the averaging procedure) of the three terms featured
in that equation, computed for $v_{\phi,0}=v_{\Sigma}$ (see equation (8)), for
three different simulations (${\mathcal{M}}=6,9,12$) listed in the figure. One
can see that in all runs the sum of the two terms in the right-hand side of
the equation (9), represented by the red dotted curve, falls right on top of
the black $\dot{M}\partial_{r}l$ curve, as expected. This demonstrates that
the relation (9) indeed holds with high accuracy in our simulations.

Figure 5: Different terms in the angular momentum balance equation (9)
extracted from our simulations. Blue ($\partial_{r}{C_{\rm S}}$), orange
(denoted as $\partial_{t}v_{\phi}$) and black ($\dot{M}\partial_{r}l$) curves
correspond to the second, third and first terms in this equation, respectively
(see the legend). Red dotted curve gives the sum of the blue and orange
curves, which should equal the black one in theory (as it does). Appropriately
averaged (see Appendix A for details) data are shown for 3 runs: (a)
${\mathcal{M}}=6$ (run M06.HR.r.a) for $t/2\pi=200-300$; (b) ${\mathcal{M}}=9$
(run M09.FR.r.a) for $t/2\pi=400-500$; (c) ${\mathcal{M}}=12$ (run
M12.FR.mix.a) for $t/2\pi=500-600$.

In most cases $\dot{M}$ (as well as $\dot{M}\partial_{r}l$, which is plotted)
is positive (meaning inflow) and largest near the star, for $r\lesssim 1.5$.
This behavior can often be directly related to that of $C_{\rm S}$, e.g. see
Figures 2Be, 3Ae,Aj, in which the positive radial derivative of $C_{\rm S}$
provides the dominant contribution to the right hand side of the equation (9)
near the star. Large $\partial_{r}C_{\rm S}>0$ leading to large $\dot{M}$
arises in the inner disk because of the deposition of the AMF carried by the
lower and resonant modes, which are trapped in the resonant cavity near the
star. On the other hand, at larger radii one can also find low-amplitude
$\dot{M}<0$, meaning outflow, see Figure 5b at $r\gtrsim 2.5$ and §5.4.

Very importantly, one can see that in all cases shown in Figure 5 the last
term in equation (9) accounting for the time-dependence of the mean rotational
profile in the disk (marked $\partial_{t}v_{\phi}$ for brevity) plays a very
significant role. In many cases its contribution (orange curve) is similar in
amplitude (but opposite in sign) to that of $\partial_{r}C_{\rm S}$, see e.g.
Figure 5a,b. At first, this may seem rather surprising since the variation of
$\Omega$ during the corresponding time intervals is quite small (see Figure
8). However, it must be also kept in mind that all terms in equation (9) are
second-order in perturbed variables and are thus small in magnitude. For that
reason, even a weak variation of $\Omega$ in time can provide a substantial
contribution to the angular momentum balance in equation (9). We will discuss
the implications of this observation in §9.1.

##  7 $C_{\rm S}$-$\dot{M}$ correlation

Figure 6: Correlation between the time-averaged values of $C_{\rm S}$ and
$\dot{M}$ for two runs: (a) ${\mathcal{M}}=9$ run M09.FR.r.a and (b)
${\mathcal{M}}=12$ run M12.FR.mix.a. Dots show the time-averaged data taken at
different moments of time, with color showing the duration of time-averaging
interval: 5 (blue) or 10 (yellow) orbits. The dotted black lines are $C_{\rm
S}=\pm 0.5\dot{M}$ to illustrate equation (20).

A classical steady-state accretion disk with no torque at the origin, which is
mediated by the local shear viscosity necessarily exhibits a linear scaling
between the instantaneous, local values of the viscous angular momentum flux
$C_{\rm S}$ and $\dot{M}$ in the form $C_{\rm S}=\dot{M}l_{K}$, where $l_{K}$
is the Keplerian angular momentum. The non-local wave-driven transport that we
observe in our BL simulations does not need to obey the same relation, for
several reasons. First, the angular momentum balance in equation (9) involves
the time-dependent term, which is generally non-zero, see §6. Second, even in
steady state equation (10) provides only a differential relation between
$C_{\rm S}$ and $\dot{M}$; solving for ${C_{\rm S}}$ would require the
knowledge of its value at the boundary $r=R_{\star}$, which is not necessarily
zero (as we show next).

Nevertheless, the overall amplitude of $\dot{M}$ should still correlate
somehow with the amplitude of $C_{\rm S}$, an expectation based on the
equation (9). To check whether this is the case we performed the measurement
of $C_{\rm S}$ and $\dot{M}$ at the stellar surface, i.e. $C_{\rm
S}(R_{\star})$ and $\dot{M}(R_{\star})$. The radius $r=R_{\star}$ is chosen so
as to minimize the impact of wave damping and non-uniformity of $\dot{M}$ in
the disc on the measurement of angular momentum and mass fluxes. This
measurement is quasi-simultaneous, since we average these quantities over $5$
or $10$ inner orbits to reduce the stochastic noise, which is naturally
present near the BL. The results are plotted in Figure 6 for the
${\mathcal{M}}=9$ and ${\mathcal{M}}=12$ runs previously discussed in
§5.2.2,5.2.3.

There are several notable features in these plots. First, $C_{\rm
S}(R_{\star})$ is almost always negative, which is consistent with the $C_{\rm
S}(r)$ behavior in Figures 2,3. Second, for moderate negative values of
$C_{\rm S}(R_{\star})$ one often observes values of $\dot{M}$, which are
roughly equal in magnitude but opposite in sign. This is most likely caused by
the oscillatory nature of the wave driven transport activity near the BL (note
that Figure 6 uses quasi-logarithmic scale).

Third, and most importantly for us, the largest negative values of $C_{\rm
S}(R_{\star})$ are very tightly correlated with the largest positive values of
$\dot{M}(R_{\star})$. We quantified this correlation by examining the
dimensionless quantity $C_{\rm S}/\dot{M}$ (in simulation units, with
$R_{\star}=1$, $v_{K}(R_{\star})=1$) in every run of our simulation suite
described in Paper I. To evaluate this ratio, we averaged $C_{\rm
S}(R_{\star})$ and $\dot{M}(R_{\star})$ separately in time over five orbits
for each simulation555We also tried using 10 orbits for averaging interval and
found that the results were consistent to within $10\%$.. Part of the reason
for this averaging666We discard some bins to remove the initial data and noisy
data in which both quantities frequently change sign. is to smooth over
possible time lags between $C_{\rm S}(R_{\star})$ and $\dot{M}(R_{\star})$
arising due to the nonlocal nature of the wave driven transport. This data is
plotted in Figure 7, and results in (median $\pm$ one quartile)

$\displaystyle\left.\dfrac{C_{\rm
S}}{\dot{M}\ell_{K}}\right|_{r=R_{\star}}=-0.51^{+0.28}_{-0.17}.$ (20)

This relation (shown with dotted lines in Figures 6,7) once again demonstrates
that mass inflow into the BL requires negative $C_{\rm S}$ at the stellar
interface. It also characterizes the efficiency with which mass accretion is
enabled by the angular momentum injection into the disk by the wave modes
excited in the BL. Finally, it provides an important boundary condition, which
can be used in constructing semi-analytical models describing the structure of
the inner regions of accretion disks fully or partly mediated by the acoustic
waves.

Figure 7: The ratio of $C_{\rm S}/\dot{M}$ at $r=R_{\star}=1$ (in simulation
units, equal to $C_{\rm S}/(\dot{M}\ell_{K})|_{r=R_{\star}}$ in physical
units) as a function of $\dot{M}$, with the color denoting the value of
$\mathcal{M}$ of the runs from the data were taken. Each colored point
represents a sample of five orbit period averages of data for the
corresponding ${\mathcal{M}}$, such as blue dots shown in Fig. 6, and the
statistical properties of that sample: in panel (a) we show means $\pm$ one
standard-deviation and in (b) we show medians $\pm$ one quartile. The data are
binned in $\dot{M}$ prior to constructing the sample, the vertical dotted grey
lines indicate the bins. Horizontal placement of points within a bin is not
meaningful; the points are spread out so that they do not visually overlap.
Points without error bars correspond to samples consisting of a single five
orbit period average. The horizontal black dotted line and shaded grey band
show the global mean/median $\pm$ one standard-deviation/quartile for (a)/(b)
respectively.

##  8 Wave-driven disk evolution

Figure 8: Evolution of $\Omega(r)$ profiles in simulations with different
values of $M$ (found in the run label). Different curves correspond to
different moments of time coded in the colorbar on the right. See text for
details. Figure 9: (a) Width of the BL as a function of ${\mathcal{M}}$;
dotted line corresponds to equation (22). (b) Deviation of the plateau from
$\Omega_{K}(R_{\star})$ as a function of ${\mathcal{M}}$; dotted line
corresponds to equation (23). (c) Radial extent of the $\Omega$ plateau as a
function of ${\mathcal{M}}$; dotted line corresponds to equation (24). The
dotted lines are by-eye fits to the relevant data. The data shown correspond
to the final output file of each simulation ($t/2\pi=600$).

In the absence of wave damping $C_{\rm S}$ would be exactly conserved in our
globally isothermal setup, implying no transfer of the wave angular momentum
to the disk and $\dot{M}=0$. But in our simulations wave modes of different
types can exchange their angular momentum with the disk fluid, causing its
surface density to evolve. Since our runs have no viscous or radiative
dissipation, this exchange is accomplished through the nonlinear dissipation
of the waves (Goodman & Rafikov, 2001; Rafikov, 2002): nonlinear steepening of
the acoustic wave profile eventually turns the wave into a shock (Landau &
Lifshitz, 1959), transferring its angular momentum to the background flow
(Rafikov, 2016). This process inevitably occurs even if the shock is weak.

Nonlinear damping always reduces wave amplitude and $|C_{\rm S}|$, regardless
of the sign of $C_{\rm S}$. As a result, in agreement with the equations
(9)-(10), non-zero $\dot{M}$ arises in the disk changing its structure. But
the direction of the mass flow also depends on a particular type of the wave
and, more specifically, on the sign of AMF $C_{\rm S}$ that it carries.

As we saw in §5.4, damping of the modes with $C_{\rm S}<0$, e.g. lower and
resonant modes (see §5.2), gives rise to $\partial C_{\rm S}/\partial
l\propto\partial C_{\rm S}/\partial r>0$, leading to $\dot{M}>0$ (see
equations (9),(10)), i.e. mass inflow through the disk. On the other hand,
dissipation of the modes with $C_{\rm S}>0$, i.e. the upper and vortex-driven
modes, can sometimes lead to $\partial C_{\rm S}/\partial l>0$. As a result,
the disk fluid gains angular momentum, resulting in $\dot{M}<0$, i.e. mass
outflow.

In practice, lower and resonant modes carrying $C_{\rm S}<0$ generally affect
disk evolution much stronger than the other types of modes. Since these modes
are trapped in the resonant cavity in the inner disk, gas accretion onto the
star driven by them clears out a substantial depression, or gap, in the inner
disk. The development of such a gap is illustrated in the top row of Figure 4.
One can see that the reduction of $\Sigma$ can be quite substantial, with the
gap depth reaching $\sim 80\%$ of the original density at $r\lesssim 1.5$ in
simulations with ${\mathcal{M}}\lesssim 9$. For higher-${\mathcal{M}}$ runs,
which often exhibit more regular patterns of the acoustic modes, the decrement
of $\Sigma$ is not so large, with the gap depths at the level of
$\sim(20-30)\%$ being more typical, see Figure 4. The radial extent of the gap
also decreases with increasing ${\mathcal{M}}$, in agreement with the smaller
radial width of the resonant cavity at higher ${\mathcal{M}}$, see §5.2.

In our simulations the deepening of the gap saturates once it reaches a
sufficient depth that depends on ${\mathcal{M}}$. The gap does not get
refilled by the matter arriving from larger radii, as would be expected in a
real accretion disk, because our runs do not have an an explicit mechanism of
the angular momentum transport in the outer disk (i.e. $\alpha$-viscosity or
MRI). In real disks one should expect the gap to be shallower than what we
find in our simulations.

###  8.1 Evolution of the $\Omega(r)$ profile

Sharp gradients of $\Sigma$ that develop around the gap region cause
substantial modification of $\Omega(r)$ profile away from the purely Keplerian
$\Omega_{K}(r)$ in the bulk of the disk. In a steady state $\Omega(r)$ is
given by

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

following from the radial momentum balance equation. Because of the initial
non-uniform profile of $\Sigma$ and the resultant radial pressure gradient,
some deviation of $\Omega$ from $\Omega_{K}$ is present even at the start of
simulations.

But after the gap develops in the disk, these $\Omega$ deviations increase as
the second term in the right hand side of the equation (21) becomes much
larger. This is illustrated in Figure 8, which shows the profiles of
$\Omega(r)$ at different moments of time in simulations with different
${\mathcal{M}}$. Using the simulation data on $\Sigma(r)$ (e.g. the ones shown
in the top row of Figure 4) we find that the $\Omega(r)$ profile agrees very
well with the equation (21). Interestingly, equation (21) describes the
$\Omega(r)$ behavior quite well even inside the BL, where $\Omega(r)$
substantially deviates from $\Omega_{K}$. This means that even within the BL
the radial velocity of the accreting matter $v_{r}$ is still quite small, so
that the inertial terms in the radial momentum balance equation can be
neglected.

In agreement with the equation (21), we find $\Omega<\Omega_{K}$ interior from
the deepest part of the gap, where $\Sigma$ and $P$ decrease with radius. In
fact, $\Omega$ tends to exhibit a relatively flat, plateau-like segment with
$\Omega(r)\approx\Omega_{\rm max}$ just outside the BL. Right outside the
deepest part of the gap, $P$ increases with $r$ (very strongly for lower
${\mathcal{M}}$), driving $\Omega$ above $\Omega_{K}$. This can be most easily
seen in Figure 8c-f. As the gap width and depth become smaller for higher
${\mathcal{M}}$, so does the deviation of $\Omega$ from $\Omega_{K}$: the
region near the BL where $\Omega(r)$ exhibits a plateau gets narrower as
${\mathcal{M}}$ grows, and the value of $\Omega$ in this region gets closer to
$\Omega_{K}$.

To characterize these $\Omega(r)$ features developing near the BL, in Figure 9
we plot the width of the BL $\delta_{\rm BL}$, defined as the radial extent
over which $\Omega$ transitions from 25% to 75% of its maximum value, i.e.
$\delta_{\rm BL}\equiv R\left(\Omega=0.75\Omega_{\rm
max}\right)-R\left(\Omega=0.25\Omega_{\rm max}\right)$ (panel (a)); the
deviation of $\Omega(r)$ in the plateau region $\delta\Omega$ from
$\Omega_{K}(R_{\star})=1$, i.e. $\delta\Omega\equiv 1-\Omega_{\rm max}$ (panel
(b)); and the radial width $\delta_{\rm plateau}$ of the plateau in $\Omega$,
defined as the radial extent over which $\Omega\geq 0.9\Omega_{\rm max}$. One
can see that all three plotted variables generally decrease with increasing
${\mathcal{M}}$, as mentioned previously. Note that some values of
${\mathcal{M}}$ for which we have multiple simulations show a spread in these
derived variables, see e.g. the variation of $\delta_{\rm BL}$ for
${\mathcal{M}}=9$ in panel (a).

Such spreads in $\delta_{\rm BL}$ are caused by the emergence of an inflection
point in $\Omega(r)$ profile inside the BL in some runs, which acts to
increase the radial range over which $\Omega$ varies. This feature is clearly
seen e.g. in Figure 8d-f (see Belyaev et al. 2012 and Dittmann 2021 for
similar observations). To circumvent the effect of this feature on the
determination of $\delta_{\rm BL}$ we focus on the smallest values of
$\delta_{\rm BL}$ for a given ${\mathcal{M}}$. To guide the eye, in Figure 9a
we run777This and other dotted lines in Figure 9 are not formal fits to the
data. a dotted line

$\displaystyle\delta_{\rm BL}\approx 2~{}{\mathcal{M}}^{-2}$ (22)

through these points that appears to provide a decent fit to the lower
envelope of the measured $\delta_{\rm BL}$ values. As previously noted in
Belyaev et al. (2012), this scaling implies that $\delta_{\rm BL}$ is equal to
a fixed (i.e. independent of ${\mathcal{M}}$) number of stellar scaleheights
(equal to ${\mathcal{M}}^{-2}$ in our notation). In our case this number is
$\approx 2$, which is significantly less than $7-8$ found in Belyaev et al.
(2012), a difference naturally explained by the fact that Belyaev et al.
(2012) did not attempt to account for the effect of the inflection point on
the $\Omega$ profile.

In Figure 9b we plot $\delta\Omega$ — the deviation of $\Omega$ from unity in
the plateau region — together with a by-eye fit

$\displaystyle\delta\Omega\approx 1.5~{}{\mathcal{M}}^{-1}.$ (23)

This rough scaling can be understood analytically using equation (21) if we
also assume that the radial scale of $\Sigma$ variation in the gap near the
stellar surface is of order of the disk scaleheight ${\mathcal{M}}^{-1}$,
which appears to be the case in our runs. The data for the radial extent of
$\Omega$ plateau $\delta_{\mathrm{plateau}}$ is shown in Figure 9c together
with a simple fit (dotted curve)

$\displaystyle\delta_{\mathrm{plateau}}\approx 0.9~{}{\mathcal{M}}^{-2/3}.$
(24)

Based on the Keplerian rotation law and equation (23), one might have expected
a steeper dependence on ${\mathcal{M}}$ for a flat $\Omega$ plateau, closer to
$\propto{\mathcal{M}}^{-1}$. However, examination of Figure 8 shows that
approximating the $\Omega(r)$ plateau shape as flat is very approximate, and
that it has a certain non-trivial radial dependence to it. Thus, the exponent
in equation (24) different from $-1$ should not come as a surprise.

Knowledge of $\delta\Omega$ and $\delta_{\mathrm{plateau}}$ is important since
we often find the plateau region to contain near-surface vortices that drive
prominent spiral arms in the disc, see Paper I. The pattern speed of such
features should be roughly equal to $1-\delta\Omega$ (as the vortices are
comoving with the fluid), while $\delta_{\mathrm{plateau}}$ informs us about
the distance out to which these near-surface vortices might be expected to be
found.

Given that the exact shape of the $\Omega$ plateau depends on the $\Sigma(r)$
profile in the gap, one may wonder if scalings (23), (24) would still hold in
the presence of mass inflow from larger radii, which is expected in real
disks. It is very likely that $\delta\Omega$ and $\delta_{\rm plateau}$
behaviors shown in Figure 9 will change in that situation. However, the
scaling (22) for the width of the BL ($\delta_{\rm BL}$) must be robust
regardless of the presence or absence of the mass inflow, since the existence
of the BL is independent of the presence or absence of a depression in
$\Sigma$ outside the star.

##  9 Discussion

Results of the previous sections illustrate the important effect of the wave
modes triggered by the supersonic shear inside the BL on the angular momentum
and mass transport in the inner disk. We now discuss some additional aspects
of the wave-driven transport processes in the vicinity of the BL.

###  9.1 Modes dominating mass transport and disk evolution

The discussion in §5.3.1-5.3.5 allows us to identify the modes most relevant
for driving the accretion onto the star. Since accretion implies $\dot{M}>0$,
we can immediately conclude that it is the lower acoustic and resonant modes
that must be responsible for driving the mass inflow onto the stellar surface.
This conclusion is supported by Figure 6 and the results presented in §7,
which show that the highest values of $\dot{M}$ require ${C_{\rm S}}<0$ ,
which is a unique feature of these two types of modes.

At the same time, the vortex-driven (and, at lower amplitude, the upper
acoustic) modes drive a weak mass outflow in the disk, which is not conducive
to accretion. In realistic disks this outflow would be suppressed by the mass
inflow caused by other sources of the angular momentum transport capable of
operating in the bulk of the disk, e.g. the MRI.

It is important that the $\dot{M}$ driven by the lower acoustic and resonant
modes is large only in radially narrow region of the disk near the stellar
surface, where the mode is trapped. This means that these modes by themselves
are not capable of maintaining constant $\dot{M}$ through the whole disk,
which would be necessary to keep it in a steady state. Other transport
mechanisms are clearly needed to ensure a steady delivery of mass towards the
BL, e.g. the MRI.

It is also important to realize that in realistic disks, inside the resonant
cavity where the $\dot{M}$-driving modes are trapped, the angular momentum
transport would be provided both by the waves, globally, and by the MRI,
locally (since $d\Omega/dr<0$ in that region, see Figure 8). This enhanced
transport is likely to result in a non-trivial structure of the surface
density in this near-BL part of the disk.

###  9.2 Global nature of the angular momentum transport

Our morphological study in Paper I demonstrates that the waves excited in the
vicinity of the BL can propagate over substantial distances in the disk. This
is clear for the vortex-driven and upper acoustic modes that can propagate
from their corotation region near the BL all the way into the outer disk where
they get gradually damped. However, even the modes which are trapped near the
BL — the lower acoustic waves and the resonant modes — still travel over a
significant radial range in the disk, comparable to $R_{\star}$. This is
especially true for the low-$m$ resonant modes, for which the corotation
radius can lie quite far out in the disk, see Paper I.

The long-range propagation of the waves makes the angular momentum transport
effected by them truly global. It gives rise to a non-local coupling between
the BL, where the accreting matter uploads its angular momentum and energy to
the waves, and the more distant parts of the disk, where the waves dissipate,
driving $\dot{M}$ and injecting their energy into the disk fluid (Rafikov,
2016). Such non-locality is typical also for the angular momentum transport by
the planet-driven density waves in the protoplanetary disks (Lunine &
Stevenson, 1982; Goodman & Rafikov, 2001; Rafikov, 2002; Rafikov & Petrovich,
2012; Petrovich & Rafikov, 2012).

Because of its non-locality, the wave-driven transport cannot be characterized
by a single $\alpha$ parameter, which is often invoked to describe local
transport (Shakura & Sunyaev, 1973; Balbus & Papaloizou, 1999). As we
demonstrate in §5.4, the values of $\alpha$ near the BL computed using stress
${C_{\rm S}}$ and $\dot{M}$ (which are the same for local transport) are
different not only in magnitude but also in sign, see Figure 4. For that
reason, the standard methods for calculating the structure and viscous
evolution of the disk based on the diffusion-type equation approach (Lynden-
Bell & Pringle, 1974), would fail to describe the disk in the vicinity of the
BL. Instead, calculations of the wave-mediated structure of the inner disk
need to explicitly account for the long-range propagation and dissipation of
the waves excited in the BL.

###  9.3 Observational implications

Non-local wave-driven transport effected by the waves has important
implications for the observational appearance of the objects accreting through
the BL. Models of the BL structure using the local $\alpha$-viscosity
inevitably predict that local energy dissipation in the layer should heat gas
in the BL to high temperatures (Popham et al., 1993; Narayan & Popham, 1993).
This would naturally give rise to a high-energy component in the spectra of
objects accreting through the BL, which is often not observed (Ferland et al.,
1982).

Energy transport by the waves changes the whole pattern of the energy
dissipation: it occurs not in the narrow BL itself, but over a larger area of
the disk, lowering the effective temperature of emission associated with the
BL and naturally alleviating the "missing boundary layer" problem (Ferland et
al., 1982). The rate $d\dot{E}/dr$ (per unit time and radius) at which thermal
energy is deposited in the disk by a wave (with pattern speed $\Omega_{P}$) is
closely related to the rate of dissipation of the wave angular momentum flux,
as $d\dot{E}/dr=(\Omega-\Omega_{P})d{C_{\rm S}}/dr$ (Lynden-Bell & Kalnajs,
1972; Goodman & Rafikov, 2001; Rafikov, 2016). Thus, the results of our study
on the radial profiles of ${C_{\rm S}}$ for different wave modes can be
directly employed to analyze the spatial distribution of the disk heating
induced by the BL-excited waves.

###  9.4 Comparison with the existing studies

In the past, Belyaev et al. (2012, 2013a, 2013b) analyzed the behavior of the
wave angular momentum flux ${C_{\rm S}}$ in the vicinity of the BL, and
related it to the $\dot{M}$ behavior and the wave amplitude in the disk. More
recently, Dittmann (2021) looked at the ${C_{\rm S}}$ behavior in the disk (as
well as the flow of mass and angular momentum into the star) for non-zero
stellar spin rates. Our present work goes beyond these studies in several
ways.

First, we analyze ${C_{\rm S}}$ behavior for a variety of individual modes
present in the disk. Coupling this with the results of our morphological study
in Paper I allows us to understand the ${C_{\rm S}}$ behavior for each mode
type — lower and upper, resonant, vortex-driven, and so on. Second, we follow
how the behavior of ${C_{\rm S}}$ changes as our simulations progress and
different types of wave modes come and go. Third, we systematically explore
the differences in the ${C_{\rm S}}$ behavior as a function of the Mach number
${\mathcal{M}}$ of our runs.

##  10 Summary

In this work we studied the angular momentum and mass transport driven by the
different types of the waves operating in the vicinity of the BL of an
accretion disk. Our analysis is based on a large suite of global, 2D
hydrodynamic simulations described in Paper I, in which a large number of
modes have been previously identified and characterized. Results of our
present study can be summarized as follows.

  * •

In wave-mediated inner region of the disk mass accretion rate $\dot{M}$ has a
significant contribution arising from the correlated radial velocity $v_{r}$
and the surface density perturbation $\delta\Sigma$.

  * •

The efficiency of angular momentum transport expressed through the effective
$\alpha$-parameter appears to not depend strongly on the Mach number of the
flow ${\mathcal{M}}$.

  * •

By examining angular momentum flux $C_{{\rm S},m}$ carried by each individual
Fourier harmonics, we were able to determine the transport properties
associated with each type of the mode for different values of the Mach number
${\mathcal{M}}$. In particular, we find that lower acoustic and resonant modes
carry negative angular momentum flux, whereas the vortex-driven and upper
acoustic modes carry positive angular momentum flux.

  * •

Nonlinear damping of the modes leads to mass inflow ($\dot{M}>0$) for the
lower acoustic and resonant modes, and weak mass outflow for the vortex-driven
and upper acoustic modes. This implies that accretion onto the central object
must be mediated by a combination of the lower acoustic and resonant modes.

  * •

Wave-driven transport is intrinsically non-local, which leads to the values of
effective $\alpha$-parameter determined through stress and $\dot{M}$ being
different not only in amplitude but also in sign. This is very different from
the conventional local transport for which these values are the same.

  * •

Despite the non-local nature of the wave-mediated transport we still find a
strong correlation, given by the equation (20), between the angular momentum
flux injected at the BL into the disk, and the mass accretion rate through the
BL.

  * •

Despite the long duration of our simulations, we find the time-dependent
contribution to the angular momentum balance caused by the variation of
$\Omega$ (see equation (9)) to play an important role (see §6).

  * •

We characterize the wave-driven evolution of the disk properties — surface
density $\Sigma$ and angular frequency $\Omega$ profile — as a function of
time and ${\mathcal{M}}$.

Our study, based on 2D hydro simulations, provides a natural foundation for
understanding the properties of the wave-driven transport in future three-
dimensional and fully MHD simulations of the BLs.

## Acknowledgements

We thank Jim Stone for making the Athena++ 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).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Arzamasskiy & Rafikov (2018) Arzamasskiy L., Rafikov R. R., 2018, The Astrophysical Journal, 854, 84
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Reviews of Modern Physics, 70, 1
  * Balbus & Papaloizou (1999) Balbus S. A., Papaloizou J. C. B., 1999, ApJ, 521, 650
  * Balsara et al. (2009) Balsara D. S., Fisker J. L., Godon P., Sion E. M., 2009, ApJ, 702, 1536
  * 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)
  * Coleman et al. (2021) Coleman M. S. B., Rafikov R. R., Philippov A. A., 2021, MNRAS,
  * Dittmann (2021) Dittmann A. J., 2021, MNRAS, 508, 1842
  * Ferland et al. (1982) Ferland G. J., Langer S. H., MacDonald J., Pepper G. H., Shaviv G., Truran J. W., 1982, ApJ, 262, L53
  * Flock et al. (2011) Flock M., Dzyurkevich N., Klahr H., Turner N. J., Henning T., 2011, ApJ, 735, 122
  * Fromang & Nelson (2006) Fromang S., Nelson R. P., 2006, A&A, 457, 343
  * Ghosh & Lamb (1978) Ghosh P., Lamb F. K., 1978, ApJ, 223, L83
  * Gnedin et al. (1995) Gnedin O. Y., Goodman J., Frei Z., 1995, AJ, 110, 1105
  * Goodman & Rafikov (2001) Goodman J., Rafikov R. R., 2001, ApJ, 552, 793
  * Hertfelder & Kley (2017) Hertfelder M., Kley W., 2017, A&A, 605, A24
  * Ju et al. (2016) Ju W., Stone J. M., Zhu Z., 2016, ApJ, 823, 81
  * Koenigl (1991) Koenigl A., 1991, ApJ, 370, L39
  * Landau & Lifshitz (1959) Landau L. D., Lifshitz E. M., 1959, Fluid mechanics. Pergamon Press
  * Larson (1984) Larson R. B., 1984, MNRAS, 206, 197
  * Larson (1990) Larson R. B., 1990, MNRAS, 243, 588
  * Lunine & Stevenson (1982) Lunine J. I., Stevenson D. J., 1982, Icarus, 52, 14
  * Lynden-Bell & Kalnajs (1972) Lynden-Bell D., Kalnajs A. J., 1972, MNRAS, 157, 1
  * Lynden-Bell & Pringle (1974) Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
  * Narayan & Popham (1993) Narayan R., Popham R., 1993, Nature, 362, 820
  * Pessah & Chan (2012) Pessah M. E., Chan C.-k., 2012, ApJ, 751, 48
  * Petrovich & Rafikov (2012) Petrovich C., Rafikov R. R., 2012, ApJ, 758, 33
  * 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
  * Rafikov (2016) Rafikov R. R., 2016, ApJ, 831, 122
  * Rafikov & Petrovich (2012) Rafikov R. R., Petrovich C., 2012, ApJ, 747, 24
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Spruit (1987) Spruit H. C., 1987, A&A, 184, 173
  * Stone et al. (2020) Stone J. M., Tomida K., White C. J., Felker K. G., 2020, ApJS, 249, 4
  * Towns et al. (2014) Towns J., et al., 2014, Computing in Science & Engineering, 16, 62
  * Velikhov (1959) Velikhov E. P., 1959, Journal of Experimental and Theoretical Physics, 36, 995

##  Appendix A Details of the time-averaging procedures

Here we summarize the technical details of various averaging and smoothing
procedures that were used in analyzing the data and producing various plots.
For all the line plots in this paper (except for Fig. 8) we used the unbinned
FFT data to measure azimuthal averages ($m=0$; all the quantities inside
angular brackets) and higher order $m>0$ modes with high cadence (exactly 20
outputs per orbit with non uniform cadence; see Paper I for more details). As
the final step of the calculation we average in time over the intervals
denoted in the figures. For Figs. 2 and 3 we chose to plot only the five modes
with the largest time (over the interval noted in the figure) and radially
(over $r>1$) integrated $|C_{{\rm s},m}|$.

We examine Eqn. (9) in Fig. 5. We rewrite and expand this equation below to be
explicit about the averaging done in the figure:

$\displaystyle\underbrace{2\pi r\left\langle\Sigma
v_{r}\right\rangle\dfrac{\partial}{\partial r}\left(r\dfrac{\left\langle\Sigma
v_{\phi}\right\rangle}{\left\langle\Sigma\right\rangle}\right)}_{\dot{M}\partial_{r}\ell}$
$\displaystyle=\underbrace{\frac{\partial C_{\rm S}}{\partial
r}}_{\partial_{r}C_{\rm S}}+\underbrace{2\pi
r^{2}\langle\Sigma\rangle\dfrac{\partial}{\partial t}\dfrac{\left\langle\Sigma
v_{\phi}\right\rangle}{\left\langle\Sigma\right\rangle}}_{``\partial_{t}v_{\phi}"},$

where $C_{\rm S}$ is computed using Eqns. (6) and (8), and the underbraces are
used to denote how the terms are labeled in Fig. 5. All the azimuthally
averaged quantities (i.e. the ones in angular brackets) are measure directly
from the $m=0$ mode of our high cadence FFT data as explained above. The
partial derivative with respect to time is computed via finite differences on
the high cadence FFT data to second order accuracy. After these computations
are completed the data is averaged in time over the time interval shown in the
figure.