Investigating Protoplanetary Disk Cooling through Kinematics: Analytical GI Wiggle

The gravitational potential of a self-gravitating protostellar disk is not exactly Keplerian: indeed, when M_d ≳ (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

$$\begin{eqnarray}&&{{\rm{\Omega }}}^{2}=\displaystyle \frac{{{GM}}_{\star }}{{r}^{3}}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial {{\rm{\Phi }}}_{d}}{\partial r},\end{eqnarray} \tag{ 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

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {{\rm{\Phi }}}_{d}}{\partial r}(r,z) & = & \displaystyle \frac{G}{r}{\displaystyle \int }_{0}^{\infty }{{dr}}^{{\prime} }\left[K(\zeta )-\displaystyle \frac{1}{4}\left(\displaystyle \frac{{\zeta }^{2}}{1-{\zeta }^{2}}\right)\right.\\ & & \left.\times \left(\displaystyle \frac{{r}^{{\prime} }}{r}-\displaystyle \frac{r}{{r}^{{\prime} }}+\displaystyle \frac{{z}^{2}}{{{rr}}^{{\prime} }}\right)E(\zeta )\right]\sqrt{\displaystyle \frac{{r}^{{\prime} }}{r}}\zeta {\rm{\Sigma }}\left({r}^{{\prime} }\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where E(ζ) and K(ζ) are the complete elliptic integrals of the first kind, ${\zeta }^{2}=4{{rr}}^{{\prime} }/[{\left(r+{r}^{{\prime} }\right)}^{2}+{z}^{2}]$, and Σ 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.

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

$$\begin{eqnarray}\begin{array}{rcl}D(\omega ,k,m) & = & {\left(\omega -m{\rm{\Omega }}\right)}^{2}-{c}_{{\rm{s}}}^{2}{k}^{2}\\ & & +2\pi G{\rm{\Sigma }}| k| -{\kappa }^{2}=0,\end{array}\end{eqnarray} \tag{ 3 }$$

where ω is the wave angular frequency, k is the radial wavenumber, m is the azimuthal wavenumber, κ(r) is the epicyclic frequency, Ω(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

$$\begin{eqnarray}&&Q=\displaystyle \frac{{c}_{s}\kappa }{\pi G{\rm{\Sigma }}},\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{k}_{\mathrm{uns}}={H}_{\mathrm{sg}}^{-1}=\displaystyle \frac{\pi G{\rm{\Sigma }}}{{c}_{s}^{2}},\end{eqnarray} \tag{ 5 }$$

where ${H}_{\mathrm{sg}}={c}_{s}^{2}/\pi G{\rm{\Sigma }}$ is the disk thickness in the self-gravitating case.

A marginally stable disk is characterized by having Q ≃ 1: in this case, the only expected excited modes would be k ≃ k_uns. Thus, Equation (3) tells us that Ω ≃ ω/m = Ω_p, where Ω_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

$$\begin{eqnarray}&&m\phi +\psi (r)=\mathrm{const},\quad \mathrm{mod}2\pi ,\end{eqnarray} \tag{ 6 }$$

where ψ(r) is the shape function and it holds that k = d ψ/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 α_p , hereafter referred to as the pitch angle; it is given by

$$\begin{eqnarray}&&\tan {\alpha }_{p}={\left|r\displaystyle \frac{\partial \phi }{\partial r}\right|}^{-1},\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{1}(r,\phi ,t)=\mathrm{Re}\left[\delta {\rm{\Sigma }}(r){e}^{i(m\phi -\omega t+\psi (r))}\right],\end{eqnarray} \tag{ 8 }$$

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

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{1}(r,\phi ,t)=-\displaystyle \frac{2\pi G}{| k| }{{\rm{\Sigma }}}_{1}(r,\phi ,t).\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{q}^{-}=-\displaystyle \frac{e}{{t}_{\mathrm{cool}}},\end{eqnarray} \tag{ 10 }$$

where e is the internal energy per unit mass and all details of the cooling are absorbed by t_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 ={\rm{\Omega }}{t}_{\mathrm{cool}}=\mathrm{const}$ (Gammie 2001). In this work, we use this prescription and hereafter we will refer to the cooling process in terms of β.

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 ≫ 1) will cool down due to radiative processes, until eventually reaching the marginally stable state Q ≃ 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 ≃ 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

$$\begin{eqnarray}&&\displaystyle \frac{\delta {\rm{\Sigma }}}{{\rm{\Sigma }}}=\chi {\beta }^{-1/2},\end{eqnarray} \tag{ 11 }$$

where χ is the proportionality factor, which is of the order of unity, and β is the adimensional cooling factor. Hence, a more efficient cooling (lower β) 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 δΣ/Σ, the β factor should be overestimated.

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

$$\begin{eqnarray}&&\displaystyle \frac{\kappa -{{\rm{\Omega }}}_{k}}{{{\rm{\Omega }}}_{k}}\lt 1,\quad \displaystyle \frac{{\rm{\Omega }}-{{\rm{\Omega }}}_{k}}{{{\rm{\Omega }}}_{k}}\lt 1,\end{eqnarray} \tag{ 21 }$$

where with the subscript k we identify the Keplerian quantities. This assumption allows us to write κ ≃ Ω ≃ Ω_k ∝ r^−3/2 and B ≃ −Ω/4. With these assumptions, the equations for the VPs are

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}\delta {u}_{r} & = & 2{im}\chi {\beta }^{-1/2}{\left(\displaystyle \frac{{M}_{d}(r)}{{M}_{\star }}\right)}^{2}{u}_{k},\\ \delta {u}_{\phi } & = & -\displaystyle \frac{i\chi {\beta }^{-1/2}}{2}\left(\displaystyle \frac{{M}_{d}(r)}{{M}_{\star }}\right){u}_{k},\end{array}\end{array}\end{eqnarray} \tag{ 22 }$$

