# Spiral Arms and a Massive Dust Disk with Non-Keplerian Kinematics: Possible Evidence for Gravitational Instability in the Disk of Elias 2-27

T. Paneque-Carreño, L. M. Pérez, M. Benisty, C. Hall, B. Veronesi, G.

Lodato, A. Sierra, J. M. Carpenter, S. M. Andrews, Jaehan Bae, et al.

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

► To cite this version:

T. Paneque-Carreño, L. M. Pérez, M. Benisty, C. Hall, B. Veronesi, et al. Spiral Arms and a Massive Dust Disk with Non-Keplerian Kinematics: Possible Evidence for Gravitational Instability in the Disk of Elias 2-27. The Astrophysical Journal, 2021, 914, 10.3847/1538-4357/abf243. insu-03705307

HAL Id: insu-03705307

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

Submitted on 2 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_2.png" /></figure>

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

# Spiral Arms and a Massive Dust Disk with Non-Keplerian Kinematics: Possible Evidence for Gravitational Instability in the Disk of Elias 2–27

T. Paneque-Carreño¹, L. M. Pérez¹, M. Benisty¹,²,³, C. Hall⁴,⁵,⁶, B. Veronesi⁷, G. Lodato⁷, A. Sierra¹, J. M. Carpenter⁸, S. M. Andrews⁹, Jaehan Bae¹⁰,<fnref n="20" />, Th. Henning¹¹, W. Kwon¹²,¹³, H. Linz¹¹, L. Loinard¹⁴,¹⁵, C. Pinte³,¹⁶, L. Ricci¹⁷, M. Tazzari¹⁸, L. Testi¹⁹, and D. Wilner⁹

¹ Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile; teresapaz.paneque@gmail.com

² Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS/INSU UMI 3386, Chile

³ Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

⁴ Department of Physics and Astronomy, The University of Georgia, Athens, GA 30602, USA

⁵ Center for Simulational Physics, The University of Georgia, Athens, GA 30602, USA

⁶ School of Physics & Astronomy, University of Leicester, University Road, Leicester, LE1 7RH, UK

⁷ Dipartimento di Fisica, Universita degli Studi di Milano, Via Celoria, 16, Milano, I-20133, Italy

⁸ Joint ALMA Observatory, Avenida Alonso de Córdova 3107, Vitacura, Santiago, Chile

⁹ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

¹⁰ Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA

¹¹ Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

¹² Department of Earth Science Education, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea

¹³ Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea

¹⁴ Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Morelia, 58089, México

¹⁵ Instituto de Astronomía, Universidad Nacional Autónoma de México, Apartado Postal 70-264, Ciudad de México 04510, México

¹⁶ School of Physics and Astronomy, Monash University, Clayton, VIC 3800, Australia

¹⁷ Department of Physics and Astronomy, California State University Northridge, 18111 Nordhoff Street, Northridge, CA 91330, USA

¹⁸ Institute of Astronomy, University of Cambridge, Madingley Road, CB3 0HA, Cambridge, UK

¹⁹ European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching bei München, Germany

Received 2020 December 30; revised 2021 March 24; accepted 2021 March 24; published 2021 June 17

## Abstract

To determine the origin of the spiral structure observed in the dust continuum emission of Elias 2–27 we analyze multiwavelength continuum ALMA data with a resolution of ~0.1/2 (~23 au) at 0.89, 1.3, and 3.3 mm. We also study the kinematics of the disk with ¹³CO and C¹⁸O ALMA observations in the *J* = 3–2 transition. The spiral arm morphology is recovered at all wavelengths in the dust continuum observations, where we measure contrast and spectral index variations along the spiral arms and detect subtle dust-trapping signatures. We determine that the emission from the midplane is cold and interpret the optical depth results as signatures of a disk mass higher than previous constraints. From the gas data, we search for deviations from Keplerian motion and trace the morphology of the emitting surfaces and the velocity profiles. We find an azimuthally varying emission layer height in the system, large-scale emission surrounding the disk, and strong perturbations in the channel maps, colocated with the spirals. Additionally, we develop multigrain dust and gas hydrodynamical simulations of a gravitationally unstable disk and compare them to the observations. Given the large-scale emission and highly perturbed gas structure, together with the comparison of continuum observations to theoretical predictions, we propose infall-triggered gravitational instabilities as the origin for the observed spiral structure.

Unified Astronomy Thesaurus concepts: Planet formation (1241); Planetary system formation (1257); Protoplanetary disks (1300); Hydrodynamical simulations (767); Gravitational instability (668)

## 1. Introduction

Protoplanetary disks around young stars have shown various structures in their thermal dust continuum emission. Observing and understanding their origin is necessary to understanding the chemical, physical, and dynamical processes that are ongoing in a protoplanetary disk. Of the structures present in dust emission, rings and gaps are the most common (e.g., ALMA Partnership et al. 2015; Andrews et al. 2018; Fedele et al. 2018; Huang et al. 2018b; Long et al. 2018), albeit arcs and spirals have also been observed in some systems (e.g., Dong et al. 2018; Huang et al. 2018c).

Observing with instruments such as the Atacama Large Millimeter/submillimeter Array (ALMA) is crucial to our understanding of planet-formation mechanisms, as we can observe at wavelengths that trace continuum emission from the

cold midplane (e.g., Testi et al. 2014), where we expect planets to be forming or have already formed. In the case of midplane spiral structures, their origin may be linked to the presence of a companion—stellar, fly-by, or planetary (Pohl et al. 2015; Bae & Zhu 2018a; Dong et al. 2018; Forgan et al. 2018b; Cuello et al. 2019; Keppler et al. 2020). Spirals may also be excited if the system is gravitationally unstable. Gravitationally instability is expected in cool and massive disks, where the disk-to-star mass ratio is larger than 0.1 (Bell et al. 1997; Gammie 2001; Lodato & Rice 2004; Cossins et al. 2009; Hall et al. 2016; Kratter & Lodato 2016; Rice 2016; Hall et al. 2019; Zhang & Zhu 2020). To date, not many spirals in dust continuum emission have a clear origin, except for those in multiple systems where the presence of spirals has been linked to stellar interactions (Kurtovic et al. 2018; Rosotti et al. 2020). On the other hand, there are disks where spirals have been reported at millimeter wavelengths and where no companion to which the

<footnote n="20">²⁰ NASA Hubble Fellowship Program Sagan Fellow.</footnote>
spiral origin may be linked has been detected yet (to date these
are Elias 27, IM Lup, WaOph 6, and MWC 758; Pérez et al.
2016; Dong et al. 2018; Huang et al. 2018c). If no companion
is detected and the disk is massive compared to the host star
mass, the gravitational instability (GI) scenario arises as a
possible explanation for the origin of the observed spirals.
Studying disks undergoing GI is important, as population
synthesis models show that GI primarily ends up forming
brown dwarf mass objects (Hall et al. 2017; Forgan et al.
2018a). It seems that giant-planet formation through GI is rare
(Rice et al. 2015), but it may still be the formation mechanism
for important systems like HR 8799 (Vigan et al. 2017).

Elias 2–27 is a young (0.8 Myr) M0 star (Andrews et al. 2009) located at a distance of $116_{-10}^{+19}$ pc (Gaia Collaboration et al. 2018) in the ρ Oph star-forming region (Luhman & Rieke 1999). It harbors an unusually massive protoplanetary disk; the disk-to-star mass ratio of Elias 2–27 is reported to be ~0.3 (Andrews et al. 2009; Pérez et al. 2016). The initial detection of two large-scale spiral arms was obtained with medium-resolution ALMA observations by Pérez et al. (2016). Due to the brightness and accessibility of the source, it became one of the Disk Substructures at High Angular Resolution Project (DSHARP; Andrews et al. 2018) targets, allowing further analysis of the dust emission at high resolution. Its distinctive morphology consists of two extended quasi-symmetric spiral arms and a gap, 14 au wide, located 69 au from the star (Huang et al. 2018b, 2018c). Due to its characteristic structure, the system has been the subject of several theoretical studies, concluding that GI is a possible origin of the spiral arms (Meru et al. 2017; Bae & Zhu 2018b; Forgan et al. 2018b; Hall et al. 2018). Though GI seems to explain the spiral morphology, it does not explain the dust gap, which could be carved by a companion of ~0.1 $M_J$ as constrained in hydrodynamical simulations by Zhang et al. (2018). Localized deviations from Keplerian motions at the location of this dust gap have been recently found, strengthening the hypothesis of a planetary-mass companion in the gap (Pinte et al. 2020). However, a lower-mass inner companion, such as the one proposed to open the gap, would not be able to excite the observed spiral arms (Meru et al. 2017).

Overall, Elias 2–27 appears to be a strong candidate to be a gravitationally unstable protoplanetary disk, but there are many tests to be done in order to determine if this is in fact the origin of the observed spirals. GI spirals will create pressure enhancements where we expect solids to be trapped and grain growth favored (Rice et al. 2004; Dipierro et al. 2015). This will not occur in the case of a companion, as companion-induced spirals will corotate with the planet at its Keplerian speed, faster than the background gas flow at their location, prohibiting dust growth and accumulation (Juhász et al. 2015). Another morphological signature is the expected symmetry of the spirals produced by GI, which should have a constant pitch angle in a logarithmic spiral model (Forgan et al. 2018b). Thus, measurements of dust-growth signatures together with symmetric, constant pitch angle, logarithmic spirals in a protoplanetary disk point toward a GI scenario.

Additionally, valuable dynamical information may be
obtained from gas observations. The presence of planets or
companions leaves distinct footprints in the kinematics and these
perturbations may be constrained by the amplitude of the gas
deviations from the expected Keplerian motion of an unper-
turbed disk. The current state-of-the art methods vary from

tracing pressure gradients (Teague et al. 2018), observing
deviations from expected isovelocity curves in the channel maps
("kinks"; Perez et al. 2015; Pinte et al. 2018a, 2019, 2020), and
using the mean velocity maps to model the velocity structure of
the disk and detect Doppler flips in the residuals (Pérez et al.
2018b, 2020). For GI spirals, Hall et al. (2020) characterize the
presence of a "GI-wiggle" that, contrary to companion-disk
interactions, will not be spatially localized but rather be a large-
scale perturbation, present throughout a wider velocity range,
colocated with the spirals. Analyzing the disk kinematics
complements the analysis of the observed dust structures and
allows us to understand and connect the various ongoing
processes.

Previously published gas observations of Elias 2–27 in the $^{12}$CO and $^{13}$CO in $J = 2–1$ transition show heavy absorption, as the star is quite embedded in its cloud (Andrews et al. 2009, 2018; Pérez et al. 2016; Pinte et al. 2020). In this study, we present $^{13}$CO and $^{C^{18}O}$observations in $J = 3–2$ transition. The higher energy transition and lower abundance of the isotopologs allow us to avoid some of the cloud contamination, while also probing closer to the midplane than previous works.

The present paper offers new observational constraints on
Elias 2–27 and is organized as follows. Section 2 provides an
overview on the calibration and imaging process of the
observations, Section 3 analyzes the spirals in the multi-
wavelength continuum data, Section 4 studies the $^{13}$CO and
$^{C^{18}O}$emission through a geometrical analysis of the moment
maps and the localization of perturbations in channel maps,
Section 5 shows the analysis of hydrodynamical simulations
computed for a GI disk using the derived observational
parameters of Elias 2–27, Section 6 discusses the results and
determines the possible origin of the spirals, and finally,
Section 7 summarizes the main findings of this work.

## 2. Observations

We present multiwavelength (Bands 3, 6, and 7) dust continuum and spectral line ($^{13}$CO $J=3-2$ and $^{C^{18}O}$ $J=3-2$) ALMA data of Elias 2–27. In the case of the Band 6 (1.3 mm) observations, the imaged data corresponds to the one presented in Pérez et al. (2016); detailed information regarding the calibration of this data set may be found in the original publication. For the Band 7 (0.89 mm) and Band 3 (3.3 mm) observations, the dust continuum data were first calibrated through the ALMA pipeline and afterward, phase and amplitude self-calibration was applied. Additionally, following the calibration procedure described in Andrews et al. (2018), we applied astrometric and flux scale alignment to correct for spatial offsets between the center of emission of different observations and relative flux scale differences. In both bands, we have short (15 m–313.7 m in Band 7, 15.1 m–2.5 km in Band 3) and long (15.1 m–1.4 km in Band 7, 21 m–3.6 km in Band 3) baseline observations in order to account for all of the emission from different spatial scales of the disk. Spatial offsets are corrected by locating the center of emission through a Gaussian fit in the image plane of each observation and adjusting phase and pointing centers with the Common Astronomy Software Applications (CASA; McMullin et al. 2007) tasks FIXVIS and FIXPLANETS, respectively. The absolute flux calibration uncertainty of ALMA data is expected to be ~10%, and we note that after self-calibration, the relative flux scales are consistent within ~5% for different observations in the same band. To adjust these amplitude scale differences, we compare the deprojected,
azimuthally averaged, visibility profiles and scale them using the short-baseline data as reference for the GAINCAL task, following the methodology in Andrews et al. (2018). For Band 7 data, we apply the spatial offset corrections and amplitude scaling before any self-calibration as relative flux scales vary 5%–10% between observation sets. In the Band 3 data set, we have a higher flux scale difference from the initial data sets (~20%), probably due to atmospheric conditions as we can visually see that the image quality is not the best. Therefore, in Band 3, we conduct self-cal of the short-baseline observation, and afterwards, with a flux scale difference of ~3%, we perform the spatial and amplitude scale corrections.

