# Investigating protoplanetary disc cooling through kinematics: analytical GI
wiggle

Cristiano Longarini Dipartimento di Fisica, Università degli Studi di Milano,
via Celoria 16, 20133 Milano, Italy  Giuseppe Lodato Dipartimento di Fisica,
Università degli Studi di Milano, via Celoria 16, 20133 Milano, Italy  Claudia
Toci Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16,
20133 Milano, Italy  Benedetta Veronesi Dipartimento di Fisica, Università
degli Studi di Milano, via Celoria 16, 20133 Milano, Italy  Univ Lyon, Univ
Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574,
F-69230, Saint-Genis,-Laval, France.  Cassandra Hall Department of Physics and
Astronomy, The University of Georgia, Athens, GA 30602, USA  Center for
Simulational Physics, The University of Georgia, Athens, GA 30602, USA
Ruobing Dong Department of Physics & Astronomy, University of Victoria, 3800
Finnerty Road, Victoria, BC V8P 5C2, Canada  Jason Patrick Terry Department of
Physics and Astronomy, The University of Georgia, Athens, GA 30602, USA
Center for Simulational Physics, The University of Georgia, Athens, GA 30602,
USA

###### Abstract

It is likely that young protostellar discs undergo a self-gravitating phase.
Such systems are characterised by the presence of a spiral pattern that can be
either in a quasi-steady state or in a non-linear unstable condition. This
spiral wave affects both the gas dynamics and kinematics, resulting in
deviations from the Keplerian rotation. Recently, a lot of attention has been
devoted to kinematic studies of planet forming environments, and we are now
able to measure even small perturbations of velocity field ($\lesssim 1\%$ of
the Keplerian speed) thanks to high spatial and spectral resolution
observations of protostellar discs. In this work, we investigate the kinematic
signatures of gravitational instability: we perform an analytical study of the
linear response of a self-gravitating disc to a spiral-like perturbation,
focusing our attention on the velocity field perturbations. We show that
unstable discs have clear kinematic imprints into the gas component across the
entire disc extent, due to the GI spiral wave perturbation, resulting in
deviations from Keplerian rotation. The shape of these signatures depends on
several parameters, but they are significantly affected by the cooling factor:
by detecting these features, we can put constraints on protoplanetary discs
cooling.

Hydrodynamics - accretion, Accretion discs - Protoplanetary discs -
Instabilities - Methods: analytical

††facilities: SPLASH: an interactive visualisation tool for SPH data (Price,
2007)

##  1 Introduction

For cold massive discs, the role of the disc self-gravity becomes dynamically
important, affecting the vertical and radial structure of the system (Kratter
& Lodato, 2016; Bertin & Lodato, 1999). In this context, gravitational
instabilities (hereafter referred as GI) often arise, determining the
evolution of the system and playing a fundamental role in the transport of
angular momentum. The development of the GI has initially been studied in the
context of galactic dynamics (Lin & Shu, 1964; Bertin & Lin, 1996; Binney &
Tremaine, 1987): as far as protostellar discs are concerned, the results are
quantitatively similar.

On the one hand, a possible outcome of GI in protostellar environments is the
fragmentation of the disc: this phenomenon can potentially lead to the
formation of low-mass stellar companions (Kratter & Matzner, 2006; Stamatellos
et al., 2007; Cadman et al., 2020), because the initial clump mass is of the
order of several Jupiter masses, too high to form a planet (Kratter & Lodato,
2016). On the other hand, GI is a very effective way to transport angular
momentum within the disc, by means of a global spiral perturbation (Lodato &
Rice, 2004, 2005).

High resolution observations with ALMA have revealed that most of the observed
protostellar discs possess substructures as rings or spirals. The origin of
rings is often explained by planets (Dipierro et al., 2015, 2018; Toci et al.,
2020b; Veronesi et al., 2020), however what causes the spirals is still
ambiguous. Indeed, super-Jupiter objects can excite spiral density waves with
azimuthal wavenumber $m\sim 1-2$ that match with good agreement the observed
structures in scattered light (Dong et al., 2015; Dong & Fung, 2017; Veronesi
et al., 2019; Rosotti et al., 2020). In addition, some spirals may also be
induced by a inner or outer stellar companion (Price et al., 2018a), or by a
flyby (Cuello et al., 2019, 2020). At the same time, large scale spiral
perturbations also characterise self-gravitating discs, with a typical $m\sim
M_{\star}/M_{d}$, where $M_{d}$ is the disc mass and $M_{\star}$ is the mass
of the star (Cossins et al., 2009). Distinguishing the origin of a spiral is
difficult, but recent high resolution observations of protostellar
environments allow us to conduct kinematic studies that might shed some light
on this issue. It is well known that the presence of a perturber inside the
disc creates a localised deviation from the Keplerian observed velocity,
called “kink” (Pinte et al., 2018; Teague et al., 2018): when the perturber is
a planet, the kink can be used as a proxy for its mass (Bollati et al., 2021).
As far as GI is concerned, Hall et al. (2020) (hereafter referred as H20)
show, based on hydrodynamical simulations, that the spiral perturbation deeply
affects the gas kinematic: in particular, it creates a global (rather than a
localised) deviation from Keplerian observed velocity - a “global kink” -
dubbed GI wiggle by H20, that is apparent in the moment one and in the channel
maps (Paneque-Carreño et al., 2021).

In this work, we present an analytical study of the response of a self-
gravitating protostellar disc to a spiral density wave in WKB regime. We focus
our attention on the velocity field perturbations (hereafter referred as VPs),
and we show that they can be written as a function of disc parameters. Thanks
to such analytical expressions, we are able to sketch the observed velocity
field, and then to make a connection to observations. The amplitude of the VPs
is linked to the cooling factor of the system, and thus we can use this
relation to investigate the cooling in protostellar discs.

The paper is organised as follows. In section 2 we summarise the theory that
we need to conduct our study, paying attention to both dynamical and
thermodynamical processes. In section 3 we obtain the analytical expression of
the VPs in WKB regime. In section 4 we discuss the observational aspects of
the wiggles, projecting the perturbed velocity field along the line of sight.
To conclude, in section 5 we discuss the shape of the perturbations, the
limitations of this work and we present a mock experiment to test our
predictions.

##  2 Self-gravitating gaseous discs

###  2.1 The disc potential

The gravitational potential of a self-gravitating protostellar disc is not
exactly Keplerian: indeed, when $M_{d}\gtrsim(H/r)M_{\star}$, where $H/r$ is
the aspect ratio of the disc, the disc contribution to the potential is not
negligible. In this case, the radial balance of forces for a cold disc is
given by

