Angular momentum transport via gravitational instability in the Elias 2-27 disc

A disc is in the gravito-turbulent regime when its angular momentum transport is driven by gravitational instability. In this context, the $R\phi$ component of the vertically integrated stress tensor can be written as:

where $g_{R},g_{\phi}$ are the radial and azimuthal component of the perturbed self-gravitating field, and the brackets indicate azimuthal averaging. To this stress, it is necessary to also add the induced Reynolds stress, given by:

where $\delta v_{R}$ and $\delta v_{\phi}$ are the perturbed fluid velocities (Balbus & Papaloizou 1999). In the classic viscous scenario, this term is responsible for angular momentum transport in the accretion disc. The Shakura and Sunyaev $\alpha-$prescription (Shakura & Sunyaev 1973) relates the stress tensor to the local disc pressure

where, for a Keplerian disc, ${\text{d}\log\Omega_{k}}/{\text{d}\log R}=-3/2$. It is possible to show that the transport of energy and angular momentum through the propagation of the GI spiral density waves can be divided into two parts (Toomre 1969; Shu 1970): a non-local term and a viscous-like term. Since the non-local term is important only for very high disc-to-star mass ratios, $M_{d}/M_{\star}\gtrsim 0.5$ (Lodato & Rice 2005), gravitationally unstable protostellar discs essentially behave as $\alpha$-discs (Lodato & Rice 2004; Forgan et al. 2011).

To characterise the transport of angular momentum through spiral density waves, we need to rely on numerical simulations of gravitationally unstable discs (Cossins et al. 2009). Usually in numerical simulations of gravitationally unstable discs, the $\beta$ cooling framework is adopted (Gammie 2001). We suppose that the disc is cooling with a rate per unit mass of:

where $e=c_{s}^{2}/\gamma(\gamma-1)$ is the internal energy per unit mass, with the adiabatic index $\gamma=5/3$, and the cooling time is defined in terms of the dynamical one: $t_{\text{cool}}=\beta\Omega^{-1}$. In the absence of external heating mechanisms, an initially stable hot disc ($Q>>1$) will cool down, eventually reaching the marginally stable state ($Q=1$). At this point, gravitational instability turns on: the disc develops a spiral structure that, by means of compression and shocks, leads to 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 (Kratter & Lodato 2016). In this regime, namely, thermal saturation, the cooling is completely balanced by heating provided by the shocks, and the amplitude of spiral perturbations saturates at a fixed value, according to Cossins et al. (2009):

where $\epsilon$ is the heating factor, $\mathcal{M}$ and $\tilde{\mathcal{M}}$ are the Mach numbers relative to the radial phase speed and Doppler-shifted radial phase speed of the wave: respectively, $m\Omega_{p}/k$ and $m(\Omega-\Omega_{p})/k$, with $\Omega_{p}$ as the spiral pattern frequency. Cossins et al. (2009) showed through numerical simulations that the relevant scaling of the spiral density perturbation is $\beta^{-1/2}$ and the other terms of the order of unity. For this reason, we introduce $\chi$, which is of the order of unity. In the thermal saturation regime, the transport of angular momentum provided by gravitational instability is described within an $\alpha-$framework(Kratter & Lodato 2016):

Combining Eqs. (7) and (6), we obtain the relationship between the amplitude of the density wave and the amount of angular momentum

where the last equation is valid for a Keplerian disc. In the following, we assume that the constant of proportionality is $\chi=1$, which is well justified by numerical simulations (Cossins et al. 2009).

Kinematics offers a unique opportunity to quantify the transport of angular momentum in a gravitationally unstable disc. Indeed, when there is a spiral density wave, also the velocity field is affected. Hall et al. (2020) predicted that a disc undergoing such instability has clear kinematic signatures in molecular line observations across the entire disc azimuth and radius, called ‘GI wiggles’. Longarini et al. (2021) provided an analytical model to describe such kinematic signatures, implemented in the publicly available code giggle¹¹1https://doi.org/10.5281/zenodo.10205110. Under the hypothesis of thin disc, in a marginally unstable regime $(Q=1)$ and in thermal saturation, the amplitude of the velocity perturbation increases with the disc-to-star mass ratio $(M_{d}/M_{\star})$, as has also been described with numerical simulations (Terry et al. 2022), and which is proportional to the cooling factor, $\beta^{-1/2}$. However, from the analysis carried out in Longarini et al. (2021), the actual quantity that determines the amplitude of the velocity perturbations $(\delta u_{R},\delta u_{\phi})$ is the amplitude of the spiral density wave, $\delta\Sigma/\Sigma$. This quantity is intrinsically linked to the efficiency of angular momentum transported by the spiral, which can be described within an $\alpha-$viscosity framework through Eq. (7). Thus, it is possible to constrain the value of $\alpha_{\text{GI}}$ from the amplitude of the wiggle. Rewriting the velocity perturbations (Eq. 22 of Longarini et al. (2021)) as a function of $\alpha_{\text{GI}}$, we obtain:

