# Investigating Protoplanetary Disk Cooling through Kinematics: Analytical GI Wiggle

Cristiano Longarini, Giuseppe Lodato, Claudia Toci, Benedetta Veronesi, Cassandra Hall, Ruobing Dong, Jason Patrick Terry

► To cite this version:

Cristiano Longarini, Giuseppe Lodato, Claudia Toci, Benedetta Veronesi, Cassandra Hall, et al. Investigating Protoplanetary Disk Cooling through Kinematics: Analytical GI Wiggle. The Astrophysical Journal Letters, 2021, 920, 10.3847/2041-8213/ac2df6. insu-03711488

HAL Id: insu-03711488

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

Submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Investigating Protoplanetary Disk Cooling through Kinematics: Analytical GI Wiggle

Cristiano Longarini<fnref n="1" />, Giuseppe Lodato<fnref n="1" />, Claudia Toci<fnref n="1" />, Benedetta Veronesi<fnref n="1" />,<fnref n="2" />, Cassandra Hall<fnref n="3" />,<fnref n="4" />, Ruobing Dong<fnref n="5" />, and Jason Patrick Terry<fnref n="3" />,<fnref n="4" />

Received 2021 August 20; revised 2021 September 28; accepted 2021 October 7; published 2021 October 21

## Abstract

It is likely that young protostellar disks undergo a self-gravitating phase. Such systems are characterized by the presence of a spiral pattern that can be either in a quasi-steady state or in a nonlinear 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 (≤1% of the Keplerian speed) thanks to high spatial and spectral resolution observations of protostellar disks. In this work, we investigate the kinematic signatures of gravitational instability: we perform an analytical study of the linear response of a self-gravitating disk to a spiral-like perturbation, focusing our attention on the velocity field perturbations. We show that unstable disks have clear kinematic imprints into the gas component across the entire disk 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 disk cooling.

*Unified Astronomy Thesaurus concepts: Astrophysical fluid dynamics (101); Stellar accretion disks (1579); Gravitational instability (668)*

## 1. Introduction

For cold massive disks, the role of the disk self-gravity becomes dynamically important, affecting the vertical and radial structure of the system (Bertin & Lodato 1999; Kratter & Lodato 2016). In this context, gravitational instabilities (GIs) 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; Binney & Tremaine 1987; Bertin & Lin 1996): as far as protostellar disks are concerned, the results are quantitatively similar.

On the one hand, a possible outcome of GI in protostellar environments is the fragmentation of the disk: 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 disk, by means of a global spiral perturbation (Lodato & Rice 2004, 2005).

High-resolution observations with the Atacama Large Millimeter/submillimeter Array (ALMA) have revealed that most of the observed protostellar disks possess substructures such 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 an 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 an 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 characterize self-gravitating disks, with a typical $m \sim M_*/M_d$, where $M_d$ is the disk mass and $M_*$ 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 disk creates a localized 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 H20) show, based on hydrodynamical simulations, that the spiral perturbation deeply affects the gas kinematics: in particular, it creates a global (rather than a localized) deviation from Keplerian observed velocity—a “global kink”—dubbed GI wiggle by H20 that is apparent in the moment-1 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 disk to a spiral density wave in the WKB regime. We focus our attention on the velocity field perturbations (VPs), and we show that they can be written as a function of disk 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 disks.

The Letter is organized as follows. In Section 2 we summarize 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 the WKB regime. In Section 4 we discuss the observational aspects of the wiggles,

<footnote n="1">¹Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, I-20133 Milano, Italy; cristiano.longarini@unimi.it</footnote>

<footnote n="2">²Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France</footnote>

<footnote n="3">³Department of Physics and Astronomy, The University of Georgia, Athens, GA 30602, USA</footnote>

<footnote n="4">⁴Center for Simulational Physics, The University of Georgia, Athens, GA 30602, USA</footnote>

