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

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.

Zoom In Zoom Out Reset image size

Table 1. Best-fit Parameters of the Logarithmic Spiral Model

Wavelength	Spiral Arm	Angle Range	R0 (au)	b	Pitch Angle
0.89 mm	NW	−251° to −5°	${249.8}_{-1.1}^{+1.2}$	0.230 ± 0.002	129 ± 01
	SE	−70° to 150°	111.5 ± 0.4	0.247 ± 0.003	139 ± 02
1.3 mm	NW	−250° to −35°	${250.2}_{-2.5}^{+2.3}$	${0.229}_{-0.004}^{+0.003}$	129 ± 02
	SE	−61° to 135°	${115.5}_{-0.9}^{+0.8}$	0.234 ± 0.007	132 ± 04
3.3 mm	NW	−260° to −23°	${249.9}_{-3.8}^{+3.9}$	0.229 ± 0.005	129 ± 03
	SE	−50° to 140°	${113.6}_{-1.1}^{+1.0}$	0.231 ± 0.010	130 ± 06

Download table as:  ASCIITypeset 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

$$\begin{eqnarray}&&r(\theta )={R}_{0}{e}^{b\theta }.\end{eqnarray} \tag{ 1 }$$

Here, θ is the polar angle in radians measured from the north and to the east (left); for reference, see coordinates in Figure 1. R₀, b are free parameters, with R₀ the radius at θ = 0, measured in astronomical units, and b relates to the pitch angle of the spiral arms (ϕ), as $\phi =\arctan (b)$. The uncertainty on the location of each measured maxima is assumed to be the astrometric error ²¹ : ${\rm{\Delta }}p={\rm{60,000}}\cdot {\left(\nu \cdot B\cdot {\rm{S}}/{\rm{N}}\right)}^{-1}$, 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, ∼129 and ∼133, 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; 157 and 164 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.

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 (Cazzoletti 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.

Zoom In Zoom Out Reset image size

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

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 (α) of the spectral energy distribution follows I_ν ∝ ν^α , where I_ν corresponds to the measured intensity, from the image plane, at frequency ν. In the optically thin regime, if α 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 especify a round, 026 × 026 beam (∼31 au). Figure 4 shows in the top panel the azimuthally averaged brightness temperature of the smoothed, equal resolution images.

Zoom In Zoom Out Reset image size

We calculate α as the best-fit slope of a linear model applied to the logarithmic space (α ∝ d lnI_ν /d lnν), using all wavelengths available in this study. From an MCMC fit we obtain a posterior distribution of α, with best-fit and uncertainties computed from the 50th, 16th, and 84th percentile values.

From the azimuthally averaged intensity profiles, we compute an α radial profile shown in the middle panel of Figure 4. We retrieve a disk average α value of 2.6 ± 0.06, with α < 2.0 in the inner ∼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_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_ν ), using Planck's law, and relate it to the measured intensity (I_ν ) through τ_ν , following ${I}_{\nu }(r)={B}_{\nu }({T}_{\mathrm{mid}}(r))(1-\exp (-{\tau }_{\nu }(r)))$. T_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 φ = 0.02 and $\mathrm{log}{L}_{* }/{L}_{\odot }$ = −0.04 ± 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_mid is ∼2 times larger than the measured T_b . In the Rayleigh–Jeans regime, as expected for millimeter emission, T_b = T_mid(1 − exp(−τ)). Then, if T_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σ 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 α values of [2.0, 2.25, 2.4] and for α uncertainties at the level of [0.01, 0.025, 0.05, 0.1] are presented. We observe that contours tracing α = 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 ∼0.1, with uncertainties of ∼0.04. Overall, we observe lower α co-located with the spiral features. At the gap location no particular behavior is observed, except for a decrease in α values inwards from the gap location. As noticed in the radial profile, in the inner disk (<50 au) α reaches values below 2.0.

Zoom In Zoom Out Reset image size

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 ∼125 au corresponding to ∼224°.

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 the azimuthal cut, calculated directly from the top panel intensity profile, as previously done for the spectral index map. In this α 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 α variations. The spatial extent of the variation in the east is resolved by our beam size (026, ∼31 au).

Zoom In Zoom Out Reset image size