The final expressions for the velocity field are:

Veronesi et al. (2021) estimated the dynamical mass of Elias 2-27 from ¹³CO and C¹⁸O rotation curves, and found $M_{\star}=0.46\text{M}_{\odot}\pm 0.03$ and $M_{d}=0.08\pm 0.04\text{M}_{\odot}$. In their fit, they used a self-similar surface density given by

where the scale radius has been fixed to $R_{c}=200\text{au}$. As for the thermal structure, they assumed a vertically isothermal disc with $T(R)=20(R/60\text{au})^{0.5}$K taken from Pérez et al. (2016), which corresponds to a sound speed of $c_{s}=281\text{m/s}(R/60\text{au})^{-0.25}$. With this information, it is possible to compute the Toomre parameter profile for Elias 2-27, as displayed in Fig. 1. Despite not being exactly $Q=1$, the Toomre profile is close enough to the critical threshold to consider gravitational instability to be significant. In addition, we are not considering uncertainties on $T(R)$, that can impact on the $Q-$parameter estimate. Finally, we know that the disc shows non-axisymmetric features in gas and dust, making the azimuthally averaged Toomre profile solely an approximation of the actual value.

In this work, we use the ¹³CO $J=3-2$ datacube presented in Paneque-Carreño et al. (2021). The images have been obtained with a robust parameter of 0.5, resulting in a beam size of $0.26^{\prime\prime}\times 0.25^{\prime\prime}$ for the ¹³CO and a spectral resolution of $\Delta v=111$m/s. Further details of the observations and reduction can be found in Paneque-Carreño et al. (2021). In this work, we use a Gaussian velocity map obtained with bettermoments (Teague & Foreman-Mackey 2018). The code also returns a map of the errors on the velocity field (Teague 2019). In the analysis, we masked the emission coming from the first two beams ($\sim 50\text{au}$).

In the analytical model for the GI wiggle of Longarini et al. (2021), the amplitude of the velocity perturbations is determined by the disc to star mass ratio and the cooling factor. There is a degeneracy between the two quantities; however, in the case of Elias 2-27 the value of the disc to star mass ratio is known (Veronesi et al. 2021). Longarini et al. (2021) showed that the amplitude of the wiggle in the PP (position-position) space scales as $\beta^{-1/2}$ (i.e. $\alpha_{\text{GI}}^{1/2}$), and they proposed this relationship as a way to constrain this unknown parameter. In this paragraph, we study the wiggle in the PV (position-velocity) space, as done in Speedie et al. (2024, submitted). We consider the observed velocity field:

where $u_{\phi}$ and $u_{R}$ are described by the model presented in Eqs. (11) and (12) and $i$ is the disc inclination. We compute the variation of velocity along the semi-minor axis of the disc (i.e. $\phi=-\pi/2$), because along this axis only the radial velocity contributes to the observed velocity field

where $\psi(R)$ is the phase function of the spiral that is given by

with $\alpha_{p}$ being the pitch angle of the spiral and $\psi_{0}$ is just a phase shift. We suppose that the pitch angle is constant over the radial extent of the disc, that is well justified for GI spirals (Cossins et al. 2009). We note that the last equation can be written as a function of $\alpha_{\text{GI}}$ as:

We then extract the PV wiggle by cutting along the semi minor axis and considering just the southern part of the disc, on account of cloud contamination (Paneque-Carreño et al. 2021). The errors on the velocity are returned by (Teague & Foreman-Mackey 2018) as shown in figure 2. The disc and spiral parameters we use are taken from literature, namely, $M_{\star}=0.46\text{M}_{\odot}$, $M_{d}=0.08\text{M}_{\odot}$, $R_{c}=200$au (Veronesi et al. 2021), $\alpha_{p}=13^{\circ}$, $i=56.2^{\circ}$, $\text{PA}=118.8^{\circ}$ (Paneque-Carreño et al. 2021), and $m=2$. Hence, the only free parameters in Eq. (17) are $\alpha_{\text{GI}}$ and $\psi_{0}$. To fit the curve to the data, we used the method of nonlinear least squares implemented in Scipy (Virtanen et al. 2020). While the spiral pattern in the continuum extends to approximately $200\text{au}$, non-axisymmetric kinematic signatures in the gas emission are observable across the entire radial extent of the disk. Hence, in the fitting procedure, we analysed the signal up to the outer edge of the disc. The best-fit values are $\alpha_{\text{GI}}=0.038\pm 0.018$ and $\psi_{0}=43^{\circ}\pm 1^{\circ}$. The value of $\alpha_{\text{GI}}$ corresponds to a $\beta=10.5$, meaning that the cooling time of the system is approximately ten times the dynamical one. We estimated the error on $\alpha_{\text{GI}}$ propagating the uncertainties on star and disc masses²²2The fitting procedure implemented in Scipy returns an error for the best fit parameters. As for the $\alpha_{\text{GI}}$ parameter, the error provided by the fit is subdominant compared to the one driven by the star and disc masses uncertainties..