$\Omega^{2}=\frac{GM_{\star}}{r^{3}}+\frac{1}{r}\frac{\partial\Phi_{d}}{\partial
r},$ (1)

where the first term is the Keplerian frequency and the second one is the disc
contribution (Bertin & Lodato, 1999). The gravitational field generated by the
disc can be written as

$\frac{\partial\Phi_{d}}{\partial
r}(r,z)=\frac{G}{r}\int_{0}^{\infty}dr^{\prime}\left[K(\zeta)-\frac{1}{4}\left(\frac{\zeta^{2}}{1-\zeta^{2}}\right)\times\left(\frac{r^{\prime}}{r}-\frac{r}{r^{\prime}}+\frac{z^{2}}{rr^{\prime}}\right)E(\zeta)\right]\sqrt{\frac{r^{\prime}}{r}}\zeta\Sigma\left(r^{\prime}\right),$
(2)

where $E(\zeta),$ $K(\zeta)$ are the complete elliptic integrals of the first
kind, $\zeta^{2}=4rr^{\prime}/[(r+r^{\prime})^{2}+z^{2}]$ and $\Sigma$ is the
disc surface density. Deviations from the Keplerian behaviour have actually
been seen in the rotation curve of Elias 2-27 (Veronesi et al., 2021), and
this can be an effective method for measuring the disc mass.

###  2.2 Gravitational instability

The linear response to gravitational instability is described by the well
known dispersion relation (Lin & Shu, 1964)

$D(\omega,k,m)=(\omega-m\Omega)^{2}-c_{\mathrm{s}}^{2}k^{2}+2\pi
G\Sigma|k|-\kappa^{2}=0,$ (3)

where $\omega$ is the wave angular frequency, $k$ the radial wavenumber, $m$
the azimuthal wavenumber, $\kappa(r)$ the epicyclic frequency, $\Omega(r)$ the
angular frequency and $c_{s}(r)$ the sound speed. This relationship has been
obtained for an infinitesimally thin disc under the WKB perturbation analysis
and in the tight winding limit (i.e. the radial wavelength is much smaller
than the azimuthal one). The stability criterion can be expressed by means of
the $Q$ parameter

$Q=\frac{c_{s}\kappa}{\pi G\Sigma},$ (4)

that contains the stabilising terms at numerator and the unstable ones at
denominator. The instability threshold is given by $Q=1$: if $Q>1$ the disc is
stable at all wavelengths while $Q<1$ identifies a range of unstable
wavelengths. For the case of an unstable disc, the most unstable wavenumber is

$k_{\text{uns}}=H_{\text{sg}}^{-1}=\frac{\pi G\Sigma}{c_{s}^{2}},$ (5)

where $H_{\text{sg}}=c_{s}^{2}/\pi G\Sigma$ is the disc thickness in the self
gravitating case.

A marginally stable disc is characterised by having $Q\simeq 1$: in this case,
the only expected excited modes would be $k\simeq k_{\text{uns}}$. Thus
equation (3) tells that $\Omega\simeq\omega/m=\Omega_{\text{p}}$, where
$\Omega_{\text{p}}$ is the spiral pattern frequency, meaning that all excited
modes are expected to be close to corotation (Cossins et al., 2009).

####  2.2.1 Spiral density waves

A spiral wave is characterised by having

$m\phi+\psi(r)=\text{const},\quad\text{mod }2\pi,$ (6)

where $\psi(r)$ is the shape function and it holds that
$k=\text{d}\psi/\text{d}r$. It is useful to introduce the radial wavenumber
$k$, that is the radial derivative of the shape function; the sign on $k$
determines wheter the spiral wave is leading ($k<0$) or trailing ($k>0$).
Another important quantity is the opening angle of the spiral $\alpha_{p}$,
hereafter referred as the pitch angle: it is given by

$\tan\alpha_{p}=\left|r\frac{\partial\phi}{\partial r}\right|^{-1},$ (7)

where the partial derivative is evaluated along (6), giving
$\tan\alpha_{p}=m/rk$.

If we consider a self-consistent spiral perturbation, we can easily link the
density perturbation to the potential one (Cossins et al., 2009). The
perturbed surface density can be written as

$\Sigma_{1}(r,\phi,t)=\text{Re}\left[\delta\Sigma(r)e^{i(m\phi-\omega
t+\psi(r))}\right],$ (8)

and it can be shown that the corresponding perturbed potential is given by

$\Phi_{1}(r,\phi,t)=-\frac{2\pi G}{|k|}\Sigma_{1}(r,\phi,t).$ (9)

####  2.2.2 Thermodynamics

So far, we have discussed only the linear growth of gravitational instability,
however, to properly understand the outcome of this process, we need to
consider the non-linear evolution. To do so, it is necessary to introduce a
parameter to capture the radiative processes of the disc.

The non-linear growth of perturbations is best understood by using numerical
simulations (Rice et al., 2005; Kratter et al., 2010; Hall et al., 2017).
However, reproducing realistically the cooling processes is not an easy task:
in the last years, a lot of effort was devoted to the modeling of realistic
thermodynamics (Johnson & Gammie, 2003; Stamatellos & Whitworth, 2009; Hirose
& Shi, 2017).

In this paper, we are not interested in the physics of cooling, but in the
relationship between the density perturbations and the rate at which the disc
cools. For this reason, we impose a prescribed cooling law

$q^{-}=-\frac{e}{t_{\text{cool}}},$ (10)

where $e$ is the internal energy per unit mass and all details of the cooling
are absorbed by $t_{\text{cool}}$. This parameter defines a typical timescale,
regardless of what process we are taking into account. Often, the ratio
between the cooling time and the dynamical one is chosen to be constant, such
that $\beta=\Omega t_{\text{cool}}=\text{const}$ (Gammie, 2001). In this work,
we use this prescription and hereafter we will refer to the cooling process in
terms of $\beta$.

The stability parameter $Q$ is proportional to the sound speed, thus to the
temperature, so that colder discs are prone to being unstable. In absence of
external heating mechanisms, an initially stable hot disc $(Q>>1)$ will cool
down due to radiative processes, until eventually reaching the marginally
stable state $Q\simeq 1$. At this point, gravitational instability turns on:
the disc develops a spiral structure that, by means of compression and shocks,
leads to an efficient energy dissipation and heating. In this sense, the
$Q-$stability condition acts as a thermostat so that heating turns on only if
the system is sufficiently cold, keeping it in a marginal stable state
$(Q\simeq 1)$. In this regime, spiral perturbations do not grow exponentially,
but their amplitude saturates at some finite value. In order to establish the
conditions under which the disc self regulates, we need to take into account
the heating processes too. As we have said, the generation of spiral density
waves leads to propagation of shocks, because of the supersonic difference of
speed between the spiral pattern and the underneath disc. The self regulation
condition can be obtained by balancing cooling and heating terms (Cossins et
al., 2009), and it gives