<footnote n="5">⁵Department of Physics & Astronomy, University of Victoria, 3800 Finnerty Road, Victoria, BC V8P 5C2, Canada</footnote>
projecting the perturbed velocity field along the line of sight. To conclude, in Section 5 we discuss the shape of the perturbations and the limitations of this work, and we present a mock experiment to test our predictions.

## 2. Self-gravitating Gaseous Disks

### 2.1. The Disk Potential

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

$$ \Omega^2 = \frac{GM_*}{r^3} + \frac{1}{r} \frac{\partial \Phi_d}{\partial r}, \quad (1) $$

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

$$ \frac{\partial \Phi_d}{\partial r}(r, z) = \frac{G}{r} \int_0^\infty dr' \left[ K(\zeta) - \frac{1}{4} \left( \frac{\zeta^2}{1 - \zeta^2} \right) \times \left( \frac{r'}{r} - \frac{r}{r'} + \frac{z^2}{rr'} \right) E(\zeta) \right] \sqrt{\frac{r'}{r}} \zeta \Sigma(r'), \quad (2) $$

where $E(\zeta)$ and $K(\zeta)$ are the complete elliptic integrals of the first kind, $\zeta^2 = 4rr'/[(r+r')^2+z^2]$, and $\Sigma$ is the disk surface density. Deviations from the Keplerian behavior 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 disk 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_s^2 k^2 + 2\pi G\Sigma|k| - \kappa^2 = 0, \quad (3) $$

where $\omega$ is the wave angular frequency, $k$ is the radial wavenumber, $m$ is the azimuthal wavenumber, $\kappa(r)$ is the epicyclic frequency, $\Omega(r)$ is the angular frequency, and $c_s(r)$ is the sound speed. This relationship has been obtained for an infinitesimally thin disk 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}, \quad (4) $$

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

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

where $H_{\text{sg}} = c_s^2/\pi G \Sigma$ is the disk thickness in the self-gravitating case.

A marginally stable disk is characterized by having $Q \simeq 1$: in this case, the only expected excited modes would be $k \simeq k_{\text{uns}}$. Thus, Equation (3) tells us that $\Omega \simeq \omega/m = \Omega_p$, where $\Omega_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 characterized by having

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

where $\psi(r)$ is the shape function and it holds that $k = d\psi/dr$. It is useful to introduce the radial wavenumber $k$, which is the radial derivative of the shape function; the sign on $k$ determines whether 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 to as the pitch angle; it is given by