Self-calibration was first conducted on the dust continuum emission short-baseline data and the result combined with the long-baseline data to be self-calibrated all together. The selected time intervals for both phase and amplitude self-calibration start with the longest available total integration period for the initial round, and afterward, each following time interval was half of the previous interval. The longest and shortest intervals used were 1500 s and 46 s for the phase calibration of the joint 0.89 mm data, and 2860 s and 178 s for the phase calibration of the joint 3.3 mm data. For the calibration of the short-baseline data, the time intervals of phase calibration were between 1080 s and 270 s for the 0.89 mm data, and 304 s and 9.5 s for the 3.3 mm data. In the case of the amplitude self-cal time intervals, in Band 7, only one round was required in both short-baseline and joint data, so the solutions were obtained from the longest time interval. In Band 3, the short-baseline data had one round of amplitude calibration and the joint data set had two rounds. Self-calibration rounds were applied until signal-to-noise ratio (S/N) improvement was less than 5% in the case of Band 7 data and until there was no S/N improvement in the case of Band 3. For the short-baseline data at 0.89 mm, we conducted three rounds of phase calibration and one round of amplitude calibration. The combined data set at 0.89 mm had seven rounds of phase calibration and one round of amplitude calibration; overall, the peak S/N in the joint data set improved 300%. The short-baseline data at 3.3 mm had six rounds of phase calibration and one round of amplitude calibration; the joint data set had five rounds of phase calibration and two rounds of amplitude calibration. In the joint data set, the S/N improvement was 6%.

For the final imaging, we used multiscale TCLEAN for all the images, using a $1\sigma$ stopping threshold (where $\sigma$ is the image rms) for the Bands 3 and 6 data, and a $2\sigma$ stopping threshold in Band 7. Robust weighting values were 1.0 in the case of Band 6, and 0.5 for Bands 3 and 7. We also set the GAIN value to 0.05 and CYCLEFACTOR parameter in TCLEAN to 2.0, to have a more detailed cleaning, by substracting a smaller fraction of the source flux from the residual image and triggering a major cycle sooner. All images have comparable beam sizes: $0''22 \times 0''17$ beam (26 au × 20 au) at 0.89 mm, $0''26 \times 0''22$ beam (30 au × 26 au) at 1.3 mm, and $0''26 \times 0''20$ (30 au × 23 au) beam at 3.3 mm. At each wavelength, the image rms is approximately 93 µJy/beam, 86 µJy/beam, and 10 µJy/beam, respectively, for 0.89 mm, 1.3 mm, and 3.3 mm.

The observations of C$^{18}$O and $^{13}$CO were obtained simultaneously with the dust continuum emission in Band 7 (0.89 mm) and in the $J=3-2$ transition. C$^{18}$O is observed at 329.330 GHz, with a spectral resolution of 0.035 MHz, and $^{13}$CO at 330.588 GHz, with a spectral resolution of 0.121 MHz. After

applying the same self-calibration solutions as to the dust continuum of Band 7, the emission was first imaged using natural weighting (robust parameter of 2.0) and not applying any uv-tapering or uv-range filtering, in order to be sensitive to large-scale emission with a beam size of $0''29 \times 0''22$ (~34 × 25 au) for $^{13}$CO and $0''30 \times 0''23$ (~35 × 27 au) for C$^{18}$O (see figures in Appendix A). In order to avoid the cloud contamination present in the disk (Pérez et al. 2016; Huang et al. 2018c), we exclude large-scale emission by considering only baselines longer than 36 m for $^{13}$CO (scales shorter than $6''4$, ~740 au) and 45 m for C$^{18}$O (scales shorter than $5''11$, ~590 au). We also applied uv-tapering of $0''2 \times 0''115$, PA = 0° for $^{13}$CO and $0''2 \times 0''0$, PA = 0° for C$^{18}$O in order to obtain a roughly round beam. The final images that trace the disk emission were obtained with a robust parameter of 0.5, resulting in a beam size of $0''26 \times 0''25$ (~29 au) for $^{13}$CO and $0''31 \times 0''29$ (~35 au) for C$^{18}$O. We imaged all channel maps in both $^{13}$CO and C$^{18}$O using a 0.111 km s$^{-1}$ spectral resolution, though the C$^{18}$O data were observed with a finer spectral resolution ($0.032$ km s$^{-1}$) and found the best compromise of S/N by using a broader spectral resolution.

We also recover gas emission from CN $v=0, N=3-2$ in $J=7/2-5/2$ and $J=5/2-3/2$; these additional molecules will be studied separately in future work. We do not recover any emission from requested observations of CN $v=0, N=3-2$ in $J=5/2-5/2$ or SO $3\Sigma v=0 J=3-2$, at the achieved sensitivity level of 3 mJy/beam.

## 3. Dust Spiral Structure

We recover the spiral structure at all wavelengths, as shown in Figure 1. Given the ~25 au beam size, we are not able to fully resolve the 69 au gap but can distinguish it as a small decrease in brightness temperature at all wavelengths in Figure 1. For all calculations throughout this study, we will assume the gap location, disk inclination, and disk position angle as derived in Huang et al. (2018b), which are $69.1 \pm 0.4$ au, $56.2 \pm 0.8$, and $118.8 \pm 0.7$, respectively.

### 3.1. Tracing the Spiral Morphology

To trace the spiral structure from our dust continuum images, we radially subtract an azimuthally averaged radial profile of the emission and trace the spiral features from these “subtracted images” (as done in Huang et al. 2018c). From the subtracted images, shown in the top panels of Figure 2, we find the radial location of maximum emission along each spiral, at azimuthal steps sampled every 9°, between a determined azimuthal angle range (values in Table 1). The radial extent where we trace the maxima of emission is determined visually. To aid our visual criteria, we consider radial locations no further than where we have a signal of 5 times the rms in the nonsubtracted image.

Previous analyses of Elias 2–27 (Pérez et al. 2016; Huang et al. 2018c) have shown that a logarithmic spiral model, with a constant pitch angle, can adequately trace the spiral morphology. Therefore, we use Markov Chain Monte Carlo (MCMC) modeling to find the best-fit parameters for a logarithmic spiral model, considering the location data of each spiral and each observed wavelength. The spiral form is given by

$$r(\theta) = R_0 e^{b\theta}. \qquad (1)$$

Here, $\theta$ is the polar angle in radians measured from the north and to the east (left); for reference, see coordinates in Figure 1. $R_0$, $b$ are free parameters, with $R_0$ the radius at $\theta=0$, measured
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Dust continuum observations of Elias 2–27 at 0.89, 1.3, and 3.3 mm. For each panel: the intensity color scale is shown on the right, the scale bar in the lower-right corner corresponds to 30 au at the distance of the star, and the ellipse in the bottom-left corner indicates the spatial resolution.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** The spiral morphology of Elias 2–27 at multiple wavelengths. Panels from left to right correspond to data from the 0.89, 1.3, and 3.3 mm observations. The northwest spiral is traced in blue; the southeast spiral is traced in red. Top panels: dust continuum maps from which the azimuthally averaged radial profile has been subtracted to highlight spiral location, red and blue points trace the maxima of emission along the arms. Middle panels: deprojected radial location of the emission maxima, as a function of the angle measured from north to the east (left). Error bars correspond to the astrometric error of each data point and gray lines show the posterior distribution of a logarithmic spiral fit with a constant pitch angle. Bottom panels: deprojection of the subtracted dust continuum observations from the top panels. The vertical line marks the dust gap location from Huang et al. (2018c), and colored lines show the best-fit logarithmic spiral model.</figcaption></figure>
<table><caption>**Table 1**

Best-fit Parameters of the Logarithmic Spiral Model</caption><thead><tr><th>Wavelength</th><th>Spiral Arm</th><th>Angle Range</th><th>R<sub>0</sub> (au)</th><th>b</th><th>Pitch Angle</th></tr></thead><tbody><tr><td rowspan="2">0.89 mm</td><td>NW</td><td>-251° to -5°</td><td>249.8<sup>+1.2</sup><sub>-1.1</sub></td><td>0.230 ± 0.002</td><td>12°9 ± 0°1</td></tr><tr><td>SE</td><td>-70° to 150°</td><td>111.5 ± 0.4</td><td>0.247 ± 0.003</td><td>13°9 ± 0°2</td></tr><tr><td rowspan="2">1.3 mm</td><td>NW</td><td>-250° to -35°</td><td>250.2<sup>+2.3</sup><sub>-2.5</sub></td><td>0.229<sup>+0.003</sup><sub>-0.004</sub></td><td>12°9 ± 0°2</td></tr><tr><td>SE</td><td>-61° to 135°</td><td>115.5<sup>+0.8</sup><sub>-0.9</sub></td><td>0.234 ± 0.007</td><td>13°2 ± 0°4</td></tr><tr><td rowspan="2">3.3 mm</td><td>NW</td><td>-260° to -23°</td><td>249.9<sup>+3.9</sup><sub>-3.8</sub></td><td>0.229 ± 0.005</td><td>12°9 ± 0°3</td></tr><tr><td>SE</td><td>-50° to 140°</td><td>113.6<sup>+1.0</sup><sub>-1.1</sub></td><td>0.231 ± 0.010</td><td>13°0 ± 0°6</td></tr></tbody></table>

in astronomical units, and *b* relates to the pitch angle of the spiral arms (*φ*), as *φ* = arctan(*b*). The uncertainty on the location of each measured maxima is assumed to be the astrometric error<fnref n="21" />: Δ*p* = 60,000 · (ν · *B* · S/N)⁻¹, where Δ*p* is the approximate position uncertainty of a feature in milliarcseconds, S/N is the peak/rms intensity ratio of the data point on the image, ν is the observing frequency in GHz, and *B* is the maximum baseline length in kilometers.

The top panel of Figure 2 shows the maxima along each spiral from the subtracted images. In the middle panels are the deprojected radial locations of each maximum, measured from the center, as a function of azimuthal angle. Gray lines show logarithmic spiral models, derived from 300 draws of the posterior values after convergence of the MCMC simulations for the best-fit parameters of each spiral arm. The horizontal line at 148 au marks the location where we observe a break in the spiral arm, at all wavelengths, clearest in the southeast spiral, also subtly present in the northwest spiral. Huang et al. (2018c) had previously noted a possible decrease in the pitch angle value outside *R* ~ 150 au. The bottom panel (Figure 2) shows the polar deprojection of the subtracted data, with a vertical line marking the dust gap location. The red/blue lines show the best-fit model for the logarithmic spirals. The median values for the parameters of the logarithmic spiral model, for each spiral and wavelength, are shown in Table 1, along with the 16th and 84th percentile uncertainties derived from the posteriors.

We note that the pitch angle values retrieved here, ~12°9 and ~13°3, for the NW and SE spiral, respectively (similar between wavelengths; see Table 1), are different from the recovered pitch angles from Huang et al. (2018c; 15°7 and 16°4 for the NW and SE spirals). When applying our method to the high-resolution data set presented in Huang et al. (2018c), we retrieve the same results as in their work. The pitch angle difference between this work and theirs is probably due to beam-smearing effects, combined with the challenge of image subtraction in lower-resolution data (as discussed in Huang et al. 2018c), as the angular resolution difference between data sets is a factor of 4–5.

## 3.2. Contrast Variations along the Spirals

We test the contrast at comparable locations over different wavelengths, for this may provide evidence of dust trapping. If indeed there is dust growth within the spiral arm or large particles are trapped at this location, then at longer wavelengths, we could expect to observe a higher contrast. The latter effect is due to larger particles being more densely packed in the spiral arm (Rice et al. 2004); this has been shown to occur

in the case of dust-trapping vortex studies (Cazzolletti et al. 2018). We compute how the contrast varies radially along each spiral arm, with respect to a fixed “interarm” region, which is expected to be a location where we can be sure there is no emission from the spiral arms. Considering the symmetry between the spirals, the interarm region of each spiral is assumed to follow the same shape of the corresponding spiral arm but rotated 90° clockwise from the spiral location. The top panels of Figure 3 show a polar projection of the dust continuum, overlaid are the best-fit model of the spirals for each wavelength in a continuous line and this same model, shifted 90° clockwise, tracing the interarm region, in dashes. From this figure, it is clear that the interarm region effectively traces zones of lower emission, compared to the spiral arm location.

The contrast is calculated, from the original images, as the ratio between the emission at the spiral and the interarm regions, at each radius, and for both spiral arms. We use the previously derived parameters and calculate the contrast along the best-fit spiral model for each wavelength, binning the data azimuthally every 15°. The standard deviation within the bin is used as the error of each averaged flux measurement, both for the interarm and spiral region. The contrast curves for each wavelength are shown in the bottom panels of Figure 3. The NW spiral remains the strongest spiral at all wavelengths, with an average contrast value 25%–33% higher than the SE spiral. Both spirals maintain a similar contrast throughout their radial extent, deviating in the 0.89 and 1.3 mm observations only ~11% from the average contrast value in the case of the SE spiral and ~17% in the case of the NW spiral. In the 3.3 mm data, the contrast variation increases and deviations from the average value increase to ~16% in the SE spiral and ~24% in the NW. For the 0.89 mm, 1.3 mm, and 3.3 mm emission, the radial distance of the maximum contrast value is, respectively, for (NW, SE) spirals, (~125 au, ~117 au), (~129 au, ~127 au), and (~129 au, ~126 au). For the minimum contrast location, in the same order, the radial distances are (~155 au, ~149 au), (~161 au, ~147 au), and (~150 au, ~142 au). We observe in our contrast curves that the radial location of the maximum contrast slightly shifts outwards with wavelength, no global shifts are detected in the minimum contrast distances; however, the minimum contrast of the SE spiral appears to shift inwards at a longer wavelength. Overall, our values are in agreement with the previous measurements of Huang et al. (2018c), where the peak contrast is located at ~123 au and low contrast at ~147 au (measured from higher angular resolution 1.3 mm data). While the form of the contrast curve in the 1.3 mm emission agrees also with the previous studies (Pérez et al. 2016; Huang et al. 2018c), at larger radial distances, we observe a decrease in the contrast of both spirals in the 0.89 mm emission and also in the SE spiral in 3.3 mm, contrary to the apparent

<footnote n="21">²¹ See Section 10.6.6 in https://almascience.nrao.edu/documents-and-tools/cycle5/alma-technical-handbook/view for further details.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Contrast of the spiral arms in Elias 2–27, panels from left to right correspond to data from the 0.89, 1.3, and 3.3 mm observations. Top panels: polar deprojection of the dust continuum emission maps, where the azimuthal angle is measured from north and to the east. Blue and red colors correspond to the northwest and southeast spirals, respectively. Continuous lines trace the spirals following the best-fit parameters found in Section 3.1, dashed lines trace the interarm region following the same best-fit parametric model, but minus 90° from the spiral. Bottom panels: calculated contrast values along each spiral arm; blue and red points correspond to the northwest and southeast spirals, respectively.</figcaption></figure>