$\frac{\delta\Sigma}{\Sigma}=\chi\beta^{-1/2},$ (11)

where $\chi$ is the proportionality factor, that is of the order of unity, and
$\beta$ is the adimensional cooling factor. Hence, a more efficient cooling
(lower $\beta$) gives rise to stronger density perturbations. Note that the
amplitude of the density perturbations is non-linear: indeed, what we did here
is to take into account non-linearities to relate the expected amplitude to
the cooling rate.

We should remember that equation (11) has been obtained under the hypothesis
that the only heating process in the system is the propagation of shocks,
induced by spiral density waves. If we consider irradiated disks, i.e. disks
that are also heated by the central star, the self regulation condition is
different. In particular, irradiation is known to reduce the amplitude of SG
perturbations: Rice et al. (2011) studied the problem of fragmentation in
irradiated disks, and they found that the fragmentation threshold could
decline by approximately a factor of two. Thus, by neglecting this effect, we
expect that in an actual system for given $\delta\Sigma/\Sigma$, the $\beta$
factor should be overestimated.

##  3 Velocity perturbations

In this section, we start from the fluid equations for an infinitesimally thin
self-gravitating accretion disc, we perturb them with a spiral-like
disturbance and we extract the velocity perturbations, following the formalism
of Binney & Tremaine (1987). We use cylindrical coordinates $(r,\phi,z)$ and
we restrict to the $z=0$ plane. The dynamics is characterised by the
continuity equation (equation 12), the two components of Euler’s equation
(equations 13), the Poisson’s equation (equations 14) and the equation of
state of the fluid (equation 15)

$\partial_{t}\Sigma+\frac{1}{r}\partial_{r}(\Sigma
u_{r}r)+\frac{1}{r}\partial_{\phi}(\Sigma u_{\phi})=0,$ (12)
$\displaystyle\partial_{t}u_{r}+u_{r}\partial_{r}u_{r}+\frac{u_{\phi}}{r}\partial_{\phi}u_{r}-\frac{u_{\phi}^{2}}{r}=-\partial_{r}(\Phi+h),$
(13a)
$\displaystyle\partial_{t}u_{\phi}+u_{r}\partial_{r}u_{\phi}+\frac{u_{\phi}}{r}\partial_{\phi}u_{\phi}+\frac{u_{\phi}u_{r}}{r}=-\frac{1}{r}\partial_{\phi}(\Phi+h),$
(13b) $\nabla^{2}\Phi=4\pi G\Sigma\delta(z),$ (14)
$dh=c_{s}^{2}\frac{d\Sigma}{\Sigma},$ (15)

where $h$ is the enthalpy of the fluid, $u_{r}$ and $u_{\phi}$ the radial and
azimuthal component of the velocity field.

For linear analysis, we assume that the spiral perturbation is small compared
to the disc background, and hence can be Fourier-decomposed in time $t$ and
azimuthal angle $\phi$. All the variables $(\Sigma,u_{\phi},u_{r},h,\Phi)$ can
be written as $X=X_{0}(r)+X_{1}(r,\phi,t)$, where $X_{0}$ refers to the basic
state and $X_{1}$ to the perturbation. The basic state of the system is given
by $\Sigma_{0}$, $u_{\phi 0}=r\Omega$, where $\Omega$ is given by equation
(1), $u_{r0}=0$, $h_{0}$.

Keeping only terms that are first order in $X_{1}$, we get the well known
equations for the velocity perturbations

$\begin{split}&u_{r1}=\frac{i}{\Delta}\left[(\omega-m\Omega)\partial_{r}(\Phi_{1}+h_{1})-\frac{2m\Omega}{r}(\Phi_{1}+h_{1})\right],\\\
&u_{\phi
1}=-\frac{1}{\Delta}\left[2B\partial_{r}(\Phi_{1}+h_{1})+\frac{m(\omega-m\Omega)}{r}(\Phi_{1}+h_{1})\right],\end{split}$
(16)

where $B(r)=-\frac{1}{2}\frac{\text{d}(\Omega r)}{\text{d}r}+\Omega$ is one of
the Oort parameter (Oort, 1927), $\Phi_{1}$ is given by (9),
$h_{1}=c_{s}^{2}{\Sigma_{1}}/{\Sigma_{0}}$ and
$\Delta=\kappa^{2}-(\omega-m\Omega)^{2}$.

Now we make some assumptions: firstly, we write the perturbation as

$X_{1}=\text{Re}\left[\delta X(r)e^{i(m\phi-\omega t+\psi)}\right],$ (17)

where $\delta X$ is exclusively a function of the radius. Secondly, we
consider a marginally stable accretion disc with $Q=1$, meaning that
$\Delta=\kappa^{2}$ and $k=k_{\text{uns}}$. Thirdly, both the potential
($\delta\Phi$) and the enthalpy ($\delta h$) perturbations are linked to the
density one that, for a self regulate state, is connected to the basic state
through the cooling rate (11). Hence, we have found a way to express the
perturbations as a function of the basic quantities and the cooling $\beta$:

$\delta\Phi=-\frac{2\pi G}{|k|}\delta\Sigma=-2c_{s}^{2}\chi\beta^{-1/2},$ (18)
$\delta
h=c_{s}^{2}\frac{\delta\Sigma}{\Sigma_{0}}=c_{s}^{2}\chi\beta^{-1/2}=-\frac{1}{2}\delta\Phi.$
(19)

Finally, the velocity perturbations become

$\begin{split}&\delta
u_{r}=\frac{2im\Omega\chi}{r\kappa^{2}}\beta^{-1/2}c_{s}^{2},\\\ &\delta
u_{\phi}=\frac{2iB\chi}{\kappa^{2}}\frac{\text{d}\psi}{\text{d}r}\beta^{-1/2}c_{s}^{2}.\end{split}$
(20)

###  3.1 Nearly Keplerian disc

In this paragraph, we write the velocity perturbations for a nearly Keplerian
disc. This regime is identified by the conditions that

$\frac{\kappa-\Omega_{k}}{\Omega_{k}}<1,\quad\frac{\Omega-\Omega_{k}}{\Omega_{k}}<1,$
(21)