where u_k is the Keplerian speed, and then the velocity field is given by ${u}_{r}=\mathrm{Re}\left[\delta {u}_{r}{e}^{i(m\phi -\omega t+\psi )}\right]$ and ${u}_{\phi }=r{\rm{\Omega }}\,+\mathrm{Re}\left[\delta {u}_{\phi }{e}^{i(m\phi -\omega t+\psi )}\right]$, where Ω is given by Equation (1). To obtain these expressions we used that Q = 1 and that d ψ/dr = k = k_uns. Note that δ u_r /δ u_ϕ ≃ 4mM_d /M_⋆: for example, when m = 2, a relatively light disk having M_d = 0.125 M_⋆ has δ u_r = δ u_ϕ .

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²/r²) (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.

So far we have considered only the case of constant β − cooling. In principle, however, self-consistent models of GI disks (Clarke 2009; Rice & Armitage 2009) show that β 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 chose any cooling law β(r, ρ(r), T(r)) and then obtain the VPs through Equation (20).

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_syst = 0, corresponds to v_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 factor β ), the opening angle α_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, ϕ) 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

$$\begin{eqnarray}&&{C}_{v}^{k}(s)=\left(\begin{array}{c}{f}_{r}^{k}(s)\\ {f}_{\phi }^{k}(s)\end{array}\right),\quad {C}_{v}^{p}(s)=\left(\begin{array}{c}{f}_{r}^{p}(s)\\ {f}_{\phi }^{p}(s)\end{array}\right),\end{eqnarray} \tag{ 24 }$$

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