growing contrast at larger radial distance observed in the 1.3 mm emission.

If we follow the peak contrast of both spirals, located at ~125 au, we see a slight variation in the contrast value, which grows larger with longer wavelength. The peak contrast value of the NW spiral increases by 11% from the 0.89 mm to the 1.3 mm emission and by 8% between the 1.3 mm and 3.3 mm emission. The SE spiral increases its peak contrast value by 5% initially and in the 3.3 mm data decreases the peak contrast by 2% with respect to the 1.3 mm emission. Minimum contrast values also shift but do not show any correlation with varying wavelengths. We note that previous studies of this source measured the contrast differently, using the ratio between the spiral flux and the minimum flux at the same radial location (Pérez et al. 2016; Huang et al. 2018c). For completeness, we test our results using the minimum flux method used in the previous studies and obtain similar contrast curves in which we also measure increasing peak contrast values toward longer wavelengths. However, we do not consider these results, for the minimum flux at a determined radial distance is found at different azimuthal angles between different wavelengths. Therefore, the contrast calculated using the minimum radial flux does not sample comparable locations between different wavelengths, something necessary in order to detect possible

signatures of dust growth at a given location, which we achieve by using the interarm region.

### 3.3. Spectral Index Analysis

The presence of dust growth throughout the disk can be inferred using the multiwavelength continuum observations to compute the disk spectral index. The spectral index ($\alpha$) of the spectral energy distribution follows $I_\nu \propto \nu^\alpha$, where $I_\nu$ corresponds to the measured intensity, from the image plane, at frequency $\nu$. In the optically thin regime, if $\alpha$ has values between 2 and 2.5, it can be an indicator toward the presence of large, up to centimeter-sized, grains, and hence, dust growth at the location of such values (e.g., Kwon et al. 2009; Testi et al. 2014; Pinilla et al. 2017; Cazzoletti et al. 2018; Macías et al. 2019). To adequately compare the emissions from all three wavelengths, the calculations in this section are done from dust images at equal resolution and centered in the peak brightness. We use the task imsmooth in CASA to specify a round, 0.26 × 0.26 beam (~31 au). Figure 4 shows in the top panel the azimuthally averaged brightness temperature of the smoothed, equal resolution images.

We calculate $\alpha$ as the best-fit slope of a linear model applied to the logarithmic space ($\alpha \propto d \ln I_\nu/d \ln\nu$), using all wavelengths available in this study. From an MCMC fit we
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Radial profiles of the brightness temperature, spectral index, and optical depth for Elias 2–27 at multiple wavelengths. Top panel: azimuthally averaged brightness temperature profiles for each wavelength, the shaded area shows the $1\sigma$ scatter at each radial bin divided by the beams spanning the angles over which the intensities are measured. Middle panel: spectral index radial profile considering all three wavelengths, the shaded area shows $1\sigma$ errors from assuming the intensity error from the azimuthally averaged intensity profile. Lightly shaded area shows the total error including a 10% flux uncertainty on the intensity measurement. Lower panel: optical depth radial profile for each wavelength, the shaded area shows $1\sigma$ uncertainty considering the $1\sigma$ errors on the stellar luminosity and from the azimuthally averaged intensity profiles. The vertical dashed line and shaded region mark the dust gap location (69 au) and width (14.3 au) as determined in Huang et al. (2018b).</figcaption></figure>

obtain a posterior distribution of $\alpha$, with best-fit and uncertainties computed from the 50th, 16th, and 84th percentile values.

From the azimuthally averaged intensity profiles, we compute an $\alpha$ radial profile shown in the middle panel of Figure 4. We retrieve a disk average $\alpha$ value of $2.6 \pm 0.06$, with $\alpha < 2.0$ in the inner $\sim 50$ au. These low values in the inner disk will be discussed in Section 5.1 and are probably related to optically thick emission where self-scattering at long wavelengths is a relevant process.

We apply the same techniques described in Huang et al. (2018b) to calculate the optical depth at our various wavelengths. The main assumption relies on approximating the midplane

temperature profile ($T_{\text{mid}}$), assuming a passively heated, flared disk in radiative equilibrium (e.g., Chiang & Goldreich 1997), and considering that the millimeter emission follows this temperature profile if it traces the midplane. From the temperature profile, we compute the expected blackbody emission ($B_{\nu}$), using Planck's law, and relate it to the measured intensity ($I_{\nu}$) through $\tau_{\nu}$, following $I_{\nu}(r) = B_{\nu}(T_{\text{mid}}(r))(1 - \exp(-\tau_{\nu}(r)))$. $T_{\text{mid}}$ will depend on the assumed flaring of the disk and the stellar luminosity (we use the DSHARP values for flaring and stellar luminosity of $\phi = 0.02$ and $\log L_*/L_{\odot} = -0.04 \pm 0.23$; Andrews et al. 2018).

The computed optical depth profiles are shown in the bottom panel of Figure 4. We note that at 0.89 mm the modeled $T_{\text{mid}}$ is $\sim 2$ times larger than the measured $T_b$. In the Rayleigh–Jeans regime, as expected for millimeter emission, $T_b = T_{\text{mid}}(1 - \exp(-\tau))$. Then, if $T_{\text{mid}}$ is overestimated it will produce lower optical depth values and the emission will appear more optically thin than it really is. The disk could be colder, or with a lower flaring value, and the emission at all wavelengths would be more optically thick. These issues will be discussed in Section 5.1.

The spectral index map for Elias 2-27 and the uncertainty on the spectral index are presented in Figure 5. We only consider emission above $5\sigma$ at all wavelengths for this calculation, and we adopt the rms of each image as the uncertainty on the intensity for each pixel. The best-fit spiral model from the 0.89 mm dust emission is overlaid for reference, along with the dust gap location (69 au). Contour lines for $\alpha$ values of [2.0, 2.25, 2.4] and for $\alpha$ uncertainties at the level of [0.01, 0.025, 0.05, 0.1] are presented. We observe that contours tracing $\alpha = 2.4$ seem to be coherent with the spiral morphology as traced by the best-fit spiral. This is specially apparent in the NW spiral location.

The variations between the spiral and interarm region at the NW spiral location are $\sim 0.1$, with uncertainties of $\sim 0.04$. Overall, we observe lower $\alpha$ co-located with the spiral features. At the gap location no particular behavior is observed, except for a decrease in $\alpha$ values inwards from the gap location. As noticed in the radial profile, in the inner disk ($<50$ au) $\alpha$ reaches values below 2.0.

Another indicator of dust growth is that smaller grains, traced by shorter wavelengths, will be less concentrated than larger grains observed at longer wavelengths. This will translate into a width difference along the dust-trapping structure, such that smaller grains will be more widely spread than larger grains within a dust trap. Such behavior has been constrained in different sources with vortex-like structures, likely tracing dust traps (Casassus et al. 2015, 2019; van der Marel et al. 2015; Cazzoletti et al. 2018). The measurement of the width variation cannot be done using a subtracted image or an azimuthally averaged profile, both options would introduce artifacts or remove relevant emission, therefore, we must use the original image. We decide to trace variations along a single azimuthal cut, choosing the angle where the highest spiral contrast is observed, which is located at $\sim 125$ au corresponding to $\sim 224^{\circ}$.

Figure 6 shows the brightness temperature profiles along the azimuthal angle of maximum spiral contrast (224°) for all three wavelengths, the spirals can be localized as a bump in the intensity curve between 100 au–150 au in the top panel. From this it can be seen that the shape of the intensity curve is similar at all three wavelengths, with no noticeable width differences. The bottom panel shows the spectral index distribution along
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: spectral index map calculated from the emission of all available wavelengths. Only emission over $5\sigma$ in the image plane is considered. Blue lines show the location of the dust gap at 69 au and the derived best-fit parametric model for the spiral arms of 1.3 mm emission. White contours correspond to $\alpha = 2.0$ and 2.25, black contours indicate $\alpha = 2.4$. Right: spectral index error map, obtained from the image rms and intensity value of each pixel. Blue lines trace the dust features and black contours correspond to $\alpha$ uncertainties at 0.01, 0.025, 0.05, and 0.1 levels.</figcaption></figure>

the azimuthal cut, calculated directly from the top panel intensity profile, as previously done for the spectral index map. In this $\alpha$ profile we see that the east side (negative radial distance) shows a small decrease in the spectral index value at the spiral location, while the west spiral does not show significant $\alpha$ variations. The spatial extent of the variation in the east is resolved by our beam size (0.26, ~31 au).

## 4. $^{13}$CO and C$^{18}$O J = 3–2 Emission Analysis

We analyze the emission from the two observed CO isotopologs in three different ways by (1) studying the presence of structures in the integrated emission and in channel maps, including large-scale emission in the entire field of view (FOV), (2) searching for velocity perturbations in channel and velocity maps, (3) constraining the vertical height of the $^{13}$CO- and C$^{18}$O-emitting layer, and analyzing the kinematics of these isotopologs.

### 4.1. Channel and Moment Maps

The channel maps for both $^{13}$CO and C$^{18}$O are shown in Appendix A. The known cloud absorption (Pérez et al. 2016; Huang et al. 2018c) affects the east side of the disk and blocks all $^{13}$CO emission from $v \ge 2.88$ km s$^{-1}$, while for C$^{18}$O some cloud absorption is present near $v = 2.55–2.66$ km s$^{-1}$. We note that the south part of the disk is the brightest. The systemic velocity is determined to be 1.95 km s$^{-1}$ based on the C$^{18}$O channel map analysis using a spectral resolution of 0.05 km s$^{-1}$ and kinematic modeling done with the EDDY package (Teague 2019).

Considering all spatial scales and imaging the whole FOV, extended, large-scale emission appears in both isotopologs around Elias 2–27 (see figures in Appendix A). As $^{13}$CO is the most abundant isotopologue, the large-scale emission appears more strongly and along a wider range of velocities than in C$^{18}$O. The extended emission is clearly identified between $v = 2.55–4.88$ km s$^{-1}$ for $^{13}$CO, and between $v = 2.66–4.21$ km s$^{-1}$ for C$^{18}$O. This large-scale emission appears to have a

striping pattern, probably due to the lack of compact baselines in the observations. The shortest baseline is 15 m, which, projected in the sky for the observed frequencies, recovers emission from $\lesssim 15''$ scales; however, the detected emission may extend beyond our $\sim 20''$ FOV. We compute integrated emission and mean velocity maps, shown in Figure 7, from the channel maps that consider the whole FOV, including emission above $3\sigma$, while in $^{13}$CO, the large-scale emission is seen through the entire FOV, the C$^{18}$O large-scale emission is more constrained and crossing only through the east side of the disk. There is no clear velocity gradient between the large-scale emission and the disk in either isotopolog.

From here on, we focus on the channel maps that trace the material closer to the disk rather than large-scale emission. These channel maps were obtained using uv-tapering and filtering of the short baselines as described in Section 2. Figure 8 shows the integrated emission (moment 0) of both CO molecules. Foreground absorption in $^{13}$CO is clear on the east side of the disk. Overlaying the continuum contours for the 0.89 mm emission and measuring along the major axis, the east side of the disk appears to have a larger extent than the west side in the gas. Figure 8 shows that the size difference is roughly 0.5 between the east and west sides of the disk, in both CO tracers. At the source distance, this size difference corresponds to 58 au between the projected emission extent of each side. To measure the east/west size variations, we trace the edge of the emission in the C$^{18}$O moment 0 map. The border of the disk is considered as the radius that encompasses 98% of the integrated emission at each sampled azimuthal angle (green points on the top panel of Figure 9). We define the center based on the emission peak from the Band 7 continuum data and deproject accordingly, assuming the inclination and position angle of the dust continuum. The deprojected radial distance as a function of azimuthal angle is shown in the middle panel of Figure 9. Errors correspond to the astrometric error, calculated as described in Section 3.1. The edge of the disk is not well described by a circle or an ellipse; it shows two local maxima and two local minima. The global minimum distance is located along the major axis on the west,
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Top: brightness temperature profile for the 0.89, 1.3, and 3 mm emission along the azimuthal angle of maximum spiral contrast (224°), measured radially. Positive radial values are measured from the center of the disk and to the west; negative values indicate the distance to the east. Vertical colored dashed lines show the location of the spiral arm according to the best-fit parametric model at each wavelength, black vertical line marks the location of the dust gap. Bottom: spectral index along the azimuthal cut, calculated from the emission of the three wavelengths.</figcaption></figure>

but the global maximum distance is shifted with respect to the east major axis. The locations of maximum and minimum border extents are not colocated with the continuum spiral features or their extension.

Another feature in the moment 0 map is the presence of an emission “gap” at large disk radii in C¹⁸O. To estimate its location, we trace the radial positions of emission minima sampling every 9°, between 185 and 300 au (this range is determined by analysis of the intensity profile). Using the mean value and standard deviation of the minima radial locations (yellow points in the top panel of Figure 9), we estimate the gap position at 241 ± 24 au. In the ¹³CO-integrated intensity map, we do not observe a gap and cannot infer one from the intensity profile, even when excluding the azimuthal angles ~65°–155°, where foreground absorption is strongest. If we trace the emission border of ¹³CO following the 98% integrated emission criteria, we obtain a similar emission border curve to the C¹⁸O, shown in Appendix A.

## 4.2. Tracing the Emitting Layer in ¹³CO and C¹⁸O

Using the method detailed in Pinte et al. (2018b), we recover the emission surface of each molecule. This is done by tracing the maxima from the upper layer of emission in the channel maps and applying geometrical relationships. The emitting layer is assumed to have a cone-like structure; this means that the height of the emitting layer should be symmetric with respect to the disk’s major axis.