where with the subscript $k$ we identify the Keplerian quantities. This
assumption allows us to write
$\kappa\simeq\Omega\simeq\Omega_{k}\propto{r}^{-3/2}$ and $B\simeq-\Omega/4$.
With these assumptions, the equations for the VPs are

$\begin{split}&\delta
u_{r}=2im\chi\beta^{-1/2}\left(\frac{M_{d}(r)}{M_{\star}}\right)^{2}u_{k}\\\
&\delta
u_{\phi}=-\frac{i\chi\beta^{-1/2}}{2}\left(\frac{M_{d}(r)}{M_{\star}}\right)u_{k},\end{split}$
(22)

where $u_{k}$ is the Keplerian speed, and then the velocity field is given by
$u_{r}=\text{Re}\left[\delta u_{r}e^{i(m\phi-\omega t+\psi)}\right]$ and
$u_{\phi}=r\Omega+\text{Re}\left[\delta u_{\phi}e^{i(m\phi-\omega
t+\psi)}\right]$, where $\Omega$ is given by equation (1). To obtain these
expressions we used that $Q=1$ and that
$\text{d}\psi/\text{d}r=k=k_{\text{uns}}$. Note that $\delta u_{r}/\delta
u_{\phi}\simeq 4mM_{d}/M_{\star}$ : for example, when $m=2$, a relatively
light disk having $M_{d}=0.125M_{\star}$ has $\delta u_{r}=\delta u_{\phi}$.

In the analysis above, we have neglected the effect of pressure gradients.
This is for two main reasons: firstly, it influences only the basic state of
the system, not the perturbations, at leat to first order. Secondly, when we
consider a marginally unstable self-gravitating disk ($Q=1$), we expect the
contribution of the pressure gradient to the velocity field to be sub-dominant
with respect to the self-gravity one. Indeed, for such a disk, the self-
gravitating contribution is of the order of $H/r$, while the pressure term is
$O(H^{2}/r^{2})$ (Kratter & Lodato, 2016; Veronesi et al., 2021). The effects
of pressure gradients are stronger when considering much lower disk/star mass
ratios (e.g., see Rosenfeld et al. 2013). For the massive disks that we
consider in this work, the pressure gradient can thus be neglected. In the
light of this, while considering the pressure gradient is critical when one
wants to explore the basic state, as done in Veronesi et al. (2021), this is
not strictly necessary in our perturbation theory.

###  3.2 Not constant cooling factor

So far we have considered only the case of constant $\beta-$cooling. In
principle, however, self consistent models of GI discs (Clarke, 2009; Rice &
Armitage, 2009) show that $\beta$ varies with the radius (Hall et al., 2016).
This happens because these models give a realistic cooling prescription, i.e.
radiative cooling, and its rate depends on the temperature of the disc at the
mid-plane and on the Rosseland opacity (Bell & Lin, 1994). If one sets the
density profile to be a power law with radius, the cooling prescription can be
written as a collection of power laws with indices $n_{i}$, depending on the
density and the temperature.

In general, we can chose any cooling law $\beta(r,\rho(r),T(r))$ and then
obtain the VPs thorugh equations (20).

##  4 Connection with observations: moment one and channel maps

In the previous section, we have computed the velocity perturbations due to
the presence of a spiral density wave. Now, we want to connect what we have
found to observations: what are the observational imprints of these
perturbations? The observed velocity field of the gas is obtained by
calculating the intensity weighted average velocity of the emission line
profile, i.e. the “moment-1” map. In this work, we do not take into account
line emission processes, and instead we simply compute the projected velocity
field along the line of sight and we study the impact of velocity
perturbations that we have just obtained.

We assume two dimensional polar system of coordinates $(r,\phi)$ centred upon
the the star, so that the velocity vector can be written as
$\mathbf{u}=(u_{r},u_{\phi})$. We consider the disc inclined with an angle
$\theta$, that in the following we take equal to $\pi/6$. Within this
framework, the projected velocity field can be written as

$v_{\text{obs}}=u_{r}\sin\phi\sin\theta+u_{\phi}\cos\phi\sin\theta+v_{\text{syst}},$
(23)

where $v_{\text{syst}}$ is the velocity of the system projected towards the
line of sight. A channel map is defined as the isovelocity contours for a
chosen observed velocity. For a purely Keplerian disc, $u_{r}=0$ and
$u_{\phi}=u_{k}=\sqrt{GM_{\star}/r}$, and we obtain the well known “butterfly
pattern”.

Hence, once we know the radial and the azimuthal components of the perturbed
velocity field (equations 22), we can sketch the projected velocity field
(moment-1 equivalent) and the channel maps, as shown in panels (a), (b) of
figure 1. As already noted in H20, the VPs due to gravitational instability
appear throughout the whole extent of the disc, rather than being localised in
position and velocity, as occurs for the kink produced by an embedded
protoplanet. This is clearly shown in panel (c) of figure 1, where we subtract
the Keplerian field to the perturbed one: an “interlocking fingers” structure
is present, as already pointed out in H20. If we look at the central channel,
the deviations from the Keplerian behaviour exactly match with the fingers
pattern.

We want to underline that in this work we are only considering the projection
of the velocity field along the line of sight, without making any assumptions
about the gas emission processes. In order to convert velocities to fluxes, it
is necessary to include the physics of the gas, specifying the selected tracer
and the emission lines observed and considering also the effect of the beam
size and the eventual presence of observational noise. To do so, we should use
radiative transfer codes, as done in H20.

Figure 1: Moment one map (a) and channel maps (b) for a self gravitating
accretion disc seen with an inclination angle of $\pi/6$ and with a systemic
velocity $v_{\text{syst}}=0$. (c) Left panel: projected map of the velocity
perturbation, after subtraction of the Keplerian field. A system of
interlocking fingers is clearly visible, as already noted by H20. Central
panel: the $v_{\text{obs}}=0$ contour (blue line) overlaid with the spiral
shape (grey line). The deviations from the Keplerian channel (that is simply a
straight line) perfectly match with the spiral pattern. Right panel: surface
density of the disc. The parameters of the disc are the following:
$r_{\text{in}}={1}\text{au}$, $r_{\text{out}}={100}\text{au}$,
$M_{\star}={1}\text{M}_{\odot}$, $M_{d}=0.3\text{M}_{\star}$, $p=-1$,
$\beta=5$, $\alpha_{p}=15^{\circ}$ and $m=2$.

##  5 Discussion and Conclusions

###  5.1 Shape of the wiggle