The channel maps for both ¹³CO and C¹⁸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 ¹³CO emission from v ≥ 2.88 km s⁻¹, while for C¹⁸O some cloud absorption is present near v = 2.55–2.66 km s⁻¹. We note that the south part of the disk is the brightest. The systemic velocity is determined to be 1.95 km s⁻¹ based on the C¹⁸O channel map analysis using a spectral resolution of 0.05 km s⁻¹ 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 ¹³CO is the most abundant isotopologue, the large-scale emission appears more strongly and along a wider range of velocities than in C¹⁸O. The extended emission is clearly identified between v = 2.55−4.88 km s⁻¹ for ¹³CO, and between v = 2.66–4.21 km s⁻¹ for C¹⁸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 ≲15'' scales; however, the detected emission may extend beyond our ∼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σ, while in ¹³CO, the large-scale emission is seen through the entire FOV, the C¹⁸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.

Zoom In Zoom Out Reset image size

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 ¹³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 05 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¹⁸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, 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.

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.

Zoom In Zoom Out Reset image size

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:

$$\begin{eqnarray}&&z(r)={z}_{0}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{\psi }+{z}_{1}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{\varphi },\end{eqnarray} \tag{ 2 }$$

where r is the radial distance from the star, as measured in the azimuth plane, and the characteristic radius r₀ 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.

Table 2. Height Model Parameters from Channel Analysis

Param.	13CO − W	13CO − E	C18O − W	C18O − E
z0 (au)	${77.6}_{-8.9}^{+8.2}$	30.0 ± 0.2	${44.0}_{-6.2}^{+7.3}$	${32.3}_{-8.5}^{+8.9}$
ψ	${1.76}_{-0.04}^{+0.03}$	1.01 ± 0.01	${2.19}_{-0.16}^{+0.12}$	${1.14}_{-0.09}^{+0.11}$
z1 (au)	$-{38.1}_{-8.1}^{+8.9}$	$-{0.3}_{-0.2}^{+0.1}$	$-{11.8}_{-7.4}^{+6.1}$	$-{15.9}_{-8.9}^{+8.6}$
φ	${2.29}_{-0.04}^{+0.08}$	$-{3.36}_{-0.52}^{+0.54}$	${3.58}_{-0.24}^{+0.46}$	${0.46}_{-0.33}^{+0.15}$