If the variations in the projected radial extent of the disk found in Section 4.1 are related to variations in the emitting surface height, then we would only expect symmetry along the disk’s semimajor axis on the west side of the disk. The latter can be determined from Figure 9, where we see that at the west major axis location (~-60°) the radial distance of the emission border grows similarly when we move toward the north (positive angles) or south (negative angles). This similarity in the radial distance increment is maintained for ~90° in each direction (north and south), which corresponds to the complete west side of the disk. On the other hand, at the location of the east major axis, we do not see this symmetry in the growth or decline of the radial distances toward the north or south. The possible lack of symmetry is an important caveat, and we

expect this to affect mostly the constraints obtained for the east side of the disk (where the variations in disk extension are larger). The emission surface we measure should be taken as a rough estimate. Details on the geometry relations can be found in the original publication (Pinte et al. 2018b).

From all available channels, we visually select those in which the top layer of gas emission can be clearly identified. For C¹⁸O, we select channels at velocities +0.77 to +1.55 km s⁻¹ and +2.44 to +3.21 km s⁻¹. For ¹³CO, we select channels +0.66 to +1.66 km s⁻¹ and +2.33 to +2.77 km s⁻¹.

The recovered height profile of the emission layer for both gas isotopologs is shown in Figure 10. Measurements obtained from channels that trace the east sides of the disk (with respect to the semiminor axis) are colored red, while those from the west are blue. The ¹³CO and C¹⁸O gas emission layer we constrain follows the expected distribution for these isotopologs in a disk: ¹³CO traces a higher layer from the midplane than C¹⁸O in both east and west sides at all radii. The continuum and C¹⁸O gas gap location are highlighted in Figure 10, together with the width for the dust gap (as reported in Huang et al. 2018b) and the uncertainty regions of each gap. No clear feature is recovered in the height profiles at the gap locations.

The height distribution is modeled using the equations for a complex flared surface presented in the EDDY package (Teague 2019), such that the altitude of the emitting layer follows a double power law:

$$z(r) = z_0 \left(\frac{r}{r_0}\right)^\psi + z_1 \left(\frac{r}{r_0}\right)^\varphi, \quad (2)$$

where r is the radial distance from the star, as measured in the azimuth plane, and the characteristic radius $r_0$ is fixed at 1", corresponding to 116 au for our system. The second power law is aimed to be a correction on the first term; therefore, we first find the best set of parameters for the emission layer characterized with a single power law and then optimize close to those parameters, to include the second power law. The best-fit parameters of the model are assumed to be the median value from the posteriors of the MCMC simulations and are shown with their uncertainties (16th and 84th percentile uncertainties derived from the posteriors) in Table 2.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The left column shows $^{13}$CO emission maps, the right row $^{C^{18}O}$emission maps. Top row: integrated emission maps considering all emission over $3\sigma$. Lower row: mean velocity maps.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Integrated emission (moment 0) maps for $^{13}$CO (left) and $^{C^{18}O}$ (right) gas emission. Contours of 0.89 mm continuum emission are overlaid on top. The white grid marks the minor and major axes of the disk, as determined by the continuum emission position angle, with ticks on these axes indicate $0.′′5$ (~58 au) intervals.</figcaption></figure>

Our modeled emission surface presents consistent differences in the elevation and morphology between the east and west sides of the disk. The height profiles have a quasilinear form in the east channels, mostly tracing lower height values

than the west. At larger radii ($>200$ au), the west channels show a decrease in the emission surface height. Residuals obtained from subtracting the model from the observations are shown in the bottom panels of Figure 10. We see that the
### 4.3. Tracing the Kinematics in $^{13}$CO and C$^{18}$O

Besides tracing the height profile of the emitting layer, Pinte et al.'s method allows us to determine the velocity profile of the traced emission layer. In a given velocity channel, we know the projected radial velocity, $v_{\text{obs}}$, together with the systemic velocity for the source, $v_{\text{syst}}$. We are interested in determining the azimuthal velocity $v$ of a parcel of gas at an azimuthal radial distance $r$ and height $h$ from the star. Using the inclination angle $i$ and the polar azimuthal angle, $\theta$, we can relate the known velocities to $v$ through $v_{\text{obs}} = v_{\text{syst}} + v \cos(\theta)\sin(i)$. To obtain $\cos(\theta)$ we apply geometrical relationships from the measurements in the channel maps as defined in Pinte et al. (2018b). The velocity profile will allow us to obtain a mass estimate for the central star and also test for super-Keplerian velocities at large radial distances, which is a characteristic expected in disks undergoing GI (Bertin & Lodato 1999; Lodato 2007).

The velocity profiles of the emitting gas for both $^{13}$CO and C$^{18}$O are shown in Figure 11. An overall velocity difference is observed between both sides of the disk, with measurements from the east side having a higher velocity. From our previous results (Figure 10), we know that the east side corresponds to the side apparently closest to the midplane. The difference in the velocity profile is in agreement with the height profile variations between sides, as being closer to the midplane results in larger velocities. We do not observe any noticeable behavior of the velocity profile at the location of the dust and gas gaps.

We fit the velocity profiles with a Keplerian model to constrain the mass of the central star. Based on comparisons to stellar evolution models in the Hertzsprung-Russell diagram, the mass of Elias 2-27 has been reported to be ~0.49$M_{\odot}$ (Andrews et al. 2009, 2018; Pérez et al. 2016). For our modeling, we incorporate the height distribution of each molecule, using the best-fit double power-law model found previously (see Table 2). The modeled velocity at a radial distance $r$ from the star will follow Equation (3), where $G$ is the gravitational constant, $M_{*}$ is the central star mass, and $h$ is the height of the gas at azimuthal radius $r$:

$$ \frac{v^2}{r} = \frac{GM_*r}{(r^2 + h^2)^{\frac{3}{2}}} \qquad (3) $$

We note that Equation (3) does not include the effects of the radial pressure gradient and the disk's self-gravity (Rosenfeld et al. 2013).

We simultaneously fit the model to the data points from the east and west sides of the disk, using MCMC simulations and taking into account the different height profiles (Figure 10). The curves for the expected Keplerian motion, considering the best-fit stellar mass and its $1\sigma$ uncertainty range, are shown over the data in Figure 11. The final masses and errors are computed from the median value and 16th and 84th percentile uncertainties derived from the posteriors. From the $^{13}$CO measurements, we constrain a stellar mass of $M_{*} = 0.5 \pm 0.01 M_{\odot}$, while from the C$^{18}$O measurements, we constrain $M_{*} = 0.46_{-0.03}^{+0.02} M_{\odot}$. Both values are similar between them and compared to the previous estimates ($M_{*} \sim 0.49 M_{\odot}$; Andrews et al. 2009, 2018; Pérez et al. 2016).

Compared to the expected Keplerian velocity profile from the fits of Figure 11, we see residuals throughout the whole radial extent. These simultaneous sub- and super-Keplerian velocities are expected if the emission layer height difference between the east and west sides was larger than what was

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Top panel: integrated emission map of the C$^{18}$O gas emission. Yellow dots trace the local minimum of emission, green circles trace the disk radial extent as the location at which 98% of the total azimuthal emission is included. Middle panel: radial distance from the center of the points tracing the border, the azimuthal angle is measured from the north to the east. Bottom panel: azimuthally averaged intensity profile of the C$^{18}$O emission, the shaded area shows the $1\sigma$ scatter at each radial bin divided by the number of beams spanning the angles over which the intensities are measured. The vertical yellow line marks the average radial location of the minimum of emission and the deviation of the data is indicated by the vertical gray region.</figcaption></figure>

largest residual scatter is found between the dust and gas gap locations, which roughly coincides with the radial extent of the dust spiral arms (80–250 au). The residuals from the $^{13}$CO emission show a “curved” pattern, indicating a more complex emitting surface. We note that in both isotopologs, the residuals at the dust gap location are mostly negative, indicating a possible decrease in the emitting surface at this radii. This is unresolved with our spatial resolution.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Top panels show the emission layer height as a function of radial distance to the star constrained from the ¹³CO (left) and C¹⁸O (right) data. Blue points correspond to measurements coming from the west side of the disk, the red points come from the east side, the colored line corresponds to the best-fit double power-law height profile for the data, and the gray area shows the uncertainty on the emitting layer as derived from the posteriors. The vertical dashed line indicates the location of the gap reported in the continuum, the dotted-dashed line corresponds to the gap location in the C¹⁸O-integrated intensity map. The gray area indicates the width of the dust gap (obtained from Huang et al. 2018b), the orange area indicates the gas gap's location uncertainty. The bottom panels show the residuals of each isotopolog after subtracting the best-fit model to the data.</figcaption></figure>

<table><caption>**Table 2**

Height Model Parameters from Channel Analysis</caption><thead><tr><th>Param.</th><th>¹³CO - W</th><th>¹³CO - E</th><th>C¹⁸O - W</th><th>C¹⁸O - E</th></tr></thead><tbody><tr><td>z₀ (au)</td><td>77.6<sup>+8.2</sup><sub>-8.9</sub></td><td>30.0 ± 0.2</td><td>44.0<sup>+7.3</sup><sub>-6.2</sub></td><td>32.3<sup>+8.9</sup><sub>-8.5</sub></td></tr><tr><td>ψ</td><td>1.76<sup>+0.03</sup><sub>-0.04</sub></td><td>1.01 ± 0.01</td><td>2.19<sup>+0.12</sup><sub>-0.16</sub></td><td>1.14<sup>+0.11</sup><sub>-0.09</sub></td></tr><tr><td>z₁ (au)</td><td>-38.1<sup>+8.9</sup><sub>-8.1</sub></td><td>-0.3<sup>+0.1</sup><sub>-0.2</sub></td><td>-11.8<sup>+6.1</sup><sub>-7.4</sub></td><td>-15.9<sup>+8.6</sup><sub>-8.9</sub></td></tr><tr><td>φ</td><td>2.29<sup>+0.08</sup><sub>-0.04</sub></td><td>-3.36<sup>+0.54</sup><sub>-0.52</sub></td><td>3.58<sup>+0.46</sup><sub>-0.24</sub></td><td>0.46<sup>+0.15</sup><sub>-0.33</sub></td></tr></tbody></table>

constrained from the analysis of Figure 10. For now, we only attempt to fit a purely Keplerian rotation profile, but given the large disk mass of Elias 2–27, fitting a self-gravitating rotation curve is warranted. This possibility will be further explored in a separate paper (Veronesi et al. 2021).

Finally, we use the constrained emission surface of each side of the disk, and a stellar mass value of 0.49M☉, to build a model of the expected mean velocity maps of each isotopolog. The models are compared to the observations through residual analysis. Our model velocity maps consider only the Keplerian motion of the upper layer of the emission surface. In the case of Elias 2–27, the disk is inclined such that emission from the lower layer appears in the southern part of the disk (Huang et al. 2018c), which may be cause for larger residuals across the southern border. Additionally, our constraints on the shape of the emission surface do not cover the whole radial extent of the emission and are extracted from the data retrieved at radial distances between ~40 and 300 au (~0′′35–3′′59) from the central star (see Figure 10). Therefore, we should also expect to have larger residuals in the inner and outer regions where the

emitting surface is not directly constrained by our method described above.

Figure 12 shows the ¹³CO and C¹⁸O velocity maps in the top and bottom rows, with observations, model, and residual velocity maps in the left, middle, and right columns. The integrated velocity maps are computed using the BETTERMOMENTS package (Teague & Foreman-Mackey 2019) to accurately constrain the line-of-sight velocity from the channel maps with 0.111 km s⁻¹ spectral resolution. Initial analysis of the observations allows us to identify marked perturbations throughout the disk, especially along the south and the west, where a distinct “distorted” pattern is observed in the outer disk. Along the major axis, the C¹⁸O data also displays “distorted” perturbations in a seemingly perpendicular form with respect to the azimuthal southern “distorted” pattern. In the model maps, we see that the west side of the disk is able to reproduce to some extent the “distorted” pattern along the south of the disk in both isotopologs, given the decrease in height in the outer disk. We do not see this pattern on the east side of the model map, as the retrieved emitting surface does not present large deviations from a linear cone-like model.

In the residuals, the radial extent that was used for determining the emission layer height profile is marked as two ellipses to define the inner and outer radial bounds (40–275 au). For larger radial distances than what was sampled, the residuals on the west side of the disk are much lower than the residuals on the east side, for both isotopologs. This means that the modeled west-side emission layer with a “dip” in the emitting surface height at larger distances from the star is necessary when extending the model to larger radii. This indicates that the east side of the disk likely also has a decrease in the emitting surface height at larger radial distances. Given our
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Top panels show the data tracing the velocity of the gas emission, as a function of radial distance to the star, from the C<sup>18</sup>O (left) and <sup>13</sup>CO (right) isotopologs. Blue points correspond to measurements coming from the west side of the disk, red points come from the east side, plotted curves correspond to the best-fit Keplerian rotation profile, and the shaded area corresponds to the stellar mass uncertainty as indicated by the 16th and 84th percentiles of the posteriors. The vertical dashed line indicates the location of the gap reported in the continuum, dotted-dashed line corresponds to the gap location in the C<sup>18</sup>O-integrated intensity map. The gray area indicates the width of the dust gap (obtained from Huang et al. 2018b), the orange area indicates the gas gap's location uncertainty. The bottom panels show the residuals of each isotopologue after subtracting the best-fit model to the data.</figcaption></figure>

limited range of sampled radial distances, we may not be sensitive to this “turning point” in the east. Within our sampled radii, marked by the ellipses, we observe negative residuals in the northeast quadrant of the <sup>13</sup>CO emission. This roughly coincides with the location of the most prominent cloud absorption (see Figure 8), so we associate these residuals with the absorption. The C<sup>18</sup>O residuals within the ellipses are much stronger and display a rough “X” shape across the center, with marked positive residuals close to the minor and major axes’ in the northwest and southeast quadrants, respectively. As was noted, the radial distances within the ellipses coincide with the radial extension of the dust spiral arms, while the “X” shape does not have a clear colocation with the spiral structure, the largest positive residuals coincide with the location where a spiral starts and the other ends. This “X”-shaped residual probably traces perturbations arising closer to the midplane, as it is not observed in <sup>13</sup>CO and the C<sup>18</sup>O emission traces a lower height layer.