So far, we have seen that deviations from the Keplerian behaviour are visible
in every channel map, however we now focus on the central velocity channel
that, in a case where $v_{\text{syst}}=0$, corresponds to $v_{\text{obs}}=0$.
As far as the Keplerian case is concerned, the central channel is simply a
straight line, because the radial velocity is zero. In the spiral-perturbed
case, the central channel shows oscillation around the Keplerian behaviour
(i.e. the wiggle): this happens because the spiral wave perturbs both the
azimuthal and the radial component. The amplitude and the radial frequency of
the wiggle depend on the strength of the spiral wave (that is related to the
cooling factor $\beta$ ), the opening angle $\alpha_{p}$ and the number of
arms $m$ (figure 3), and on the structure of the underlying disc, in
particular its mass.

It is possible to quantify the amplitude of the VP considering the integrated
geometrical distance between the perturbed and the unperturbed channel.
Mathematically, a channel map $C_{v}$ is a 2D curve defined parametrically
from a one-dimensional interval $\mathcal{I}$ to a two-dimensional space
$\mathcal{R}^{2}$. In our case, the two-dimensional space is the cylindrical
space $(r,\phi)$ and the interval $\mathcal{I}$ depends on the
parameterization we choose: for simplicity, we take $\mathcal{I}=[0,1]$. For a
given channel velocity $v$, we call the Keplerian channel map $C_{v}^{k}$ and
the perturbed one $C_{v}^{p}$: mathematically speaking, the two channel maps
can be written parametrically as

$C_{v}^{k}(s)=\left(\begin{smallmatrix}f_{r}^{k}(s)\\\ \\\ f_{\phi}^{k}(s)\\\
\end{smallmatrix}\right),\quad
C_{v}^{p}(s)=\left(\begin{smallmatrix}f_{r}^{p}(s)\\\ \\\ f_{\phi}^{p}(s)\\\
\end{smallmatrix}\right),$ (24)

where $s$ is a parameterer in the interval $\mathcal{I}$, in our case
$s\in[0,1]$. The amplitude of the perturbation is then computed as

$\mathcal{A}_{v}=\left[\int_{0}^{1}\text{d}s||C_{v}^{p}-C_{v}^{k}||^{2}\right]^{1/2},$
(25)

that is the length of the curve $C_{v}^{p}-C_{v}^{k}$. Panel (a) of figure 2
schematically shows the quantities involved in equation (25). The amplitude is
determined by both the cooling factor and the disc mass: a smaller $\beta$
generates a bigger deviation from the Keplerian case, because the amplitude of
the density perturbations is inversely proportional to $\beta$ (equation 11).
On panel (b) of figure 2 we show the amplitude of the wiggle
($v_{\text{obs}}=0$) as a function of the cooling parameter $\beta$: it is
possible to describe the relation between $\mathcal{A}$ and $\beta$ as a power
law, with an index of $-1/2$, and this can be easily seen from equations (22).
The role of the cooling factor acts as the planet mass in the case of
planetary kinks: indeed, the amplitude of the kink is determined by the mass
of the embedded protoplanet, and it follows the relation $\mathcal{A}\propto
M_{p}^{1/2}$ (Bollati et al., 2021). In addition, the amplitude of the wiggle
is also determined by the disc mass: in particular, it affects the perturbed
velocities because it is related to the sound speed. Indeed, with the
hypothesis of $Q\simeq 1$ we get $c_{s}=u_{k}{H}/{r}\simeq
u_{k}{M_{d}}/{M_{\star}}$. The disc mass is directly proportional to the sound
speed, and then bigger $c_{s}$ means faster propagation of density waves.
Panel (c) of figure 2 shows how the amplitude of the wiggle depends on the
disc mass. The trend is easily explained by looking at the equations of the
VPs (22): in the radial perturbation, the disc mass appears with a quadratic
dependence, while in the azimuthal one, it appears linearly. The amplitude of
the channel (25) is proportional to the root sum of squares, thus
$\mathcal{A}\propto\sqrt{c_{1}M_{d}^{2}+c_{2}M_{d}^{4}}$, where $c_{1},c_{2}$
are two constants that depend on the other parameters, as the cooling.

Thanks to the analytical expression for the VPs, we can infer in what regimes
the radial or azimuthal perturbation dominates the wiggle. As we have already
noted, for larger $M_{d}/M_{\star}$ the radial perturbation dominates. This is
a crucial point because a perturbation of the channel
$v_{\text{obs}}=v_{\text{syst}}$ is visible only if there is a radial
perturbation: indeed, as in the Keplerian case, if the velocity is purely
azimuthal, the central channel is a straight line.

Figure 2: Schematic view of how the amplitude of the wiggle is computed (a):
the red line is the Keplerian channel map, while the blue line is the
perturbed one. The amplitude of the perturbation is computed using equation
25. Amplitude of the central channel as a function (b) of the cooling
parameter $\beta$ for a disk mass $M_{d}/M_{\star}=0.3$ and (c) of the mass of
the disc for a cooling parameter $\beta=5$. The black dots are the results of
the analytical model so far described.

Unfortunately, the mass of the disc and the cooling $\beta$ are degenerate
parameters when we consider the shape of the wiggle: as a matter of fact, they
both contribute to its amplitude. Is it possible to break the degeneracy? This
can be done if an independent method to measure $M_{d}$ is available. For
example, using again the gas kinematics, one could measure deviations from
Keplerianity in the azimuthally averaged rotation curve, and compute the mass
of the disc through equation (2), breaking the degeneracy. Veronesi et al.
(2021) did it for the system Elias 2-27, giving a dynamical estimate of the
mass of the disc $M_{d}\simeq 0.08M_{\odot}\simeq 0.17M_{\star}$. Conversely,
when an approach like this is not possible, we can give a rough estimate of
the disc mass through dust emission. Indeed, from dust thermal emission it is
possible to measure the dust mass of the disc and then, assuming a dust to gas
ratio, we can estimate also the total mass of the disc.

Interestingly, if we look at the perturbed velocity field, there could also be
a purely kinematical way to break the degeneracy between disk mass and
cooling. Indeed if we write the observed velocity field for $\phi=\pi/2$
(semi-minor axis of the disk) we get

$v_{\text{obs}}=u_{r}\sin\theta=\text{Re}\left[\delta
u_{r}(r)e^{i(m\pi/2+\psi)}\right]\sin\theta=f_{1}(r)\beta^{-1/2}\left(\frac{M_{d}}{M_{\star}}\right)^{2}u_{k},$
(26)

where $f_{1}$ is a known function of radius; conversely, for $\phi=0$ (semi-
major axis of the disk) we get