$$\begin{eqnarray}&&{{ \mathcal A }}_{v}={\left[{\int }_{0}^{1}{ds}| | {C}_{v}^{p}-{C}_{v}^{k}| {| }^{2}\right]}^{1/2},\end{eqnarray} \tag{ 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 β generates a bigger deviation from the Keplerian case, because the amplitude of the density perturbations is inversely proportional to β (Equation (11)). In panel (b) of Figure 2 we show the amplitude of the wiggle (v_obs = 0) as a function of the cooling parameter β: it is possible to describe the relation between ${ \mathcal A }$ and β 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 ${ \mathcal 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 ≃ 1 we get c_s = u_k H/r ≃ 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 ${ \mathcal A }\propto \sqrt{{c}_{1}{M}_{d}^{2}+{c}_{2}{M}_{d}^{4}}$, where c₁ and c₂ are two constants that depend on the other parameters, such as the cooling.

Zoom In Zoom Out Reset image size

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_obs = v_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 β 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 ≃ 0.08 M_⊙ ≃ 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 ϕ = π/2 (semiminor axis of the disk) we get

$$\begin{eqnarray}\begin{array}{c}\begin{array}{c}\begin{array}{rcl}{v}_{{\rm{obs}}} & = & {u}_{r}\sin \theta ={\rm{Re}}\left[\delta {u}_{r}(r){e}^{i(m\pi /2+\psi )}\right]\\ \sin \theta & = & {f}_{1}(r){\beta }^{-1/2}{\left(\displaystyle \frac{{M}_{d}}{{M}_{\star }}\right)}^{2}{u}_{k},\end{array}\end{array}\end{array}\end{eqnarray} \tag{ SPH26 }$$

where f₁ is a known function of radius; conversely, for ϕ = 0 (semimajor axis of the disk) we get

$$\begin{eqnarray}\begin{array}{c}\begin{array}{c}\begin{array}{rcl}{v}_{{\rm{obs}}} & = & {u}_{\phi }\sin \theta =\left(r{\rm{\Omega }}+{\rm{Re}}\left[\delta {u}_{\phi }(r){e}^{i\psi }\right]\right)\\ \sin \theta & = & {f}_{2}(r){\beta }^{-1/2}\displaystyle \frac{{M}_{d}}{{M}_{\star }}{u}_{k}+r{\rm{\Omega }}\sin \theta ,\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 27 }$$

where f₂ 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 β, 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 α_p and m. We clearly see that the frequency of the wiggle is bigger when decreasing α_p and increasing m.

Zoom In Zoom Out Reset image size

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 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 ≠ 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

$$\begin{eqnarray}&&{ \mathcal J }=m\displaystyle \frac{\pi G{\rm{\Sigma }}}{r{\kappa }^{2}}\displaystyle \frac{4{\rm{\Omega }}}{\kappa }{\left|\displaystyle \frac{d\mathrm{ln}{\rm{\Omega }}}{d\mathrm{ln}r}\right|}^{1/2}\approx \sqrt{6}m\displaystyle \frac{{M}_{d}}{{M}_{* }}.\end{eqnarray} \tag{ 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

$$\begin{eqnarray}\begin{array}{rcl}{D}_{\mathrm{thick}}(\omega ,k,m) & = & {\left(\omega -m{\rm{\Omega }}\right)}^{2}-{c}_{s}^{2}{k}^{2}\\ & & -{\kappa }^{2}+2\pi G{\rm{\Sigma }}\displaystyle \frac{| k| }{1+| k| H}.\end{array}\end{eqnarray} \tag{ 29 }$$

In this case, the stability criterion only changes slightly, becoming Q ≳ 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, ¹²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 ¹³CO or C¹⁸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.

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, Σ ∝ r⁻¹, M_d = 0.5, and M_⋆, M_⋆ = 1 M_⊙. The cooling factor β has been set to β = 8 and the two parameters that control the artificial viscosity to α^AV = 0.1, β^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 ∝ r^−0.25. However, this profile is rapidly modified by the cooling. For this simulation, we used N = 5 × 10⁵ particles of gas. The simulation shows a predominance of the m = 2 azimuthal mode, and the computed pitch angle is α_p ∼ 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_ϕ (r, ϕ). 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_⋆, as expected. We have broken the degeneracy between the mass and the cooling: thus, we are now able to constrain the cooling β 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 ⁶ is ${{ \mathcal A }}_{\mathrm{sim}}=0.11\,\mathrm{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 }}_{\mathrm{in}}{\beta }^{-0.5}$, where ${{ \mathcal A }}_{\mathrm{in}}=0.33\mathrm{rad}$. The estimated cooling factor is then ${\beta }_{\mathrm{sim}}={\left({{ \mathcal A }}_{\mathrm{sim}}/{{ \mathcal A }}_{\mathrm{in}}\right)}^{-2}\simeq 9$, which is in good agreement with the real value β = 8. The overestimated value of β 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.

Zoom In Zoom Out Reset image size

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.

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 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 β 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.