4.4. Features in the Channel Maps of <sup>13</sup>CO and C<sup>18</sup>O

In the CO channel maps (see figures in Appendix A), we observe several perturbations, which do not follow the expected Keplerian velocity field, and we refer to them as “kinks.” In the following figures, we overlay several previously characterized features of Elias 2–27 to use as reference: the continuum spirals and their extension, the location of the dust continuum, and C<sup>18</sup>O gaps, together with the expected isovelocity curves for each channel. We note that the isovelocity curves seem perturbed as they are obtained from the model velocity emission shown in Figure 12, considering the constrained emission layer

of each isotopolog and disk side. We only show the isovelocity curves of the top layer of emission because the geometrical constraints have been derived for this layer and do not adequately trace the bottom layer, which we can visually identify. In the bottom layer, the emission seems to be coming from a layer at a farther distance from the midplane than what we trace with the geometry of the top layer. This difference should be studied in future work and could be caused by the layers tracing different sectors of the disk due to temperature effects (Pinte et al. 2018b), optical depth, or some other asymmetry in the vertical disk structure.

We observe two types of perturbations: inner kinks, at roughly the location of the spiral arms (but outside the dust gap at ~70 au), and outer kinks, beyond the extent of the continuum emission at ~250 au. Perturbations are strongly present in the central velocity channels of both CO tracers, shown in Figure 13. We observe the inner kink (marked with a yellow arrow; Figure 13) close to the spiral on the south side of the disk and along several channels; it appears strongest at channels +1.77 to +1.55 km s<sup>-1</sup>. This feature is colocated with the NW spiral. A large outer “C” shape (marked with a green arrow; Figure 13) can be seen beyond the gas gap in the south where the emission is brightest. This “C”-shape feature is strongest in <sup>13</sup>CO, along channels 1.55–2.33 km s<sup>-1</sup>, but in channels 1.66–1.88 km s<sup>-1</sup> of C<sup>18</sup>O, we can also recognize it in the southern part of the disk. We suggest that the “C” is not a deviation from Keplerian motion but rather is the projected emission from the upper and lower sides of the disk, connected by material bridging in the center between both sides.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Velocity maps, model, and residuals for $^{13}$CO (top row) and C$^{18}$O (bottom row). In each row, the first column shows the integrated emission velocity map (moment 1). The second column shows the velocity map model, computed using the constraints found for the emission surface and stellar mass. The third row shows the residuals calculated by subtracting the model map from the observations. The spiral arms show the best-fit parametric model from the 0.89 mm dust emission; inner and outer ellipses indicate radial limits for the data used to derive the emission surface geometry and stellar mass values.</figcaption></figure>

Besides these features located near the disk systemic velocity, we see more subtle deviations in high-velocity channels, for both east and west sides of the disk in the C$^{18}$O channel maps. To highlight these perturbations, we show the expected isovelocity curves for these velocity channels, along with the dust spirals in Figure 14. The top and bottom panels of Figure 14 show the west and east C$^{18}$O emission, imaged with a finer spectral resolution that is available only for the C$^{18}$O data. The deviations are most visible at 0.95–1.0 km s$^{-1}$ in the west, where the top layer of the disk emission does not precisely follow the isovelocity curve (blue line) and appears perturbed at the spiral arm location. On the east side, at 2.95–2.90 km s$^{-1}$, similar deviations are apparent in the top emission layer of the disk, roughly colocated with the SE spiral. These kinks are not clearly observed in the $^{13}$CO maps; however, it is expected that perturbations due to the spirals should be more apparent in C$^{18}$O than in $^{13}$CO, as the C$^{18}$O traces a layer closer to the midplane, where the spirals reside. The deviations are better discerned in the west channels, possibly due to the lack of cloud absorption at these velocities, but also because, if the kinks are caused by the spiral arms, the highest-contrast spiral is in the west side of the disk.

Recently, Pinte et al. (2020) reported the presence of a kink in the northern side of Elias 2–27, at the location of the dust gap, which was signaled as a possible indicator of a planetary companion. We do not recover this feature, possibly given our lower spatial resolution: the DSHARP data have 4–5 times better angular resolution than this work and is sensitive to spatial scales down to ~6 au in this system. This work analyzes data with higher spectral resolution (3–6 times better) and less affected by cloud absorption than previously published studies. The latter makes us sensitive to the perturbations reported in

this work; however, we are not able to detect small spatial scale perturbations due to our angular resolution.

## 5. Hydrodynamic Simulations of a Disk Undergoing GI

Elias 2–27 has been subject to different modeling approaches in order to explain the origin of the observed spiral substructure (Meru et al. 2017; Tomida et al. 2017; Bae & Zhu 2018b; Forgan et al. 2018b; Hall et al. 2018; Cadman et al. 2020), with most of the modeling efforts oriented toward a gravitationally unstable disk as the disk-to-star mass ratio ($q$) has been estimated to have values around 0.2–0.3 (Andrews et al. 2009; Isella et al. 2009; Ricci et al. 2010; Pérez et al. 2016) and GI's are expected when disk-to-star mass ratios are > 0.1 (Kratter & Lodato 2016). The number of spiral arms ($m$) excited in a GI scenario will depend inversely on the disk-to-star mass ratio ($m \sim M_{*}/M_{d}$). Given the $m=2$ spiral mode observed in Elias 2–27, previous simulations (Meru et al. 2017; Tomida et al. 2017; Forgan et al. 2018b; Hall et al. 2018) have aimed at producing the system using higher disk mass estimates ($M_{d}/M_{*} \sim 0.5$) than those derived from the observations ($M_{d}/M_{*} \sim 0.2–0.3$, Pérez et al. 2016). Recent work by Cadman et al., however, has shown that a disk-to-star mass ratio of 0.27 may also reproduce the observations. While this low disk-to-star mass ratio predicts a high amount of spirals, it has been shown that ALMA sampling can make a disk of $q=0.25$ appear as an $m=2$ system (Dipierro et al. 2014).

We performed a total of 10 three-dimensional, dusty, gaseous hydrodynamical simulations using the SPH code PHANTOM (Price et al. 2018). To accurately compare the simulations to our multiwavelength observations, we use the multigrain setup considering five different grain sizes, ranging from 1 µm to 1 cm in five logarithmically spaced size bins,
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Selected central channels of $^{13}$CO (top) and C$^{18}$O emission (bottom). The white continuous line shows the dust features: the inner gap at 69 au and the spirals as traced from the 0.89 continuum emission. The dotted white line traces the C$^{18}$O gas gap location at 241 au. Dashed lines show how the spirals traced in the dust would extend further outside of the continuum emission. The blue curve traces the expected isovelocity curve of each channel, following the constrained emission layer geometry of the top layer. The velocity of each channel map is indicated in the top-right corner of each panel; the beam size is in the bottom-left corner. Green arrows mark the outer perturbation; yellow arrows mark the inner perturbation.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Selected high-velocity channels of C$^{18}$O emission. White lines trace the spirals detected in the 0.89 mm continuum emission, blue lines indicate the isovelocity curves expected at each channel velocity, indicated in the top-right corner of each panel, following the constrained emission layer geometry of the top layer. Arrows indicate were deviations from expected isovelocity curves ("kinks") are observed.</figcaption></figure>

assuming a size distribution of $dn/da \propto a^{-3.5}$. Multiple grain sizes are necessary because, while the most efficient emission at wavelength $\lambda$ comes from dust grains of size $a \sim \lambda/2\pi$ (e.g., Draine 2006), there is an overall contribution from all grains. The dust is modeled self-consistently with the gas, using the multigrain “one-fluid” approach, where we limit dust flux using the Ballabio switch (Ballabio et al. 2018). Because the disk is massive and self-gravitating, the dust remains in the strongly coupled regime ($St \lesssim 1$) out to grain sizes of several centimeters, so we do not need to use the two-fluid approach. In this regime, the dust exerts a force back on the gas (back-reaction) that is

significant (Dipierro et al. 2018), so we include this effect on the gas.

In all 10 simulations, we used 1 million SPH particles and assumed a central stellar mass of $0.5M_{\odot}$, represented by a sink particle (Bate et al. 1995), with the accretion radius set to 1 au. We set the initial inner and outer disk radii to 5 au and 300 au, respectively. Different simulations vary in disk-to-star mass ratio and density profile. The total dust mass in the system is kept constant at 0.001 $M_{\odot}$, because this is observationally constrained (Pérez et al. 2016). As the total dust mass is a fixed value in our simulations to obtain different disk-to-star mass
ratios we vary the gas-to-dust ratio ($\epsilon$), which corresponds to $\epsilon = 100, 151$, and $252$ for $q = 0.2, 0.3$, and $0.5$, respectively. We do not sample lower gas-to-dust ratios because we do not expect to recover an $m = 2$ spiral arm morphology for $q < 0.2$. While sampling the emission with ALMA can make disks with high $m$ values appear as $m = 2$ systems, when going to very low $q$, this effect does not hold (Dipierro et al. 2014). The sound speed profile was set as $c_s \propto R^{-0.25}$, and we used two surface density profiles: either a simple power law,

$$ \Sigma(R) = \Sigma_0 \left(\frac{R}{R_0}\right)^{-p}, \qquad (4) $$

where $\Sigma_0$ is the surface density at the inner edge of the disk and $p$ either $1.3$ or $1.5$, or an exponentially tapered power law,

$$ \Sigma(R) = \Sigma_c \left(\frac{R}{R_0}\right)^{-p} \exp\left[-\left(\frac{R}{R_c}\right)^{2-p}\right], \qquad (5) $$

where $R_c$ is the characteristic radius of the profile, which we set to $R_c = 200$ au and $p$ either $0.7$ or $1.0$. In both surface density profiles, $R_0$ is the reference radius and is set to $R_0 = 10$ au.

We used a polytropic equation of state and assumed that the disk cooled through the $\beta$ cooling prescription (Gammie 2001), where the cooling timescale, $t_c$, is related to the dynamical timescale, such that $t_c = \beta t_{dyn}$. The dynamical timescale is the rotation period, $2\pi/\Omega$, and we set $\beta = 15$. Finally, each simulation is computed for 10 orbital periods at the outer radius (300 au), from which we receive outputs every 0.1 fraction of an orbital period. The details of the model parameters are shown in Table 3.

While SPH simulations portray the overall dynamic behavior of a system, in order to accurately compare the model to an observation, it is necessary to produce radiative transfer calculations of the SPH outputs, and then simulate mock observations using the same observing conditions (uv-coverage) as the actual observations. Radiative transfer is necessary because it accounts for the multiwavelength emission the grains will have given the stellar characteristics and grain distribution. Sampling the radiative transfer output with the same observing configuration is crucial, for the antenna distribution and observation time will result in a given angular resolution and will form an image according to the sampled uv-coverage.

## 5.1. Dust Simulations

The Monte Carlo radiative transfer MCFOST code (Pinte et al. 2006, 2009) was used to compute the disk thermal structure and synthetic continuum emission maps at each observed wavelength (0.89, 1.3, and 3.3 mm). We assumed $T_{gas} = T_{dust}$, and used $10^7$ photon packets to calculate $T_{dust}$. We set the parameters for the central star to match those of the Elias 2–27 system (Andrews et al. 2009; Pérez et al. 2016), with temperature $T = 3850$ K, $M = 0.5 M_\odot$, and $R_* = 2.3 R_\odot$.

To create the density structure as input into the MCFOST calculation, each SPH simulation underwent Voronoi tessellation such that each SPH particle corresponds to one MCFOST cell. We assumed the dust is a mixture of silicate and amorphous carbon (Draine & Lee 1984), and optical properties were calculated using Mie theory. The grain population consists of 100 logarithmic bins ranging in size from 0.03 µm to 1 mm. The dust density of a grain size $a_i$ was obtained by interpolating from the SPH dust sizes in each cell in the model. We assume that grains smaller than half the

<table><caption>**Table 3**
SPH Model Parameters</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>$M_*(M_\odot)$</td><td>0.5</td></tr><tr><td>$R_{in}$ (au)</td><td>5</td></tr><tr><td>$R_{out}$ (au)</td><td>300</td></tr><tr><td>$M_{dust}$ ($M_\odot$)</td><td>0.001</td></tr><tr><td>Gas-to-dust mass ratio</td><td>100, 151, 252</td></tr><tr><td>Min. grain size (cm)</td><td>10<sup>-4</sup></td></tr><tr><td>Max. grain size (cm)</td><td>1</td></tr></tbody></table>

smallest SPH grain size (0.5 µm) are perfectly coupled to the gas distribution. We normalized the dust size distribution by integrating over all grain sizes, where a power-law relation between grain size $a$ and number density of dust grains $n(a)$ was assumed such that $dn(a) \propto a^{-3.5}da$.

The radiative transfer emission map of each wavelength was sampled with the same uv-coverage as the observations at each corresponding wavelength, using GALARIO (Tazzari et al. 2018) to create mock ALMA visibilities that were afterward processed using the same deconvolution procedures as with the observations (described in Section 2) to obtain the final mock ALMA images.

For each simulation image, we subtract its azimuthally averaged radial profile of emission, following the same procedure described to trace the spiral morphology on the multiwavelength observations (Section 3.1). We find that most of our models are able to accurately reproduce the $m = 2$ large-scale spiral morphology. This is expected in the ALMA images, even at lower disk-to-star mass ratios (where we expect a larger number of spiral arms), as was shown by Dipierro et al. (2014). To select the simulation that best resembles our observations, we measure the spiral’s pitch angle and compare them to the observational values. The pitch angles measured for each spiral and model setup are shown in Table 4 in Appendix B. The model pitch angle ($\phi_{model}$) and the observational pitch angle ($\phi_{obs}$) are compared considering the difference between the values of each spiral and wavelength and weighing by the error of the observational pitch angle ($\sigma_{obs}$), following a likelihood parameter determined by

$$ \sum_{\lambda} \sqrt{\left[ \left( \frac{\phi_{\text{model}}^{\text{NW}} - \phi_{\text{obs}}^{\text{NW}}}{\sigma_{\text{obsNW}}} \right)^2 + \left( \frac{\phi_{\text{model}}^{\text{SE}} - \phi_{\text{obs}}^{\text{SE}}}{\sigma_{\text{obsSE}}} \right)^2 \right]}_{\lambda}. \qquad (6) $$