$v_{\text{obs}}=u_{\phi}\sin\theta=\left(r\Omega+\text{Re}\left[\delta
u_{\phi}(r)e^{i\psi}\right]\right)\sin\theta=f_{2}(r)\beta^{-1/2}\frac{M_{d}}{M_{\star}}u_{k}+r\Omega\sin\theta,$
(27)

where $f_{2}$ is a known function of the radius. Since the perturbed
velocities scale differently with disk mass, if we could measure accurately
the ratio of the two components of the perturbed velocity, we could in
principle directly obtain a measurement of the disk mass. However, we note
that it is challenging to extract these information from an actual
observation.

Breaking the degeneracy allows us to constrain the cooling parameter $\beta$,
which gives important information about the physical processes that are
happening in the protoplanetary environment, and on the tendency of the disc
to fragment into bound clumps.

So far we have described what determines the amplitude of the wiggle: as far
as its frequency is concerned, it is determined by the pitch angle and by the
number of spiral arms. In figure 3 we show the shape of the wiggle for
different values of $\alpha_{p}$ and $m$. We clearly see that the frequency of
the wiggle is bigger when decreasing $\alpha_{p}$ and increasing $m$.

Figure 3: Shape of the wiggle varying (from left to right) the cooling factor
$\beta$, the mass of the disc $M_{d}$, the pitch angle $\alpha_{p}$ in degrees
and the azimuthal wavenumber $m$. The reference disc parameters are
$r_{\text{in}}={1}\text{au}$, $r_{\text{out}}={50}\text{au}$,
$M_{\star}={1}{\text{M}_{\odot}}$, $M_{d}=0.3\text{M}_{\star}$, $p=-1$,
$\beta=5$, $\alpha_{p}=15^{\circ}$ and $m=2$.

###  5.2 Limitations

As can be seen in figure 3, the number of spiral arms slightly influences the
amplitude of the wiggle: all the calculations have been made under the WKB
assumption, that requires $m/rk<<1$. Thus, for high $m$, this assumption is
not valid anymore: in fact, we are not considering the relations between the
number of spiral arms and the mass of the disc, or its thickness; the only way
to take into account these non-linear effects is by means of numerical
simulations (Terry et al., 2021). This argument is better understood when
looking at equation 3: in the tight winding approximation (WKB), $m$ does not
enter explicitly in the dispersion relation (except in the Doppler-shifted
frequency), thus both axisymmetric ($m=0$) and non-axisymmetric ($m\neq 0$)
perturbations have the same instability threshold. However, it is well known
(Ostriker & Peebles, 1973) that massive discs are subject to large scale non-
axisymmetric instabilities even though $Q>1$. Indeed, a local dispersion
relation can also be obtained in the case of open spiral structures (Lau &
Bertin, 1978): it is a cubic rather than quadratic expression in $k$ and
depends on a new dimensionless parameter

$\mathcal{J}=m\frac{\pi
G\Sigma}{r\kappa^{2}}\frac{4\Omega}{\kappa}\left|\frac{\mathrm{d}\ln\Omega}{\mathrm{d}\ln
r}\right|^{1/2}\approx\sqrt{6}m\frac{M_{\mathrm{d}}}{M_{*}}.$ (28)

Being a function of $m$ and $M_{d}$, the $\mathcal{J}$ criterion takes into
account how massive the disc is and the number of spiral arms: massive discs
are prone to exhibit low-$m$ modes instability. Indeed, there is a link
between the number of spiral arms and the mass of the disc (Dipierro et al.,
2014a).

Another important point to stress is that we constructed a 2D model of the
disc, neglecting its height. Thus, we are basically considering only what
happens in the disc mid-plane. The main effect of the disc thickness is to
“dilute” the gravity field, and this can be incorporated into the quadratic
dispersion relation

$D_{\text{thick}}(\omega,k,m)=(\omega-m\Omega)^{2}-c_{s}^{2}k^{2}-\kappa^{2}+2\pi
G\Sigma\frac{|k|}{1+|k|H}.$ (29)

In this case, the stability criterion only changes slightly, becoming
$Q\gtrsim 0.647$ (Kratter & Lodato, 2016). Observationally speaking, the disc
thickness is important when we take into account the molecular line emission
of CO isotopologues. As a matter of fact, 12CO, that is the most abundant
isotopologue, is not a good tracer of the disc mid-plane, because it becomes
optically thick at the disc surface. On the contrary, other less abundant CO
isotopologues as 13CO or C18O have more optically thin lines and as a
consequence they trace better the disc mid-plane (Miotello et al., 2014).

In addition, our analysis takes into account a single spiral mode $m$, while
it is well known that for small disk-to-star mass ratio, there could be a
superposition of modes. However this is not an actual problem: indeed, we know
that after filtering out through the ALMA response (that we do not do in this
paper), as shown in Dipierro et al. (2014b), only the dominant mode appears.
This makes our single-mode analysis still valid.

###  5.3 Constraining the cooling factor - mock test

So far, we have seen that the cooling factor deeply influences the shape of
the channel maps. Thanks to this property, we propose a method with which we
can constrain effectively the cooling factor of observed systems. In order to
verify the accuracy of the calculations we have made, we apply the method just
described to a numerical simulation. We perform an SPH simulation using the
code PHANTOM (Price et al., 2018b). This code is widely used in astrophysical
community to study gas and dust dynamics in accretion discs (Toci et al.,
2020a; Ragusa et al., 2020; Veronesi et al., 2020); in this work, we used the
so-called “one fluid” method and we simulated a gas only disc, neglecting the
dust component. The initial conditions of the disc are
$r_{\text{in}}={1}\text{au}$, $r_{\text{out}}={50}\text{au}$, $\Sigma\propto
r^{-1}$, $M_{d}=0.5\text{M}_{\star}$, $M_{\star}={1}{\text{M}_{\odot}}$. The
cooling factor $\beta$ has been set to $\beta=8$ and the two parameters that
control the artificial viscosity to $\alpha^{AV}=0.1$, $\beta^{AV}=0.2$, in
order to reduce as much as possible the effects of artificial dissipation
(Lodato & Rice, 2004). The initial sound speed profile follows a simple power
law $c_{s}\propto r^{-0.25}$. However, this profile is rapidly modified by the
cooling. For this simulation, we used $N=5\times 10^{5}$ particles of gas. the
simulation shows a predominance of the $m=2$ azimuthal mode and the computed
pitch angle is $\alpha_{p}\sim 13^{\circ}$.