$$ \tan \alpha_p = \left| \frac{r \partial \phi}{\partial r} \right|^{-1}, \quad (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}[\delta\Sigma(r)e^{i(m\phi - \omega t + \psi(r))}], \quad (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). \quad (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 nonlinear evolution. To do so, it is necessary to introduce a parameter to capture the radiative processes of the disk.

The nonlinear 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 past several years, a lot of effort was devoted to the modeling of realistic thermodynamics (Johnson & Gammie 2003; Stamatellos & Whitworth 2009; Hirose & Shi 2017).

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

$$ q^{-} = -\frac{e}{t_{\text{cool}}}, \quad (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 disks are prone to being unstable. In the absence of external heating mechanisms, an initially stable hot disk ($Q \gg 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 disk 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 marginally 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 disk 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 disk. 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}, \qquad (11) $$

where $\chi$ is the proportionality factor, which 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 nonlinear: indeed, what we did here is to take into account nonlinearities 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 2. 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 disk, 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 characterized by the continuity equation (Equation (12)), the two components of Euler's equation (Equation (13a)), the Poisson's equation (Equation (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, \qquad (12) $$

$$ \begin{aligned} \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), \end{aligned} \qquad (13a) $$

$$ \begin{aligned} \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), \end{aligned} \qquad (13b) $$

$$ \nabla^2 \Phi = 4\pi G \Sigma \delta(z), \qquad (14) $$

$$ dh = c_s^2 \frac{d\Sigma}{\Sigma}, \qquad (15) $$

where $h$ is the enthalpy of the fluid and $u_r$ and $u_\phi$ are the radial and azimuthal components of the velocity field.

For linear analysis, we assume that the spiral perturbation is small compared to the disk 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{aligned} 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{aligned} \qquad (16) $$

where $B(r) = -\frac{1}{2} \frac{d(\Omega r)}{dr} + \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: first, we write the perturbation as

$$ X_1 = \mathrm{Re} [\delta X(r) e^{i(m\phi - \omega t + \psi)}], \qquad (17) $$

where $\delta X$ is exclusively a function of the radius. Second, we consider a marginally stable accretion disk with $Q=1$, meaning that $\Delta = \kappa^2$ and $k = k_{\mathrm{uns}}$. Third, both the potential ($\delta\Phi$) and the enthalpy ($\delta h$) perturbations are linked to the density one that, for a self-regulated 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}, \qquad (18) $$

$$ \delta h = c_s^2 \frac{\delta \Sigma}{\Sigma_0} = c_s^2 \chi \beta^{-1/2} = -\frac{1}{2} \delta \Phi. \qquad (19) $$

Finally, the velocity perturbations become

$$ \begin{aligned} \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{d\psi}{dr} \beta^{-1/2} c_s^2. \end{aligned} \qquad (20) $$
### 3.1. Nearly Keplerian Disk

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

$$ \frac{\kappa - \Omega_k}{\Omega_k} < 1, \quad \frac{\Omega - \Omega_k}{\Omega_k} < 1, \qquad (21) $$

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

$$ \delta u_r = 2im\chi\beta^{-1/2} \left( \frac{M_d(r)}{M_*} \right)^2 u_k, \qquad (22) $$

$$ \delta u_{\phi} = -\frac{i\chi\beta^{-1/2}}{2} \left( \frac{M_d(r)}{M_*} \right) u_k, $$

where $u_k$ is the Keplerian speed, and then the velocity field is given by $u_r = \text{Re}[\delta u_r e^{i(m\phi - \omega t + \psi)}]$ and $u_{\phi} = r\Omega + \text{Re}[\delta u_{\phi} e^{i(m\phi - \omega t + \psi)}]$, where $\Omega$ is given by Equation (1). To obtain these expressions we used that $Q=1$ and that $d\psi/dr = k = k_{\text{uns}}$. Note that $\delta u_r/\delta u_{\phi} \approx 4mM_d/M_*$: for example, when $m=2$, a relatively light disk having $M_d=0.125 M_*$ 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: first, it influences only the basic state of the system, not the perturbations, at least to first order. second, 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 disks (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 disk at the midplane 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 choose any cooling law $\beta(r, \rho(r), T(r))$ and then obtain the VPs through Equation (20).

## 4. Connection with Observations: Moment-1 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)$ centered upon the star, so that the velocity vector can be written as $\mathbf{u} = (u_r, u_{\phi})$. We consider the disk 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}}, \qquad (23) $$

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

Hence, once we know the radial and the azimuthal components of the perturbed velocity field (Equation (22)), we can sketch the projected velocity field (moment-1 equivalent) and the channel maps, as shown in panels (a) and (b) of Figure 1. As already noted in H20, the VPs due to gravitational instability appear throughout the whole extent of the disk, rather than being localized 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 behavior 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.

## 5. Discussion and Conclusions

### 5.1. Shape of the Wiggle

So far, we have seen that deviations from the Keplerian behavior 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 behavior (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 (which is related to the cooling
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Moment-1 map (a) and channel maps (b) for a self-gravitating accretion disk seen with an inclination angle of π/6 and with a systemic velocity v<sub>syst</sub> = 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<sub>obs</sub> = 0 contour (blue line) overlaid with the spiral shape (gray 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 disk. The parameters of the disk are the following: r<sub>in</sub> = 1 au, r<sub>out</sub> = 100 au, M* = 1 M<sub>☉</sub>, M<sub>d</sub> = 0.3M*, p = −1, β = 5, α<sub>p</sub> = 15°, and m = 2.</figcaption></figure>

factor $\beta$), the opening angle $\alpha_p$, and the number of arms $m$ (Figure 3), and on the structure of the underlying disk, 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
<figure><img src="image_4.png" /><figcaption>**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_*$ = 0.3 and (c) of the mass of the disk for a cooling parameter $\beta$ = 5. The black dots are the results of the analytical model so far described.</figcaption></figure>

speaking, the two channel maps can be written parametrically as

$$C_v^k(s) = \begin{pmatrix} f_r^k(s) \\ f_\phi^k(s) \end{pmatrix}, \quad C_v^p(s) = \begin{pmatrix} f_r^p(s) \\ f_\phi^p(s) \end{pmatrix}, \qquad (24)$$

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

$$A_v = \left[ \int_0^1 ds ||C_v^p - C_v^k||^2 \right]^{1/2}, \qquad (25)$$

which 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 disk 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)). In 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 $A$ and $\beta$ as a power law, with an index of $-1/2$, and this can be easily seen from Equation (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 $A \propto M_p^{1/2}$ (Bollati et al. 2021). In addition, the amplitude of the wiggle is also determined by the disk 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_*$. The disk 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 disk mass. The trend is easily explained by looking at the equation of the VPs (22): in the radial perturbation, the disk 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 $A \propto \sqrt{c_1 M_d^2 + c_2 M_d^4}$, where $c_1$ and $c_2$ are two constants that depend on the other parameters, such 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_*$ 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.

Unfortunately, the mass of the disk 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 disk 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 disk $M_d \simeq 0.08 M_\odot \simeq 0.17 M_*$. Conversely, when an approach like this is not possible, we can give a rough estimate of the disk mass through dust emission. Indeed, from dust thermal emission it is possible to measure the dust mass of the disk and then, assuming a dust-to-gas ratio, we also can estimate the total mass of the disk.

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$ (semiminor axis of the disk) we get

$$\begin{align}
v_{\text{obs}} &= u_r \sin \theta = \text{Re}[\delta u_r(r) e^{i(m\pi/2 + \psi)}] \nonumber \\
\sin \theta &= f_1(r) \beta^{-1/2} \left(\frac{M_d}{M_*}\right)^2 u_k, \tag{SPH26}
\end{align}$$

where $f_1$ is a known function of radius; conversely, for $\phi = 0$ (semimajor axis of the disk) we get

$$\begin{align}
v_{\text{obs}} &= u_{\phi} \sin \theta = (r\Omega + \text{Re}[\delta u_{\phi}(r) e^{i\psi}]) \nonumber \\
\sin \theta &= f_2(r) \beta^{-1/2} \frac{M_d}{M_*} u_k + r\Omega \sin \theta, \tag{27}
\end{align}$$
<figure><img src="image_5.png" /><figcaption>Figure 3. Shape of the wiggle varying (from left to right) the cooling factor $\beta$, the mass of the disk $M_d$, the pitch angle $\alpha_p$ in degrees, and the azimuthal wavenumber $m$. The reference disk parameters are $r_{in} = 1$ au, $r_{out} = 50$ au, $M_* = 1 M_\odot$, $M_d = 0.3 M_*$, $p = -1$, $\beta = 5$, $\alpha_p = 15^\circ$, and $m = 2$.</figcaption></figure>

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 this 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 disk 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$.

## 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 \ll 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 disk, or its thickness; the only way to take into account these nonlinear 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 nonaxisymmetric ($m \neq 0$) perturbations have the same instability threshold. However, it is well known (Ostriker & Peebles 1973) that massive disks are subject to large-scale nonaxisymmetric 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{d \ln \Omega}{d \ln r} \right|^{1/2} \approx \sqrt{6} m \frac{M_d}{M_*} \quad (28) $$

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

Another important point to stress is that we constructed a 2D model of the disk, neglecting its height. Thus, we are basically

considering only what happens in the disk midplane. The main effect of the disk 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} \quad (29) $$

In this case, the stability criterion only changes slightly, becoming $Q \gtrsim 0.647$ (Kratter & Lodato 2016). Observationally speaking, the disk thickness is important when we take into account the molecular line emission of CO isotopologues. As a matter of fact, $^{12}$CO, which is the most abundant isotopologue, is not a good tracer of the disk midplane, because it becomes optically thick at the disk surface. On the contrary, other less abundant CO isotopologues such as $^{13}$CO or $^{18}$O have more optically thin lines and as a consequence they trace better the disk midplane (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 (which we do not do in this paper), as shown in Dipierro et al. (2014), 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 a smoothed particles hydrodynamics (SPH) simulation using the code PHANTOM (Price et al. 2018b). This code is widely used in the astrophysical community to study gas and dust dynamics in accretion disks (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 disk, neglecting the dust component. The initial conditions of the disk are $r_{in} = 1$ au, $r_{out} = 50$ au, $\Sigma \propto r^{-1}$, $M_d = 0.5$, and $M_*$, $M_* = 1 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
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison between analytical (black line) and numerical (green line) perturbation for different viewing angles of the simulated GI disk. Each picture is rotated along the z-axis of an angle of π/2.</figcaption></figure>

shows a predominance of the *m* = 2 azimuthal mode, and the computed pitch angle is αₚ ~ 13°.

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 disk mass that best describes the curve using Equation (2): the best value corresponds to $M_d = 0.5M_\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 wiggle<fnref n="6" /> is $A_{sim} = 0.11$ rad. The amplitude of the wiggle as a function of the cooling factor, for the parameters previously reported, is described by $A(\beta) = A_{in}\beta^{-0.5}$, where $A_{in} = 0.33$rad. The estimated cooling factor is then $\beta_{sim} = (A_{sim}/A_{in})^{-2} \simeq 9$, which 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 behavior is visible in the comparison with the numerical simulations, hence the simulated perturbation appears less wide than the analytical one.

Figure 4 shows a comparison between the analytical and the numerical wiggle, for different viewing angles of the spiral structure (the observer is assumed to be along a vertical line at the bottom of the images). There is a very good agreement in panels (a) and (c) for which the line of sight intercepts the largest extent of the spiral, while the perturbation is overestimated in panels (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 disk (where the perturbation is however smaller), our model overestimates the perturbation in the outer disk, 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 disk 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-1 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 disk. The kinematics deviations, called GI wiggles, depend on the

<footnote n="6">⁶ Here, we refer to panel (a) of Figure 4.</footnote>
structure of the spiral density wave, namely, its amplitude (connected to the cooling and to the disk 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 disks: 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 disk gravitational potential, meaning that the effects of disk self-gravity are not negligible.

Higher-resolution observations of systems like those will make it possible to investigate the cooling of protoplanetary disks: 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.

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 program under the Marie Skłodowska-Curie grant agreement No. 823823 (DUSTBUSTERS RISE project). B.V. acknowledges funding from the ERC CoG project PODCAST No 864965.

**Facility:** SPLASH: an interactive visualization tool for SPH data (Price 2007).

## ORCID iDs

Cristiano Longarini https://orcid.org/0000-0003-4663-0318

Giuseppe Lodato https://orcid.org/0000-0002-2357-7692

Claudia Toci https://orcid.org/0000-0002-6958-4986

Benedetta Veronesi https://orcid.org/0000-0002-2611-7931

Cassandra Hall https://orcid.org/0000-0002-8138-0425

Ruobing Dong https://orcid.org/0000-0001-9290-7846

Jason Patrick Terry https://orcid.org/0000-0002-8590-7271

## References

* Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987

* Bertin, G., & Lin, C. C. 1996, Spiral Structure in Galaxies a Density Wave Theory (Cambridge, MA: MIT Press)

* Bertin, G., & Lodato, G. 1999, A&A, 350, 694

* Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, NJ: Princeton Press)

* Bollati, F., Lodato, G., Price, D. J., & Pinte, C. 2021, MNRAS, 504, 5444

* Cadman, J., Rice, K., Hall, C., Haworth, T. J., & Biller, B. 2020, MNRAS, 492, 5041

* Clarke, C. J. 2009, MNRAS, 396, 1066

* Cossins, P., Lodato, G., & Clarke, C. J. 2009, MNRAS, 393, 1157

* Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114

* Cuello, N., Louvet, F., Mentiplay, D., et al. 2020, MNRAS, 491, 504

* Dipierro, G., Lodato, G., Testi, L., & de Gregorio Monsalvo, I. 2014, MNRAS, 444, 1919

* Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73

* Dipierro, G., Ricci, L., Pérez, L., et al. 2018, MNRAS, 475, 5296

* Dong, R., & Fung, J. 2017, ApJ, 835, 38

* Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJL, 809, L5

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

* Hall, C., Dong, R., Teague, R., et al. 2020, ApJ, 904, 148

* Hall, C., Forgan, D., Rice, K., et al. 2016, MNRAS, 458, 306

* Hall, C., Forgan, D., & Rice, K. 2017, MNRAS, 470, 2517

* Hirose, S., & Shi, J.-M. 2017, MNRAS, 469, 561

* Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018, ApJL, 869, L43

* Johnson, B. M., & Gammie, C. F. 2003, ApJ, 597, 131

* Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271

* Kratter, K. M., & Matzner, C. D. 2006, MNRAS, 373, 1563

* Kratter, K. M., Murray-Clay, R. A., & Youdin, A. N. 2010, ApJ, 710, 1375

* Lau, Y. Y., & Bertin, G. 1978, ApJ, 226, 508

* Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646

* Lodato, G., & Rice, W. K. M. 2004, MNRAS, 351, 630

* Lodato, G., & Rice, W. K. M. 2005, MNRAS, 358, 1489

* Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96

* Oort, J. H. 1927, BAN, 3, 275

* Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467

* Paneque-Carreño, T., Pérez, L. M., Benisty, M., et al. 2021, ApJ, 914, 88

* Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Sci, 353, 1519

* Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJL, 860, L13

* Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJL, 890, L9

* Price, D. J. 2007, PASA, 24, 159

* Price, D. J., Cuello, N., Pinte, C., et al. 2018a, MNRAS, 477, 1270

* Price, D. J., Wurster, J., Tricco, T. S., et al. 2018b, PASA, 35, e031

* Ragusa, E., Alexander, R., Calcino, J., Hirsh, K., & Price, D. J. 2020, MNRAS, 499, 3362

* Rice, W. K. M., & Armitage, P. J. 2009, MNRAS, 396, 2228

* Rice, W. K. M., Armitage, P. J., Mamatsashvili, G. R., Lodato, G., & Clarke, C. J. 2011, MNRAS, 418, 1356

* Rice, W. K. M., Lodato, G., & Armitage, P. J. 2005, MNRAS, 364, L56

* Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner, D. J., & Qi, C. 2013, ApJ, 774, 16

* Rosotti, G. P., Benisty, M., Juhász, A., et al. 2020, MNRAS, 491, 1335

* Stamatellos, D., Hubber, D. A., & Whitworth, A. P. 2007, MNRAS, 382, L30

* Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 400, 1563

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

* Terry, J. P., Hall, C., Longarini, C., et al. 2021, arXiv:2108.11282

* Toci, C., Lodato, G., Christiaens, V., et al. 2020a, MNRAS, 499, 2015

* Toci, C., Lodato, G., Fedele, D., Testi, L., & Pinte, C. 2020b, ApJL, 888, L4

* Veronesi, B., Lodato, G., Dipierro, G., et al. 2019, MNRAS, 489, 3758

* Veronesi, B., Paneque-Carreño, T., Lodato, G., et al. 2021, ApJL, 914, L27

* Veronesi, B., Ragusa, E., Lodato, G., et al. 2020, MNRAS, 495, 1913