Figure 2 shows the comparison between the extracted PV wiggle from ¹³CO data and the model that best describes the data. By comparing the model and the data, we observe that the overall shape of the perturbation is well reproduced. However, it appears that there is a radial shift between the two curves. This effect can be attributed to the fact that in the analytical model, we assume that a single spiral mode is present, with $m=2$, and we also suppose that the perturbation wave-number, $k,$ is the most unstable $(k_{\text{uns}})$. While the $m=2$ mode, with $k=k_{\text{uns}}$ may be the dominant one, determining the overall morphology of the spiral, this does not prevent the presence of additional lower amplitude modes that will interfere with the dominant one and thus create a more complex pattern than the purely sinusoidal one that we assume here. The scope of this work is to find the amount of stress generated by gravitational instability; hence, we are interested in the amplitude of the wiggle, not in the whole shape. For this reason, the additional lower amplitude modes are not a concern in terms of the amplitude of the perturbation.

Figure 3 shows the analytical observed velocity field with the aforementioned parameters. Despite its simplicity, the analytical model matches the shape of the different channels of the data. In particular, the bending of the isovelocity contours is described very well by the analytical model, especially in the blue-shifted region, where the effect of cloud contamination are negligible. The analytical models are produced using the publicly available code giggle. It is important to point out that the model presented in Longarini et al. (2021) is limited to two dimensions and does not consider the vertical extent of the disc. There is evidence suggesting that the ¹³CO emission in Elias 2-27 is optically thick (Paneque-Carreño et al. 2022); consequently, the received signal does not originate from the midplane, but from a layer at $z\neq 0$. The analytical model we employ is not able to reproduce this effect, and, for a proper comparison, hydrodynamical and radiative transfer simulations would be needed. In any case, assuming that $z/R<0.3$ for the ¹³CO (Paneque-Carreño et al. 2022), the difference in inclination induced by the finite height of the emitting layer is $\delta i=\arcsin(0.3)\simeq 0.29$. The corresponding geometrical error is roughly $1-\sin(i+\delta i)/\sin(i)\simeq 15\%$, that is subdominant compared to the one driven by the mass estimate.

For a gravitationally unstable disc, kinematics offers a unique opportunity to quantify the transport of angular momentum of a gravitationally unstable disc. We show that gravitational instability in Elias 2-27 transports angular momentum and the equivalent $\alpha-$viscosity parameter is $\alpha_{\text{GI}}=0.038\pm 0.018$. Since the effective viscosity is responsible for the accretion process, it is useful to predict the expected accretion rate onto the central object, and compare it to the observed one.

According to the self-similar solution (Lynden-Bell & Pringle 1974), the surface density and the accretion rate of the disc can be written as:

where we have supposed that $\nu=\nu_{c}(R/R_{c})$, where the subscript $c$ means that the quantity is evaluated at the scale radius $R_{c}$. Within a $\alpha-$viscosity framework, the kinematic viscosity, $\nu,$ is

We note that within our assumption that $T\propto R^{0.5}$, for a constant $\alpha$ the kinematic viscosity scales as $\nu\propto R$. In this way, we can write the accretion rate onto the central object as the limit for $R\to 0$ of Eq. (19):

where $\Omega_{c}=\Omega(R_{c})=\sqrt{GM_{\star}/R_{c}^{3}}$.

Using $\alpha_{\text{GI}}=0.038$ and the disc parameters described in the previous paragraph, it is possible to compute the accretion rate onto the central object by using Eq. (21). Thus, we obtain:

where the error has been computed through propagation from the errors in $\alpha_{\text{GI}}$. The model for the accretion rate of Elias 2-27 reproduces the one measured by Natta et al. (2006) very well, namely, $\log_{10}\dot{M}_{\star}[\text{M}_{\odot}/\text{yr}]=-7.2\pm 0.5$. They used J and K-band spectra to derive the mass accretion rate of objects in the $\rho-$Ophiuchi star forming region from the intensity of the hydrogen recombination lines. More recently, a new estimate of Elias 2-27 accretion rate was provided by Testi et al. (2022), being $\log_{10}\dot{M}_{\star}[\text{M}_{\odot}/\text{yr}]=-7.3$, which is consistent with the previous measurement.