The simulation that best reproduces the observations follows a dust density profile of an exponentially tapered power law with index 1.0 and a disk-to-star mass ratio of 0.3. These parameters are close to the previously published observational constraints for this disk (Pérez et al. 2016) and similar to the result derived in Cadman et al. The simulation is shown in Figure 15. Besides reproducing the pitch angle values, the radial extension of the spirals and their overall morphology in the simulation are similar to the observations. We note that the comparison to each simulation set was made for a specific time step within all simulation outputs, and the selection was based on the output that showed a clear two-spiral-arm feature after at least six outer orbits. For the case of the best-fit simulation, the selected time step was after 6.4 outer orbits (at 300 au from the star). It is important to state that, as expected, the spirals arising from GI are constantly excited and de-excited throughout the
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Panels from left to right correspond to data from the 0.87, 1.3, and 3 mm simulated observations for an exponentially tapered dust density profile with index 1.0 and disk-to-star mass ratio *q* = 0.3. Top: images of simulated emission with the azimuthally averaged intensity profile subtracted, blue and red dots trace the maxima location along the spirals. Bottom: blue and red dots correspond to the deprojected radial location of the traced spirals in the simulated emission. Colored solid lines show the constant pitch angle logarithmic spiral fit, dashed colored lines extend the fit to lower radii. Black points are the deprojected radial location of the spirals from the observations (Section 3) and their astrometric error. The pitch angle likelihood parameter is indicated in the bottom-right panel.</figcaption></figure>

time-lapse of our simulation. This means that for the same set of parameters, there may be several time steps that accurately reproduce the morphology but others where no spirals are seen.

Even though there is a simulation that reproduces the observations better than the rest, the pitch angles of the different SPH simulations are in most cases comparable to those of the observations (~12°.9 and ~13°.2, for the NW and SE spirals, respectively, with small variations between wavelengths). This shows that it is possible to reproduce the grand-design spirals even at lower disk-to-star mass ratios than previously tested in this system (Meru et al. 2017; Hall et al. 2018), with stellar mass, disk dust mass, and density profile values comparable to the observational constraints. The likelihood value of the rejected models is shown in Table 4 in Appendix B. Additionally, we see that not all GI spirals, recovered from the models and measured with our method (Section 3.1), are perfectly symmetric (same pitch angle), which is a property that has been predicted in other works for GI-excited spiral arms (Forgan et al. 2018b). This is specifically observed in the models with *q* = 0.2.

We compare our best-fit model with the high-angular-resolution DSHARP (Andrews et al. 2018) data, at 1.3 mm. The comparison of the subtracted, residual images is shown in Figure 16. From the visual comparison, we clearly see that the internal structure of the spirals is different. Even at high angular resolution, the observations show thicker and wider, continuous spirals. On the other hand, the simulation shows thinner and discontinuous spirals, which are made from a superposition of filaments. At larger radii, we note that the spirals of the observation get wider, while in the simulation they remain thin. The causes for these differences are discussed in Section 6.3.

We measure the pitch angle value, sampling the same angular extent as in the work by Huang et al. (2018c). We retrieve a value of 15°.56 ± 0°.06 in the NW spiral and 12°.76 ± 0°.06 in the SE spiral. These values differ from the constraint from the lower angular resolution simulations, showing that the beam smearing does impact the pitch angle measurements, as proposed in Section 3.1. The pitch angle values for the DSHARP observations are 15°.7 ± 0°.2 in the NW spiral and 16°.4 ± 0°.2 in the SE spiral. The main difference is from the SE spiral, possibly related to the effects of the morphology differences and the thickness of the spiral from the observation.

## 5.2. Gas Simulations

From the simulation that best reproduces the pitch angle of the spiral arms (*q* = 0.3, density profile with a tapered power law of index 1.0), we compute the simulated channel maps. As with the dust simulations, we use MCFOST (Pinte et al. 2006, 2009), with the same parameters as before, to compute the disk thermal structure and synthetic $^{13}$CO $J=3-2$ and $C^{18}$O $J=3-2$ line maps. The molecule abundances, relative to local H$_2$ are set to $7 \times 10^{-7}$ for $^{13}$CO (as done in Hall et al. 2020) and $2 \times 10^{-7}$ for $C^{18}$O (following the estimate from Frerking et al. 1982). The spectral resolution for the kinematic simulations is set to 0.111 km s$^{-1}$ to match the observations, and we additionally compute $C^{18}$O channel maps with 0.05 km s$^{-1}$ resolution to compare with the finer spectral resolution data.

The results for two representative channel maps are shown in Figure 17. While in the radiative output we observe the characteristic “GI-wiggle” shown by Hall et al. (2020) at the spiral
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Residual emission after subtracting the azimuthally averaged intensity profile. Left: DSHARP (Andrews et al. 2018) high-angular-resolution emission at 1.3 mm. Right: simulated emission of a GI disk with an exponentially tapered dust density profile of index 1.0 and disk-to-star mass ratio *q* = 0.3. Both data sets have the same angular resolution, indicated by the beam in the bottom left of each panel.</figcaption></figure>

arm's location, when sampling the data with the observation's uv-coverage, the “GI-wiggle” is not visible in any isotopolog. In the work by Hall et al. (2020), the “GI-wiggle” remains visible even after convolving with a Gaussian beam. Compared to Hall et al. (2020), our simulated ALMA images have a ~4 times lower spectral resolution (0.111 km s⁻¹ compared with 0.03 km s⁻¹) and ~3 times lower angular resolution (0″/3 compared to 0″/1). While the higher-spectral-resolution C¹⁸O observations have a comparable channel width (0.05 km s⁻¹) with the analysis of Hall et al. (2020), the angular resolution smears the “GI-wiggle” features (see Figure 18). We note that no north/south brightness asymmetry is present in the simulated channel maps. Additional channel maps computed with a 90° inclination do not show any significant vertical difference between the east/west sides.

# 6. Discussion

## 6.1. Spiral Structure and Multiwavelength Dust Continuum Emission

The morphology of the dust spiral structure can be a key indicator of the origin of the spirals. We measure symmetric spirals, present in all three wavelengths, with similar radial extension and pitch angles. Symmetric spirals with constant pitch angles are predicted in the case of GI (Forgan et al. 2018b), rather than for companion perturbations where asymmetric, variable-pitch-angle spirals are expected (Bae & Zhu 2018b). Additionally, we measure similar contrasts for both spirals (~30% difference in contrast between the NW and SE spirals), which is also consistent with GI predictions, as companion-induced spirals are expected to show a clear primary spiral (Bae & Zhu 2018b). The latter was already shown for the emission at 1.3 mm (Pérez et al. 2016; Huang et al. 2018c). In this study, we extend the finding to 0.89 and 3.3 mm. Additionally, the spectral index map shows spiral morphology with slightly lower alpha values along the spirals.

This coincides with the prediction for dust trapping, expected for GI (Dipierro et al. 2015).

We measure the optical depth profiles, which appear to show similar values in all three wavelengths. However, the temperature profile used for deriving the optical depth (computed using the flaring and stellar luminosity values from Andrews et al. 2018) results in a midplane temperature ~2 times higher than the 0.89 mm brightness temperature at all radii. Most likely, the disk is much colder and optically thick than what we derive, and our optical depths must then be taken as lower limits.

Several works have shown that when scattering from dust grains is considered, optically thick disks can display lower intensities and be categorized as more optically thin disks (Liu 2019; Zhu et al. 2019; Sierra & Lizano 2020). When scattering is not present and the emission is optically thin, then $\alpha$ cannot reach values below 2.0. When dust scattering is included, a region of high optical depth can have $\alpha$ < 2.5 (and even attain a spectral index below 2.0) if the albedo decreases with wavelength (Zhu et al. 2019; Sierra & Lizano 2020), something observed in the innermost regions of TW Hya (Tsukagoshi et al. 2016; Huang et al. 2018a). Furthermore, while $\alpha$ also depends on other dust properties, such as the grain size distribution (e.g., Testi et al. 2014), values of $\alpha$ < 3 are not expected when the optical depth is low, and even in presence of centimeter-sized grains, $\alpha$ should not attain values below ~2.5 (Zhu et al. 2019).

From our spectral index profiles $\alpha$ < 2.0 in the inner ~40 au and $\alpha$ < 2.5 inside ~200 au, outside this region $\alpha$ grows, reaching a maximum value of 3.0 in the outer disk. This indicates that the outer disk is probably optically thin with grains of 0.1–10 cm, favoring a dust distribution $n(s) \propto s^{-3.5}$ (see Figure 9 in Zhu et al. 2019). However, in the region between ~70–200 au the spectral index increases slowly between ~2.2–2.5. This scenario favors a dust distribution $n(s) \propto s^{-2.5}$ and can be explained by either optically thin
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Individual emission of channels at velocities +1.73 km s⁻¹ (left) and +1.62 km s⁻¹ (right) for simulated ¹³CO (top two rows) and C¹⁸O (bottom two rows) emission. From top to bottom, the first and third rows correspond to the MCFOST output emission. The second and fourth rows show the emission after applying uv-coverage as in observations and processing with CASA. White lines trace the spirals from the best simulation ($q = 0.3$, exponentially tapered dust density profile index 1.0). The beam for the simulated ALMA images is shown in the bottom left of each corresponding panel.</figcaption></figure>

emission and 0.1–10 cm or optically thick emission with maximum grain sizes ~0.1 cm. For the inner disk (~70 au), the spectral index reaches values below 2.0. Most likely, the emission is optically thick and dust scattering is at work, even when the maximum optical depth we constrain (under standard assumptions) is $τ \sim 0.5$ at all wavelengths.

For the DSHARP sample, it has been shown that the optical depths of $τ \sim 0.6$ at 1.3 mm in bright rings can be obtained from optically thick regions with a scattering albedo of $ω_ν \sim 0.89$ (Zhu et al. 2019). For Elias 2–27, we measure an average $τ \sim 0.45$ at 1.3 mm inside 70 au, which could be obtained with a scattering albedo of $ω_ν \sim 0.93$ (see Equations (14) and (15) from Zhu et al. 2019). This albedo value is sufficient to mask optically thick 1.3 mm emission ($τ_{real} \sim 1–5$)

in the inner disk, which we would be inferring as optically thin 1.3 mm emission ($τ_{obs} \sim 0.45$) in the standard assumption of only absorption opacity. As was previously discussed, our optical depth values could be underestimated due to the difference between the measured brightness temperatures and model midplane temperature. If the midplane temperature was a factor of ~1.5 lower, the optical depth in the inner 70 au of the disk, at 1.3 mm would be $τ \sim 0.99$, which could be obtained with a scattering albedo $ω_ν \sim 0.74$. Similar analysis can be done with the other wavelengths from the measured $τ$ values of the inner disk. This results in a scattering albedo of 0.92 and 0.95 for 0.89 mm and 3.3 mm, respectively.

Together with the effect of scattering at long wavelengths, the observed low spectral index values can also occur due to the disk’s temperature. For cold (<30 K) systems, such as Elias 2–27, spectral indices below 2 are expected due to the displacement of the peak blackbody radiation to the sub-millimeter range (Sierra & Lizano 2020). Additionally, if the emission is optically thick and there are fluctuations in the vertical temperature structure, this could also result in $α < 2.0$ (see Figure 5 in Sierra & Lizano 2020).

The underestimation of the optical depth in Elias 2–27 when not including the effect of scattering will impact its solid mass estimate. This effect is larger for inclined disks and when the emission area is compact (Zhu et al. 2019), which is the case of our source in the inner regions. Considering the large disk extent (up to ~250 au), most of the dust mass resides in the optically thin outer disk, thus, the disk mass could be underestimated by up to a factor of ~2 (Zhu et al. 2019). The latter implicates that the previously constrained disk-to-star mass ratio of 0.1–0.3 (Andrews et al. 2009; Ricci et al. 2010, using standard assumptions such as a gas-to-dust ratio of 100) is a lower bound. If the disk mass was higher by up to a factor of 2, the resulting disk-to-star mass ratio would make GIs a likely cause for the spiral structure.

## 6.2. Asymmetries and Perturbations in the Gas

Our analysis of new ¹³CO and C¹⁸O (see Figures 19 and 20) allows us to constrain a highly perturbed morphology for the emitting gas layer in Elias 2–27 is a new characteristic for this system and offers new insight into the ongoing dynamic processes. The asymmetric structure of the ¹³CO- and C¹⁸O-emitting layer, as well as the dust spiral arms present in Elias 2–27, could be in principle caused by fly-by interaction (e.g., Cuello et al. 2019) or an external companion. But if this were the case, we also expect a strong kinematical perturbation in the integrated emission and velocity maps of Elias 2–27, such as those reported by Kurtovic et al. (2018) for disks with known external companions. Furthermore, observations in near-infrared of Elias 2–27 have not found any companion (Cieza et al. 2009; Zurlo et al. 2020).

With no outer perturbation, the emission layer height of the gas should follow hydrostatic equilibrium, and its structure depends on the gas temperature (Armitage 2015). In this case, we expect a cone-like or flared emission layer, with height increasing at a larger radial distance (Rosenfeld et al. 2013). We will discuss two possible origins for the observed asymmetries in the gas: ongoing infall of material from its surrounding cloud/envelope, or a warped inner disk causing azimuthal temperature variations.

Three-dimensional simulations of circumstellar disks with ongoing accretion show that the vertical structure of the disk
separated by the semiminor axis, we have discussed that this is
probably by chance and that we cannot determine an exact
symmetry angle. The method used to derive the emitting
surface assumes symmetry with respect to the semimajor axis
and this biases our results. We also show that, when tracing the
deprojected border of emission (see Figures 9 and 23), there are
two local maxima and minima radial extensions, suggesting
that there is not in fact a single symmetry axis (we would have
expected only one minimum and maximum extent, roughly
symmetrically opposed for the case of a warp; Juhász &
Facchini 2017). More complex processes, or a combination of
effects, are occurring.