We now constrain the cooling factor using the method described previously, and
make a comparison with the actual value set in the simulation. First of all we
compute the rotation curve, azimuthally averaging $u_{\phi}(r,\phi)$. Then, we
find the value of the disc mass that best describe the curve using equation
(2): the best value corresponds to $M_{d}=0.5\text{M}_{\star}$, as expected.
We have broken the degeneracy between the mass and the cooling: thus, we are
now able to constrain the cooling $\beta$ through the wiggle amplitude. Figure
4 shows the central channel of the projected velocity field of the simulation
compared with the one from the analytical model, for different spiral angles:
as expected, a wiggle is present. The amplitude of the simulated
wiggle111Here, we refer to panel (a) of Figure 4. is
$\mathcal{A}_{\text{sim}}={0.11}\text{rad}$. The amplitude of the wiggle as a
function of the cooling factor, for the parameters previously reported, is
described by $\mathcal{A}(\beta)=\mathcal{A}_{\text{in}}\beta^{-0.5}$, where
$\mathcal{A}_{\text{in}}={0.33}\text{rad}$. The estimated cooling factor is
then
$\beta_{\text{sim}}=(\mathcal{A}_{\text{sim}}/\mathcal{A}_{\text{in}})^{-2}\simeq
9$, that is in good agreement with the real value $\beta=8$. The overestimated
value of $\beta$ can actually be interpreted by the lack of viscosity in our
analytical model. Indeed, because of its dissipative nature, we expect it to
damp GI-driven perturbations, resulting in a less pronounced wiggle. This
behaviour is visible in the comparison with the numerical simulations, hence
the simulated perturbation appears less wide than the analytical one.

Figure 4: Comparison between analyical (black line) and numerical (green line)
perturbation for different viewing angle of the simulated GI disk. Each
picture is rotated along the $z$-axis of an angle of $\pi/2$.

Figure 4 shows a comparison between the analytical and the numerical wiggle,
for different viewing angle of the spiral structure (the observer is assumed
to be along a vertical line on the bottom of the images). There is a very good
agreement in panel (a) and (c) for which the line of sight intercept the
largest extent of the spiral, while the perturbation is overestimated in panel
(b) and (d) in which the prominent spiral structure lies on a line
perpendicular to the line of sight. For the former case, the agreement between
our model and the simulation is remarkable. For the two other orientations,
while there is good agreement in the inner disc (where the perturbation is
however smaller), our model overestimates the perturbation in the outer disc,
where in the simulation the spiral structure vanishes. This suggests that in
actual observations, our analysis is most reliable when a density spiral (for
example traced by the dust continuum) is also visible superimposed to the
kinematical wiggle.

###  5.4 Conclusions

In this work we have analytically studied the velocity perturbations in a
self-gravitating disc caused by the presence of a spiral density wave in the
WKB regime. We then applied this result to obtain the projected velocity field
(moment one equivalent) and the channel maps, studying their deviations from
the Keplerian case. We found what H20 have already seen from numerical
simulations, that deviations from Keplerian rotation are a global phenomenon,
resulting in velocity “kinks” across the entire radial and azimuthal extent of
the disc. The kinematics deviations, called GI wiggles, depend on the
structure of the spiral density wave, namely its amplitude (connected to the
cooling and to the disc mass) and its radial frequency (connected to the pitch
angle and to the azimuthal wavenumber).

Pinte et al. (2020) found nine deviations from Keplerian rotation pattern in
the DSHARP circumstellar discs: three of them, Elias 2-27 (Pérez et al., 2016;
Paneque-Carreño et al., 2021), IM Lup and WaOph6 show also spiral structures
in the millimetric continuum emission. These three systems are believed to be
self-gravitating (Huang et al., 2018), thus their deviations from Keplerian
rotation may be interpreted as wiggles. In addition, Veronesi et al. (2021)
have shown that the rotation curve of Elias 2-27 is better described adding
the contribution of the disc gravitational potential, meaning that the effects
of disc self gravity are not negligible.

Higher resolution observations of systems like those will make it possible to
investigate the cooling of protoplanetary discs: indeed, the degeneracy
between mass and cooling can be broken by means of the rotation curve, and
thus the cooling parameter $\beta$ can be constrained effectively through the
wiggle’s amplitude, as we have shown in section 5.3. Knowing more about the
cooling will give us insights about the gravitational instability process.

## Acknowledgements

The authors thank the anonymous referee for insightful comments and
suggestions, and Richard Booth, Cathie Clarke and Pietro Curone for useful
discussions. This project has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No 823823 (DUSTBUSTERS RISE project). BV acknowledges
funding from the ERC CoG project PODCAST No 864965.