To constrain $\alpha_{\text{GI}}$ in Elias 2-27, we made the strong hypothesis that this quantity, or equivalently the cooling $\beta$, is constant throughout the disc. We are aware that, realistically, this is not the case. However, we point out that our disc model is self consistent. Indeed, assuming a self similar solution for the surface density with $\Sigma\propto R^{-1}$, we are imposing that the kinematic viscosity scales as $\nu\propto R$. Since $\nu=\alpha c_{s}H$, with $c_{s}\propto R^{-0.25}$ (Pérez et al. 2016) and $H\propto R^{1.25}$, the viscosity coefficient, $\alpha,$ should be constant with radius. We also point out that, in order to obtain an estimate of the amount of angular momentum transported within the disc from the GI wiggle, an assumption should be made with respect to how $\alpha$ varies with the radius. A more realistic disc model for Elias 2-27 is not within the scope of this work. We stress that the choice of a self-similar profile for the surface density was made to maintain consistency with the work of Veronesi et al. (2021).

Elias 2-27 is a young system, and its interactions with the surrounding environment are possibly perturbing the disc. Specifically, the disc is partially embedded within the molecular cloud, which absorbs ¹²CO and a portion of the ¹³CO emission (Pérez et al. 2016; Paneque-Carreño et al. 2021), then feeds the disc with mass. In particular, infall is an alternative way to trigger gravitational instability (Kratter & Matzner 2006; Kratter et al. 2008; Kratter & Lodato 2016). Gravitational instability occurs when the Toomre parameter is of the order of unity, and this threshold can be reached by cooling the disc (i.e. decreasing the sound speed) or by adding mass (i.e. increasing the surface density). When GI is triggered by infall, there is a mechanism akin to the thermal saturation (Kratter & Lodato 2016).

Our $\alpha_{\text{GI}}$ estimate remains agnostic to the source of self-regulation, whether it arises from cooling or from the addition of mass to the disc. Indeed, what we are measuring through the wiggle is the amplitude of the surface density perturbation, $\delta\Sigma/\Sigma$. It is intrinsically linked to the efficiency of angular momentum that is transported by the spiral, regardless of the origin of the instability.

The value of $\beta$ we measure is much higher than the threshold for disc fragmentation into bound gas clumps (Gammie 2001; Deng et al. 2017), implying that (while the disc is gravitationally unstable) direct planet formation through gravitational instability is unlikely. However, many studies (Paardekooper 2012; Young & Clarke 2015) showed that stochastic fragmentation of gas in spiral arms can happen for high values of $\beta$. As for the solid component, Longarini et al. (2023a, b); Rowther et al. (2024) investigated the possibility of forming planetary cores in gravitationally unstable discs through dust collapse. They found that for a sufficiently long cooling time, $\beta>10,$ and high disc-to-star mass ratio, $M_{d}/M_{\star}\sim 0.2$, dust efficiently collects inside spiral arms and its dispersion velocity is so low to induce collapse into bound objects with a mass of $\sim 10\text{M}_{\oplus}$. The inferred disc-to-star mass ratio (Veronesi et al. 2021) and the cooling time for Elias 2-27 make it a perfect candidate for planet formation through dust collapse.

Huang et al. (2018) characterised annular substructures in the discs within the DSHARP sample and found that Elias 2-27 has a gap at $R_{g}=69.1\pm 0.4$au with a width of $\Delta=14.3\pm 1.1$au. Although several mechanisms can explain the origin of gaps in protoplanetary discs, a common explanation is planet disc interaction. Under the planetary interpretation, the width of the gap scales as the Hill radius of the planet, defined as:

where $M_{p}$ is the mass of the protoplanet. The last relation has been obtained by averaging results from hydrodynamical simulations. Following Lodato et al. (2019), the relation between the gap width and the Hill radius is $\Delta=5.5R_{h}$, that translates into

Using the gap width and location of Huang et al. (2018), and the star mass of Veronesi et al. (2021), the inferred mass of the protoplanet is $M_{p}=24\pm 6\text{M}_{\oplus}$. This result is in good agreement with the mass range of Longarini et al. (2023a, b). Another element that points towards the dust collapse is the value of the Toomre parameter. As shown in Longarini et al. (2023a), when the gravitational instability is driven by the cold component (dust in this case), the critical value of the Toomre parameter is $>1$, as observed in Elias 2-27 (see Fig. 1).