Download table as:  ASCIITypeset image

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 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 ¹³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.

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_obs, together with the systemic velocity for the source, v_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, θ, we can relate the known velocities to v through ${v}_{\mathrm{obs}}={v}_{\mathrm{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 ¹³CO and C¹⁸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.

Zoom In Zoom Out Reset image size

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.49M_⊙ (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:

$$\begin{eqnarray}&&\displaystyle \frac{{v}^{2}}{r}=\displaystyle \frac{{{GM}}_{* }r}{{\left({r}^{2}+{h}^{2}\right)}^{\tfrac{3}{2}}}.\end{eqnarray} \tag{ 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σ 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 ¹³CO measurements, we constrain a stellar mass of M_* = 0.5 ± 0.01 M_⊙, while from the C¹⁸O measurements, we constrain ${M}_{* }={0.46}_{-0.03}^{+0.02}$ M_⊙. Both values are similar between them and compared to the previous estimates (M_* ∼ 0.49M_⊙; 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 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 (∼035–359) 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.

Zoom In Zoom Out Reset image size

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 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 ¹³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¹⁸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 ¹³CO and the C¹⁸O emission traces a lower height layer.

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¹⁸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⁻¹. 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 ¹³CO, along channels 1.55–2.33 km s⁻¹, but in channels 1.66–1.88 km s⁻¹ of C¹⁸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.

Zoom In Zoom Out Reset image size

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¹⁸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¹⁸O emission, imaged with a finer spectral resolution that is available only for the C¹⁸O data. The deviations are most visible at 0.95–1.0 km s⁻¹ 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⁻¹, 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 ¹³CO maps; however, it is expected that perturbations due to the spirals should be more apparent in C¹⁸O than in ¹³CO, as the C¹⁸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.

Zoom In Zoom Out Reset image size

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.

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⁷ 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.5M_⊙ , and R_* = 2.3 R_⊙.

To create the density structure as input into the mcfost calculation, each SPH simulation underwent Voronoi tesselation 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 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 (ϕ_model) and the observational pitch angle (ϕ_obs) are compared considering the difference between the values of each spiral and wavelength and weighing by the error of the observational pitch angle (σ_obs), following a likelihood parameter determined by

$$\begin{eqnarray}&&\displaystyle \sum _{\lambda }\sqrt{{\left[{\left(\displaystyle \frac{{\phi }_{\mathrm{model}}^{\mathrm{NW}}-{\phi }_{\mathrm{obs}}^{\mathrm{NW}}}{{\sigma }_{\mathrm{obsNW}}}\right)}^{2}+{\left(\displaystyle \frac{{\phi }_{\mathrm{model}}^{\mathrm{SE}}-{\phi }_{\mathrm{obs}}^{\mathrm{SE}}}{{\sigma }_{\mathrm{obsSE}}}\right)}^{2}\right]}_{\lambda }}.\end{eqnarray} \tag{ 6 }$$

Table 4. Pitch Angle Values for Simulations

Simulation	NW0.89 mm	SE0.89 mm	NW1.3 mm	SE1.3 mm	NW3.3 mm	SE3.3 mm	Likelihood ϕ
Tapered: 0.7
q = 0.2	1346 ± 008	1496 ± 008	1086 ± 006	1527 ± 008	985 ± 007	137 ± 01	29.31
q = 0.3	1252 ± 006	1154 ± 006	1224 ± 006	1082 ± 005	1209 ± 006	106 ± 006	24.04
q = 0.5	⋯	⋯	⋯	⋯	⋯	⋯	⋯
Tapered: 1.0
q = 0.2	⋯	⋯	⋯	⋯	⋯	⋯	⋯
q = 0.3	1275 ± 006	1142 ± 005	1234 ± 006	1133 ± 005	1308 ± 006	1119 ± 005	21.04
q = 0.5	⋯	⋯	⋯	⋯	⋯	⋯
Power law: 1.3
q = 0.2	139 ± 008	1116 ± 011	1475 ± 007	1040 ± 008	1482 ± 007	1048 ± 011	36.19
q = 0.3	1488 ± 007	1566 ± 007	1553 ± 006	1536 ± 007	1534 ± 007	1592 ± 006	45.32
q = 0.5	1486 ± 007	1808 ± 007	1447 ± 006	1774 ± 008	142 ± 006	1882 ± 007	53.08
Power law: 1.5
q = 0.2	1094 ± 009	734 ± 006	1238 ± 009	58 ± 007	1149 ± 008	824 ± 007	66.08
q = 0.3	1203 ± 006	1195 ± 007	1155 ± 006	104 ± 006	1194 ± 006	1105 ± 007	27.36
q = 0.5	1104 ± 006	1076 ± 005	1062 ± 006	1074 ± 005	1042 ± 005	105 ± 005	46.58

Download table as:  ASCIITypeset image

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

Zoom In Zoom Out Reset image size

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 (∼129 and ∼132, 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.

Zoom In Zoom Out Reset image size

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 1556 ± 006 in the NW spiral and 1276 ± 006 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 157 ± 02 in the NW spiral and 164 ± 02 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.

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 ¹³CO J = 3−2 and C¹⁸O J = 3−2 line maps. The molecule abundances, relative to local H₂ are set to 7 × 10⁻⁷ for ¹³CO (as done in Hall et al. 2020) and 2 × 10⁻⁷ for C¹⁸O (following the estimate from Frerking et al. 1982). The spectral resolution for the kinematic simulations is set to 0.111 km s⁻¹ to match the observations, and we additionally compute C¹⁸O channel maps with 0.05 km s⁻¹ 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 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 (03 compared to 01). 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 α cannot reach values below 2.0. When dust scattering is included, a region of high optical depth can have α < 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 α also depends on other dust properties, such as the grain size distribution (e.g., Testi et al. 2014), values of α < 3 are not expected when the optical depth is low, and even in presence of centimeter-sized grains, α should not attain values below ∼2.5 (Zhu et al. 2019).

From our spectral index profiles α < 2.0 in the inner ∼40 au and α < 2.5 inside ∼200 au, outside this region α 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) ∝ 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) ∝ s^−2.5 and can be explained by either optically thin 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 τ ∼ 0.5 at all wavelengths.

For the DSHARP sample, it has been shown that the optical depths of τ ∼ 0.6 at 1.3 mm in bright rings can be obtained from optically thick regions with a scattering albedo of ω_ν ∼ 0.89 (Zhu et al. 2019). For Elias 2–27, we measure an average τ ∼ 0.45 at 1.3 mm inside 70 au, which could be obtained with a scattering albedo of ω_ν ∼ 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 ∼ 1–5) in the inner disk, which we would be inferring as optically thin 1.3 mm emission (τ_obs ∼ 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 τ ∼ 0.99, which could be obtained with a scattering albedo ω_ν ∼ 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 submillimeter 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.

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 will become asymmetric, as the accreting gas shocks the disk from above or below, along the z plane, as described by cylindrical coordinates (see Figure 8 in Hennebelle et al. 2017). Furthermore, simulations for ongoing infall predict the appearance of spiral structures in the surface, generated by the infall process and shocks, both in 3D simulations (Harsono et al. 2011; Hennebelle et al. 2017) and 2D simulations (Lesur et al. 2015). Infall-triggered spirals may have been observed in the 343 GHz ALMA emission of VLA 1, a Class I source with active envelope infall (Lee et al. 2020). From our observations, we have shown the presence of large-scale emission surrounding the disk (see Figures 7, 21, and 22). Our data lack appropriate uv-coverage to accurately sample the whole FOV, which extends for 20''; this is possibly the cause of the striped pattern seen in the channel maps (Figures 21 and 22). We do not recover velocity gradients connecting the large-scale emission to the disk emission, this could be due to the lack of sensitivity or angular resolution; nevertheless, it opens the option for infall to be ongoing. While this is not expected for a Class II source, it has been proposed that Elias 2–27 could be a very young Class II disk (Tomida et al. 2017).

An azimuthal variation of the temperature in the disk could also explain an azimuthally varying emission layer height, which can be expected when a disk is warped. A warp will affect the disk illumination, depending on the position and characteristic angles of the warp itself (Nixon & Pringle 2010). Warps are generally detected through their shadowing effects in scattered light (Marino et al. 2015; Benisty et al. 2018) and their distinct kinematical signatures (Juhász & Facchini 2017; Walsh et al. 2017; Pérez et al. 2020). Given its high extinction, no scattered-light observations are available for Elias 2–27, so we can only compare our observations with the kinematic predictions for a warp.

In ALMA observations, a misaligned disk can be inferred through the kinematical signature in the velocity map: if a non-misaligned disk is modeled then positive and negative residuals appear, opposite with respect to a "symmetry" axis (e.g., as seen in the velocity field residuals for HD 100546 or HD 143006; Walsh et al. 2017; Pérez et al. 2018a). We do not observe such residuals in our data (see Figure 12), rather we see an "X"-shaped residual in the C¹⁸O emission. Other signatures related to warps are asymmetric illumination, deviations in line profiles, and twisted features in the channel, and integrated emission maps (Juhász & Facchini 2017; Facchini et al. 2018). The gas emission of Elias 2–27 is characterized by a strong illumination asymmetry between the north–south sides of the disk. Additionally, we have shown the presence of "curved" and "wavy" features and deviations in both the integrated velocity maps and velocity cube (see Figures 12 and 13, respectively). These features may relate to an inclined disk, warped such that the disk bends perpendicular to the line of sight (see Case C in Juhász & Facchini 2017). The latter configuration produces channel maps with curved structures and asymmetric illumination, such as observed in our data, while it also shows further structure in the integrated intensity map and asymmetric line profiles (Juhász & Facchini 2017). We may not be sensitive to all these features given our moderate spatial resolution (simulations are done with resolution ∼01; we have ∼03) and the effects of cloud contamination in the system.

If indeed there is a warp, we should be able to roughly trace its location through the temperature variations it will cause. We expect temperature variations to affect the emitting layer's height and should therefore be able to trace these variations in the projected emission. While we do observe an elevation difference in the east and west sides of the disk, apparently 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).

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_J (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 (562), 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 (1188) 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.

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 ∼10M_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.