## References

  * Bell & Lin (1994) Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987, doi: 10.1086/174206
  * Bertin & Lin (1996) Bertin, G., & Lin, C. C. 1996, Spiral structure in galaxies a density wave theory (MIT Press)
  * Bertin & Lodato (1999) Bertin, G., & Lodato, G. 1999, A&A, 350, 694.  https://arxiv.org/abs/astro-ph/9908095
  * Binney & Tremaine (1987) Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton Press)
  * Bollati et al. (2021) Bollati, F., Lodato, G., Price, D. J., & Pinte, C. 2021, arXiv e-prints, arXiv:2103.00016.  https://arxiv.org/abs/2103.00016
  * Cadman et al. (2020) Cadman, J., Rice, K., Hall, C., Haworth, T. J., & Biller, B. 2020, MNRAS, 492, 5041, doi: 10.1093/mnras/staa187
  * Clarke (2009) Clarke, C. J. 2009, MNRAS, 396, 1066, doi: 10.1111/j.1365-2966.2009.14774.x
  * Cossins et al. (2009) Cossins, P., Lodato, G., & Clarke, C. J. 2009, MNRAS, 393, 1157, doi: 10.1111/j.1365-2966.2008.14275.x
  * Cuello et al. (2019) Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114, doi: 10.1093/mnras/sty3325
  * Cuello et al. (2020) Cuello, N., Louvet, F., Mentiplay, D., et al. 2020, MNRAS, 491, 504, doi: 10.1093/mnras/stz2938
  * Dipierro et al. (2014a) Dipierro, G., Lodato, G., Testi, L., & de Gregorio Monsalvo, I. 2014a, MNRAS, 444, 1919, doi: 10.1093/mnras/stu1584
  * Dipierro et al. (2014b) —. 2014b, MNRAS, 444, 1919, doi: 10.1093/mnras/stu1584
  * Dipierro et al. (2015) Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73, doi: 10.1093/mnrasl/slv105
  * Dipierro et al. (2018) Dipierro, G., Ricci, L., Pérez, L., et al. 2018, MNRAS, 475, 5296, doi: 10.1093/mnras/sty181
  * Dong & Fung (2017) Dong, R., & Fung, J. 2017, ApJ, 835, 38, doi: 10.3847/1538-4357/835/1/38
  * Dong et al. (2015) Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJ, 809, L5, doi: 10.1088/2041-8205/809/1/L5
  * Gammie (2001) Gammie, C. F. 2001, ApJ, 553, 174, doi: 10.1086/320631
  * Hall et al. (2017) Hall, C., Forgan, D., & Rice, K. 2017, MNRAS, 470, 2517, doi: 10.1093/mnras/stx1244
  * Hall et al. (2016) Hall, C., Forgan, D., Rice, K., et al. 2016, MNRAS, 458, 306, doi: 10.1093/mnras/stw296
  * Hall et al. (2020) Hall, C., Dong, R., Teague, R., et al. 2020, ApJ, 904, 148, doi: 10.3847/1538-4357/abac17
  * Hirose & Shi (2017) Hirose, S., & Shi, J.-M. 2017, MNRAS, 469, 561, doi: 10.1093/mnras/stx824
  * Huang et al. (2018) Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018, ApJ, 869, L43, doi: 10.3847/2041-8213/aaf7a0
  * Johnson & Gammie (2003) Johnson, B. M., & Gammie, C. F. 2003, ApJ, 597, 131, doi: 10.1086/378392
  * Kratter & Lodato (2016) Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271, doi: 10.1146/annurev-astro-081915-023307
  * Kratter & Matzner (2006) Kratter, K. M., & Matzner, C. D. 2006, MNRAS, 373, 1563, doi: 10.1111/j.1365-2966.2006.11103.x
  * Kratter et al. (2010) Kratter, K. M., Murray-Clay, R. A., & Youdin, A. N. 2010, ApJ, 710, 1375, doi: 10.1088/0004-637X/710/2/1375
  * Lau & Bertin (1978) Lau, Y. Y., & Bertin, G. 1978, ApJ, 226, 508, doi: 10.1086/156635
  * Lin & Shu (1964) Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646, doi: 10.1086/147955
  * Lodato & Rice (2004) Lodato, G., & Rice, W. K. M. 2004, MNRAS, 351, 630, doi: 10.1111/j.1365-2966.2004.07811.x
  * Lodato & Rice (2005) —. 2005, MNRAS, 358, 1489, doi: 10.1111/j.1365-2966.2005.08875.x
  * Miotello et al. (2014) Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96, doi: 10.1051/0004-6361/201424712
  * Oort (1927) Oort, J. H. 1927, Bull. Astron. Inst. Netherlands, 3, 275
  * Ostriker & Peebles (1973) Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467, doi: 10.1086/152513
  * Paneque-Carreño et al. (2021) Paneque-Carreño, T., Pérez, L. M., Benisty, M., et al. 2021, ApJ, 914, 88, doi: 10.3847/1538-4357/abf243
  * Pérez et al. (2016) Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519, doi: 10.1126/science.aaf8296
  * Pinte et al. (2018) Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13, doi: 10.3847/2041-8213/aac6dc
  * Pinte et al. (2020) —. 2020, ApJ, 890, L9, doi: 10.3847/2041-8213/ab6dda
  * Price (2007) Price, D. J. 2007, PASA, 24, 159, doi: 10.1071/AS07022
  * Price et al. (2018a) Price, D. J., Cuello, N., Pinte, C., et al. 2018a, MNRAS, 477, 1270, doi: 10.1093/mnras/sty647
  * Price et al. (2018b) Price, D. J., Wurster, J., Tricco, T. S., et al. 2018b, PASA, 35, e031, doi: 10.1017/pasa.2018.25
  * Ragusa et al. (2020) Ragusa, E., Alexander, R., Calcino, J., Hirsh, K., & Price, D. J. 2020, MNRAS, 499, 3362, doi: 10.1093/mnras/staa2954
  * Rice & Armitage (2009) Rice, W. K. M., & Armitage, P. J. 2009, MNRAS, 396, 2228, doi: 10.1111/j.1365-2966.2009.14879.x
  * Rice et al. (2011) Rice, W. K. M., Armitage, P. J., Mamatsashvili, G. R., Lodato, G., & Clarke, C. J. 2011, MNRAS, 418, 1356, doi: 10.1111/j.1365-2966.2011.19586.x
  * Rice et al. (2005) Rice, W. K. M., Lodato, G., & Armitage, P. J. 2005, MNRAS, 364, L56, doi: 10.1111/j.1745-3933.2005.00105.x
  * Rosenfeld et al. (2013) Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16, doi: 10.1088/0004-637X/774/1/16
  * Rosotti et al. (2020) Rosotti, G. P., Benisty, M., Juhász, A., et al. 2020, MNRAS, 491, 1335, doi: 10.1093/mnras/stz3090
  * Stamatellos et al. (2007) Stamatellos, D., Hubber, D. A., & Whitworth, A. P. 2007, MNRAS, 382, L30, doi: 10.1111/j.1745-3933.2007.00383.x
  * Stamatellos & Whitworth (2009) Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 400, 1563, doi: 10.1111/j.1365-2966.2009.15564.x
  * Teague et al. (2018) Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018, ApJ, 860, L12, doi: 10.3847/2041-8213/aac6d7
  * Terry et al. (2021) Terry, J. P., Hall, C., Longarini, C., et al. 2021, arXiv e-prints, arXiv:2108.11282.  https://arxiv.org/abs/2108.11282
  * Toci et al. (2020a) Toci, C., Lodato, G., Christiaens, V., et al. 2020a, MNRAS, 499, 2015, doi: 10.1093/mnras/staa2933
  * Toci et al. (2020b) Toci, C., Lodato, G., Fedele, D., Testi, L., & Pinte, C. 2020b, ApJ, 888, L4, doi: 10.3847/2041-8213/ab5c87
  * Veronesi et al. (2019) Veronesi, B., Lodato, G., Dipierro, G., et al. 2019, MNRAS, 489, 3758, doi: 10.1093/mnras/stz2384
  * Veronesi et al. (2021) Veronesi, B., Paneque-Carreño, T., Lodato, G., et al. 2021, ApJ, 914, L27, doi: 10.3847/2041-8213/abfe6a
  * Veronesi et al. (2020) Veronesi, B., Ragusa, E., Lodato, G., et al. 2020, MNRAS, 495, 1913, doi: 10.1093/mnras/staa1278