Finally, if a warp is responsible for the kinematic effects,
there is the question of its origin. On one hand, warps are
thought to arise from close-in binary interactions or inclined
planetary orbits (Nealon et al. 2018; Aly & Lodato 2020).
Warps in very young systems, possibly produced by the infall
of material, have been predicted (Bate et al. 2010) and also
reported (Sakai et al. 2019). In the case of Elias 2–27, we have
discussed that infall may be occurring, given the detection of
large-scale emission at velocities close to the disk velocities. If
the disk is indeed warped, this effect could be caused either by
a planet in an inclined orbit or infall, both options require
further investigation.

Aided by the isovelocity curves computed from the
constrained emission surface of the disk, we detect multiple
deviations from Keplerian motion in the velocity cube of ¹³CO
and C¹⁸O, colocated with the spiral structure. The isovelocity
curves themselves show a perturbed nature, given the complex
emission surface. The colocation of the perturbations with the
spiral structures and the strong deviations could indicate a
connection between the spirals and the emitting surface
morphology. Previously reported deviations for planetary
companions have been around 15% (Pinte et al. 2018a) with
respect to the channel velocity, while the perturbations we
observe are much larger, reaching up to 80% of the channel
velocity for some of the features in the central channel maps.
Such large perturbations increase the likelihood that they relate
to the spirals, rather than to a companion. Large planetary
deviations require a high perturber mass, and in that case, we
would expect to see its effect on the dust emission.
Furthermore, we expect deviations caused by a planet to be
spatially localized (Pinte et al. 2018a, 2019, 2020) and our
observations show deviations on both the east and west sides of
the disk, present along several channels. The detected
perturbations agree with the predictions of Hall et al. (2020);
they are colocated with the spirals and the morphology of the
kink is similar to what was predicted in their work. If indeed
the disk is warped or suffering considerable infall of material,
as previously discussed, the observed “kinks” could be the
combination of both kinematical deviations induced by a warp
and perturbations due to the spirals.

In this study, we attempted a purely Keplerian fit to the
expected super-Keplerian velocity profile when GI is the
governing process (Bertin & Lodato 1999; Lodato 2007). The
Keplerian rotation curve is able to fit the observed velocity
profile, but we do note that the east channels, especially from
the C¹⁸O emission, show super-Keplerian velocities. Further
analysis is needed to check if a self-gravitating rotation curve
may be a better description of this data. This is currently being
studied and will be presented in a future publication (Veronesi
et al. 2021).

Finally, regarding the gap in the C¹⁸O gas emission, we do
not see evidence of it being produced by a physical perturber.
The gap does not appear to be colocated with any perturbation
in the channel maps (see Figure 13), which we would expect if
the gap origin was planetary (e.g., Pinte et al. 2018a; Teague
et al. 2018). It may be due to chemical processes of the gas or
optical depth effects (see discussion in Guzmán et al. 2018).

6.3. *Comparison with SPH Simulations*

We find that GIs can accurately reproduce the spiral
morphology at multiple wavelengths, with parameters close
to the observational constraints, as shown in Figure 15.
However, we see considerable morphological differences in the
comparison to high-angular-resolution data. The width and
morphology of spiral arms in GI environments can be regulated
by varying the cooling parameters (see figures from Rice et al.
2003; Boss 2017). Future SPH simulations should attempt to
sample this parameter space to further understand the cooling
prescription of the disk. We note that Cadman et al. also
analyzes the spiral structure of Elias 2–27 at DSHARP angular
resolution, and they obtain thicker spirals in the subtracted
image. These spirals are constructed with a semianalytical
model, considering grain growth, and while they are thicker
than what we recover, they are still not able to reproduce the
precise morphology of the observations, as they are too smooth
(see subtracted images in Figure 16 from Cadman et al. 2020).

GI does not explain the dust gap, if the dust gap was carved
by a planet, estimates indicate at best a mass of 0.1M<sub>j</sub> (Zhang
et al. 2018), which is much smaller than planets expected to be
formed in a GI environment. Though not shown in this work,
we produced additional SPH simulations of an unstable disk
with different planetary-mass companions. In all cases, after a
few orbits, the planet migrated onto the star. This is expected,
according to the predictions of Type I migration for planets
orbiting several astronomical units from the star (Baruteau et al.
2014). Not allowing migration could allow a gap to form (see
Meru et al. 2017); however, it would not be a realistic scenario.

In our spectral line observations, we recover perturbations
that are colocated with the spirals and span several channels,
sharing some similarities to those predicted by GI disks (Hall
et al. 2020). Our gas simulations from the best-fit GI
parameters do not show the perturbation features from Hall
et al. (2020) when sampled and imaged with the uv-coverage of
the observations. This implies that the inner perturbation
observed in the ¹³CO and C¹⁸O channel maps is indeed very
large, as it appears at our lower angular resolution. The
perturbation is seen across several channels, which is in
agreement with the predicted perturbations of a GI disk (Hall
et al. 2020). The fact that we do not see the perturbation in the
channel maps from the simulated emission could be due to the
decisions regarding the chosen time stamp for the simulation
and the position angle. As we select a specific time frame, it is
possible that at other evolutionary stages, the strength of the
gas perturbations would have been larger. Also, it has been
shown by Hall et al. (2020) that the strength of the “GI-wiggle”
will vary depending on the position angle of the emission.
While the inclination of the simulation matches the value of the
observations (56°2), the position angle is set before inclining
the disk. When producing the ALMA mock image, the position
angle of the emission is determined by the value that allows
a good visual comparison of the final, inclined, simulated
emission image to the observations. Due to the computational
expense of sampling various position angles, we adopted the observational value (118°.8) for the position angle of the simulations, varying only by 90° or 180° in some cases. These shifts were decided using a visual criteria. We therefore note that there could be a range of position angle values that allow a good comparison with the observations, while also producing stronger kinematic perturbations, this is not tested in this work.

## 6.4. Spiral Structure Origin

The hypothesis for the origin of the spiral arms observed in the dust emission of Elias 2–27 are either GIs or perturbation by a companion. From the observational constraints presented in this work, there are some key features that may help us define which scenario fits best the Elias 2–27 disk. To begin, however, we must state that neither option can accurately predict, on its own, both the spirals and the dust gap. In the case of a disk undergoing GI, gaps are not expected features even if a planet was formed, given the fast migration even massive planets will have under these conditions (Baruteau et al. 2011). In order to reproduce the spiral arms, a companion would have to be located beyond the spiral extent (Meru et al. 2017). The possibility of the spirals being formed by a companion seems unlikely, as has been discussed in previous studies, given the contrast, symmetry, and extent of the spirals (Bae & Zhu 2018b; Forgan et al. 2018b). Furthermore, the possibility of an external perturber, such as a stellar companion or a fly-by causing the spiral structure, is also unlikely, due to the lack of a clear kinematical signature in the data and the nondetection of any object nearby (Ratzka et al. 2005; Cieza et al. 2009; Launhardt et al. 2020). Launhardt et al. (2020) conducted a NACO/VLT survey in search of planetary companions around 200 stars, including Elias 2–27. For this system, they reach a 50% probability of detecting a $2M_J$ companion outside 100 au or a ~10$M_J$ companion at 40–50 au (Launhardt et al. 2020, R. Launhardt 2021, private communication). We note that the detection of an external perturber is made more challenging by the extinction that affects this region.

Elias 2–27 has been shown to have a large disk-to-star mass ratio (Andrews et al. 2009; Ricci et al. 2010), in this work, we additionally discuss the possibility of the disk mass being up to a factor of ~2 higher if scattering is a relevant process (Zhu et al. 2019). Even if scattering was not relevant, we show that the disk is probably more optically thick than the reported values, which also points toward a larger disk mass. The values of disk-to-star mass ratio are sufficient to excite GIs, and we can accurately reproduce the spiral morphology using SPH models of a self-gravitating disk in medium resolution at least. Additionally, we detect tentative dust-trapping signatures in the continuum observations; the contrast curves show variations with increasing wavelength and lower spectral index values along the spirals. We also measure strong kinematic perturbation colocated with the spirals over multiple channels. The high disk mass, together with the strong deviations from Keplerian motion, consistent with the kinematical prediction for a GI disk (Hall et al. 2020), leads us to signal the origin of the spirals to be GIs, rather than a companion.

Infall of material would explain the high disk-to-star mass ratio of the system and the excitation of spiral structures due to GI (Hennebelle et al. 2017; Tsukamoto et al. 2017). While infall mechanisms are expected to be present in younger Class 0/I systems than Elias 2–27, which is classified as a Class II

disk (Andrews et al. 2009), it has been proposed that Elias 2–27 could be an extremely young Class II object to explain the spiral structure (Tomida et al. 2017). Though GI can explain the spiral structure traced in the dust, there is also a clear dust gap (Huang et al. 2018b), which emphasize is not a feature predicted or explained by GI. On the other hand, infall may explain on its own the perturbed morphology of the gas emission layer, but there is also a marked brightness asymmetry, which could be related to the presence of a warp in the disk (Nixon & Pringle 2010). If the disk is warped, due to infalling material breaking the disk, this could possibly explain the dust gap observed at high angular resolution by Huang et al. (2018b), if the separation between inner and outer misaligned disks were located at ~70 au from the central star. Certainly, further observations on the source are required to sample shorter baselines and adequately study the dynamics of the large-scale emission and search for infall signatures, and also observations at a higher spatial resolution, to analyze the origin of the dust gap and better constrain the kinematic perturbations in the system.

## 7. Summary

We have presented and analyzed multiwavelength dust continuum observations of the protoplanetary disk around Elias 2–27. We also studied the gas emission for $^{13}$CO and C$^{18}$O $J=3-2$. This provides new observational constraints on this source, which allowed us to study the origin of the prominent spiral structure. Our findings are as follows:

1. The spiral substructure is present in dust observations at multiple wavelengths, from 0.89 to 3.3 mm, and shows a higher contrast at longer wavelengths. These signs are possible indicators of grain growth and dust trapping at the spiral arm location.

2. From the spectral index analysis we trace a spiral morphology with lower spectral index values along the spiral location. This is expected for dust trapping, which is a key signature of gravitational instabilities not observed before in other systems with spiral morphology. The spectral index values also indicate the presence of large grains in the outer disk. Inwards of ~70 au the spectral index drops to values even lower than 2. This can be explained if this region has high optical depth and high albedo in the presence of dust scattering. We discuss that our optical depth estimates are lower limits, given the low brightness temperature measured in the observations. The presence of scattering and higher optical depths implies that the solids mass estimated under standard assumptions is likely a lower limit. The mass of the disk in Elias 2-27 is presumably higher than previously estimated.

3. We compute SPH simulations of a gravitationally unstable disk with parameters as those in Elias 2–27, and are able to replicate the spiral morphology at the three different wavelengths we study at ~0.′′2 resolution. Discrepancies at high angular resolution could be due to the cooling prescription used.

4. Observations show that the gas emission is not azimuthally symmetric in the vertical direction, i.e., the disk has a larger emission layer height in the west than in the east at most radial distances. Additionally, at larger radial distances, the kinematic data indicate that the emission layer height decreases. This is the first time we
observe an azimuthal emission layer height difference in a protoplanetary disk, and it does not appear in predictions for either a GI disk or the presence of a planetary companion.

5. Tracing the different heights of the ¹³CO and C¹⁸O emission layers, we show that ¹³CO comes from a higher layer than C¹⁸O, with velocities consistent with Keplerian rotation. The stellar mass we constrain (~0.46–0.5 M☉), is in agreement with the literature value (0.49M☉). Gas emission depletion (a gap) is observed in the distribution of C¹⁸O at a radius of ~240 au. This gap does not appear to be colocated with the main perturbations we recognize in the channel maps.

6. We see “kinks” or perturbations in the channel maps of both CO tracers that appear colocated with the spiral features. These kinks are stronger and present across a wide velocity range, making it unlikely that they have a planetary origin. The characteristics of these perturbations are similar to what has been predicted in a GI disk by Hall et al. (2020).

7. Based on observations that show large-scale emission surrounding and connecting to the disk, we propose the infall of material from the surrounding cloud is responsible for exciting GI in the disk, which in turn causes the dust spiral arm features. Infall of material can also explain the perturbed emission layer constrained from the gas tracers. Additionally, if infall warped the disk, this could explain the brightness asymmetry in the channel and integrated emission maps. Depending on the warp location, it could also explain the dust gap observed at higher angular resolution. Further observations are necessary to effectively detect the presence of a warp and confirm ongoing infall.

This paper makes use of the following ALMA data: #2013.1.00498.S, #2016.1.00606.S, and #2017.1.00069.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. L.M.P. acknowledges support from ANID project Basal AFB-170002 and from ANID FONDECYT Iniciación project #11181068. M.B. acknowledges funding from ANR of France under contract number ANR-16-CE31-0013 (Planet Forming Disks). C.H. was a Winton Fellow, and this research was supported by Winton Philanthropies/The David and Claudia Harding Foundation. A.S. acknowledges support from ANID/CONICYT Programa de Astronomía Fondo ALMA-CONICYT2018 31180052. J.M.C. acknowledges support from the National Aeronautics and Space Administration under grant No. 15XRP15_20140 issued through the Exoplanets Research Program. S.M.A. acknowledges funding support from the

National Aeronautics and Space Administration under grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. J.B. acknowledges support by NASA through the NASA Hubble Fellowship grant \#HST-HF2-51427.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. Th.H. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. L.L. acknowledges the financial support of DGAPA, UNAM (project IN112820), and CONACYT, México. M.T. has been supported by the UK Science and Technology research Council (STFC) via the consolidated grant ST/S000623/1. L.T. acknowledges support from the Italian Ministero dell'Istruzione, Università e Ricerca through the grant Progetti Premiali 2012—iALMA (CUP C52I13000140001), by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Ref No. FOR 2634/1 TE 1024/1-1, and the DFG cluster of excellence Origins (www.origins-cluster.de). 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 (DUSTBUS-TERS) and from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02). This research used the ALICE2 High Performance Computing Facility at the University of Leicester. This research also used the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. This work was partially supported by the University of Georgia Office of Research and the Department of Physics and Astronomy.

*Software:* CASA (McMullin et al. 2007), eddy (Teague 2019), bettermoments (Teague & Foreman-Mackey 2019), PHANTOM (Price et al. 2018), mcfost (Pinte et al. 2006, 2009), galario (Tazzari et al. 2018), frankenstein (Jennings et al. 2020), Astropy (Astropy Collaboration et al. 2013, 2018), Matplotlib (Hunter 2007), emcee (Foreman-Mackey et al. 2013).

## Appendix A
### Additional Gas Analysis Figures

The channel map emission used for constraining the vertical disk structure (Figures 19 and 20) and the complete FOV where large-scale emission is apparent (Figures 21 and 22) are presented. Figure 23 shows the radial extent of the integrated emission for both tracers.
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Individual emission of each channel for $^{13}$CO emission; corresponding velocities are written in the top right of each panel.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Individual emission of each channel for C$^{18}$O emission; corresponding velocities are written in the top right of each panel.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Individual emission of each channel for $^{13}$CO emission; corresponding velocities are written in the top right of each panel. Channels correspond to velocities where large-scale emission is observed. Imaging enhances large-scale structure; parametric model of the spiral arms is plotted for reference.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 22.** Individual emission of each channel for C$^{18}$O emission; corresponding velocities are written in the top right of each panel. Channels correspond to velocities where large-scale emission is observed. Imaging enhances large-scale structure; parametric model of the spiral arms is plotted for reference.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 23.** Radial extension of the integrated emission map for $^{13}$CO (in blue) and C$^{18}$O (in red), obtained for each azimuthal angle as the radius that encloses 98% of the azimuthally integrated emission.</figcaption></figure>

## Appendix B
### Additional SPH Model Details

Table 4 shows the results for the complete set of simulations, varying disk-to-star mass ratio and density distribution. The retrieved values for the pitch angle of each spiral at each observed wavelength are presented. Additionally, the likelihood parameter is shown in the last column of the table.

<table><caption>**Table 4**
Pitch Angle Values for Simulations</caption><thead><tr><th>Simulation</th><th>NW<sub>0.89 mm</sub></th><th>SE<sub>0.89 mm</sub></th><th>NW<sub>1.3 mm</sub></th><th>SE<sub>1.3 mm</sub></th><th>NW<sub>3.3 mm</sub></th><th>SE<sub>3.3 mm</sub></th><th>Likelihood φ</th></tr></thead><tbody><tr><td>Tapered: 0.7</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>q = 0.2</td><td>13°46 ± 0°08</td><td>14°96 ± 0°08</td><td>10°86 ± 0°06</td><td>15°27 ± 0°08</td><td>9°85 ± 0°07</td><td>13°7 ± 0°1</td><td>29.31</td></tr><tr><td>q = 0.3</td><td>12°52 ± 0°06</td><td>11°54 ± 0°06</td><td>12°24 ± 0°06</td><td>10°82 ± 0°05</td><td>12°09 ± 0°06</td><td>10°6 ± 0°06</td><td>24.04</td></tr><tr><td>q = 0.5</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Tapered: 1.0</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>q = 0.2</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>q = 0.3</td><td>12°75 ± 0°06</td><td>11°42 ± 0°05</td><td>12°34 ± 0°06</td><td>11°33 ± 0°05</td><td>13°08 ± 0°06</td><td>11°19 ± 0°05</td><td>21.04</td></tr><tr><td>q = 0.5</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Power law: 1.3</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>q = 0.2</td><td>13°9 ± 0°08</td><td>11°16 ± 0°11</td><td>14°75 ± 0°07</td><td>10°40 ± 0°08</td><td>14°82 ± 0°07</td><td>10°48 ± 0°11</td><td>36.19</td></tr><tr><td>q = 0.3</td><td>14°88 ± 0°07</td><td>15°66 ± 0°07</td><td>15°53 ± 0°06</td><td>15°36 ± 0°07</td><td>15°34 ± 0°07</td><td>15°92 ± 0°06</td><td>45.32</td></tr><tr><td>q = 0.5</td><td>14°86 ± 0°07</td><td>18°08 ± 0°07</td><td>14°47 ± 0°06</td><td>17°74 ± 0°08</td><td>14°2 ± 0°06</td><td>18°82 ± 0°07</td><td>53.08</td></tr><tr><td>Power law: 1.5</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>q = 0.2</td><td>10°94 ± 0°09</td><td>7°34 ± 0°06</td><td>12°38 ± 0°09</td><td>5°8 ± 0°07</td><td>11°49 ± 0°08</td><td>8°24 ± 0°07</td><td>66.08</td></tr><tr><td>q = 0.3</td><td>12°03 ± 0°06</td><td>11°95 ± 0°07</td><td>11°55 ± 0°06</td><td>10°4 ± 0°06</td><td>11°94 ± 0°06</td><td>11°05 ± 0°07</td><td>27.36</td></tr><tr><td>q = 0.5</td><td>11°04 ± 0°06</td><td>10°76 ± 0°05</td><td>10°62 ± 0°06</td><td>10°74 ± 0°05</td><td>10°42 ± 0°05</td><td>10°5 ± 0°05</td><td>46.58</td></tr></tbody></table>
ORCID iDs

T. Paneque-Carreño https://orcid.org/0000-0002-4044-8016

L. M. Pérez https://orcid.org/0000-0002-1199-9564

M. Benisty https://orcid.org/0000-0002-7695-7605

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

A. Sierra https://orcid.org/0000-0002-5991-8073

J. M. Carpenter https://orcid.org/0000-0003-2251-0602

S. M. Andrews https://orcid.org/0000-0003-2253-2270

Jaehan Bae https://orcid.org/0000-0001-7258-770X

Th. Henning https://orcid.org/0000-0002-1493-300X

W. Kwon https://orcid.org/0000-0003-4022-4132

L. Loinard https://orcid.org/0000-0002-5635-3345

C. Pinte https://orcid.org/0000-0001-5907-5179

L. Ricci https://orcid.org/0000-0001-8123-2943

M. Tazzari https://orcid.org/0000-0003-3590-5814

L. Testi https://orcid.org/0000-0003-1859-3070

D. Wilner https://orcid.org/0000-0003-1526-7587

References

ALMA Partnership, Brogan, C. L., Pérez, L. M., et al. 2015, ApJL, 808, L3

Aly, H., & Lodato, G. 2020, MNRAS, 492, 3306

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

Andrews, S. M., & Williams, J. P. 2007a, ApJ, 671, 1800

Andrews, S. M., & Williams, J. P. 2007b, ApJ, 659, 705

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502

Armitage, P. J. 2015, arXiv:1509.06382

Astropy Collaboration, Price-Whelan, A. M., & Sipőcz, B. M. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bae, J., & Zhu, Z. 2018a, ApJ, 859, 118

Bae, J., & Zhu, Z. 2018b, ApJ, 859, 119

Ballabio, G., Dipierro, G., Veronesi, B., et al. 2018, MNRAS, 477, 2766

Baruteau, C., Crida, A., Paardekooper, S. J., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 667

Baruteau, C., Meru, F., & Paardekooper, S.-J. 2011, MNRAS, 416, 1971

Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362

Bate, M. R., Lodato, G., & Pringle, J. E. 2010, MNRAS, 401, 1505

Bell, K. R., Cassen, P. M., Klahr, H. H., & Henning, T. 1997, ApJ, 486, 372

Benisty, M., Juhász, A., Facchini, S., et al. 2018, A&A, 619, A171

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

Boss, A. P. 2017, ApJ, 836, 53

Cadman, J., Hall, C., Rice, K., Harries, T. J., & Klaassen, P. D. 2020, MNRAS, 498, 4256

Casassus, S., Marino, S., Lyra, W., et al. 2019, MNRAS, 483, 3278

Casassus, S., Wright, C. M., Marino, S., et al. 2015, ApJ, 812, 126

Cazzoletti, P., van Dishoeck, E. F., Pinilla, P., et al. 2018, A&A, 619, A161

Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

Cieza, L. A., Padgett, D. L., Allen, L. E., et al. 2009, ApJL, 696, L84

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

Cuello, N., Dipierro, G., Ment REPL, D., et al. 2019, MNRAS, 483, 4114

Dipierro, G., Laibe, G., Alexander, R., & Hutchison, M. 2018, MNRAS, 479, 4187

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

Dipierro, G., Pinilla, P., Lodato, G., & Testi, L. 2015, MNRAS, 451, 974

Dong, R., Liu, S.-y., Eisner, J., et al. 2018, ApJ, 860, 124

Draine, B. T. 2006, ApJ, 636, 1114

Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

Faccini, S., Juhász, A., & Lodato, G. 2018, MNRAS, 473, 4459

Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Forgan, D. H., Hall, C., Meru, F., & Rice, W. K. M. 2018a, MNRAS, 474, 5036

Forgan, D. H., Ilee, J. D., & Meru, F. 2018b, ApJL, 860, L5

Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

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

Guzmán, V. V., Huang, J., Andrews, S. M., et al. 2018, ApJL, 869, L48

Hall, C., Dong, R., Rice, K., et al. 2019, ApJ, 871, 228

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

Hall, C., Rice, K., Dipierro, G., et al. 2018, MNRAS, 477, 1004

Harsono, D., Alexander, R. D., & Levin, Y. 2011, MNRAS, 413, 423

Hennebelle, P., Lesur, G., & Fromang, S. 2017, A&A, 599, A86

Huang, J., Andrews, S. M., Cleeves, L. I., et al. 2018a, ApJ, 852, 122

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018b, ApJL, 869, L42

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

Hunter, J. D. 2007, CSE, 9, 90

Isella, A., Carpenter, J. M., & Sargent, A. I. 2009, ApJ, 701, 260

Jennings, J., Booth, R. A., Tazzari, M., Rosotti, G. P., & Clarke, C. J. 2020, MNRAS, 495, 3209

Juhász, A., Benisty, M., Pohl, A., et al. 2015, MNRAS, 451, 1147

Juhász, A., & Facchini, S. 2017, MNRAS, 466, 4053

Keppler, M., Penzlin, A., Benisty, M., et al. 2020, A&A, 639, A62

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

Kurtovic, N. T., Pérez, L. M., Benisty, M., et al. 2018, ApJL, 869, L44

Kwon, W., Looney, L. W., Mundy, L. G., Chiang, H.-F., & Kemball, A. J. 2009, ApJ, 696, 841

Launhardt, R., Henning, T., Quirrenbach, A., et al. 2020, A&A, 635, A162

Lee, C.-F., Li, Z.-Y., & Turner, N. J. 2020, NatAs, 4, 142

Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9

Liu, H. B. 2019, ApJL, 877, L22

Lodato, G. 2007, NCimR, 30, 293

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

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17

Luhman, K. L., & Rieke, G. H. 1999, ApJ, 525, 440

Macías, E., Espaillat, C. C., Osorio, M., et al. 2019, ApJ, 881, 159

Marino, S., Perez, S., & Casassus, S. 2015, ApJL, 798, L44

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Meru, F., Juhász, A., Ilee, J. D., et al. 2017, ApJL, 839, L24

Nealon, R., Dipierro, G., Alexander, R., Martin, R. G., & Nixon, C. 2018, MNRAS, 481, 20

Nixon, C. J., & Pringle, J. E. 2010, MNRAS, 403, 1887

Pérez, L. M., Benisty, M., Andrews, S. M., et al. 2018a, ApJL, 869, L50

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

Pérez, S., Casassus, S., & Benítez-Llambay, P. 2018b, MNRAS, 480, L12

Pérez, S., Casassus, S., Hales, A., et al. 2020, ApJL, 889, L24

Perez, S., Dunhill, A., Casassus, S., et al. 2015, ApJL, 811, L5

Pinilla, P., Pérez, L. M., Andrews, S., et al. 2017, ApJ, 839, 99

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967

Pinte, C., Ménard, F., Duchêne, G., et al. 2018b, A&A, 609, A47

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

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

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

Pinte, C., van der Plas, G., Ménard, F., et al. 2019, NatAs, 3, 1109

Pohl, A., Pinilla, P., Benisty, M., et al. 2015, MNRAS, 453, 1768

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

Ratzka, T., Köhler, R., & Leinert, C. 2005, A&A, 437, 611

Ricci, L., Testi, L., Natta, A., & Brooks, K. J. 2010, A&A, 521, A66

Rice, K. 2016, PASA, 33, e012

Rice, K., Lopez, E., Forgan, D., & Biller, B. 2015, MNRAS, 454, 1940

Rice, W. K. M., Armitage, P. J., Bate, M. R., & Bonnell, I. A. 2003, MNRAS, 339, 1025

Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2004, MNRAS, 355, 543

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

Sakai, N., Hanawa, T., Zhang, Y., et al. 2019, Natur, 565, 206

Sierra, A., & Lizano, S. 2020, ApJ, 892, 136

Tazzari, M., Beaujean, F., & Testi, L. 2018, MNRAS, 476, 4527

Teague, R. 2019, JOSS, 4, 1220

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

Teague, R., & Foreman-Mackey, D. 2019, bettermoments: Line-of-sight velocity calculation v1.5, Astrophysics Source Code Library, ascl:1901.009
Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 339

Tomida, K., Machida, M. N., Hosokawa, T., Sakurai, Y., & Lin, C. H. 2017, ApJ, 835, L11

Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J. 2017, ApJ, 845, 44

Tsukagoshi, T., Nomura, H., Muto, T., et al. 2016, ApJL, 829, L35

Tsukamoto, Y., Okuzumi, S., & Kataoka, A. 2017, ApJ, 838, 151

van der Marel, N., Pinilla, P., Tobin, J., et al. 2015, ApJL, 810, L7

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

Vigan, A., Bonavita, M., Biller, B., et al. 2017, A&A, 603, A3

Walsh, C., Daley, C., Facchini, S., & Juhász, A. 2017, A&A, 607, A114

Zhang, S., & Zhu, Z. 2020, MNRAS, 493, 2287

Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJL, 869, L47

Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJL, 877, L18

Zurlo, A., Cieza, L. A., Pérez, S., et al. 2020, MNRAS, 496, 5089