© The Authors 2023

1. Introduction

Planet formation appears to be a robust and efficient process, occurring around both single and multiple stellar systems (Kostov et al. 2016) in protoplanetary disks. The advent of high resolution imaging facilities has demonstrated that nearly all bright and extended disks show substructures, in particular in the small (micron-sized) and large (millimeter-sized) dust tracers seen through scattered and thermal light, respectively (e.g., Andrews et al. 2018; Long et al. 2018; Rich et al. 2022; Benisty et al. 2022; Bae et al. 2022a). Such high resolution studies applied to gas tracers allow the overall physical conditions in the disk to be probed, such as its temperature structure, its surface height (Rich et al. 2021; Law et al. 2021), and pressure variations (Teague et al. 2018a,b; Rosotti et al. 2020). Studies of the disk density and the velocity structure reveal a great complexity, including localized non-Keplerian features that can be attributed to embedded massive protoplanets (e.g., Pinte et al. 2022; Wölfer et al. 2023). Such perturbations from smooth density and velocity distributions can directly constrain planet formation since they are expected to leave clear signatures on the disk structure (e.g., Perez et al. 2015; Yun et al. 2019). For example, the mapping of spiral wakes (Calcino et al. 2022) and the detections of so-called Doppler flips (change of sign in the non-Keplerian feature; e.g., Casassus & Pérez 2019; Norfolk et al. 2022), meridional flows within dust-depleted gaps (Teague et al. 2019a), and a velocity perturbation associated with a circumplanetary disk candidate (Bae et al. 2022b) enable us to take a closer look at planet–disk interaction processes. While most localized kinematical perturbations are analyzed empirically, statistical methods for quantifying their significance have been developed and led to the detection of localized signatures possibly associated with unseen planets (Izquierdo et al. 2021, 2022). Prime targets to search for protoplanets still embedded in their birth environment are the so-called transition disks. As in PDS70 (Keppler et al. 2019) or AB Aur (Tang et al. 2017), these disks host a dust-depleted cavity that has possibly been cleared by massive companions (Zhu et al. 2011).

In this Letter we focus on RXJ1604.3−2130 A (d = 144.6 pc, 1.24 M_⊙; Gaia Collaboration 2023; Manara et al. 2020, respectively), hereafter J1604, one of the brightest protoplanetary disks of the Upper Scorpius Association in the millimeter regime (Barenfeld et al. 2016), which exhibits a prominent cavity in the dust continuum and CO line emission (Zhang et al. 2014; Dong et al. 2017; van der Marel et al. 2021). J1604 has a stellar companion located at ∼2300 au, itself a binary with separation 13 au (Köhler et al. 2000). The outer disk of J1604 was resolved with the Atacama Large Millimeter/submillimeter Array (ALMA; Mayama et al. 2018) and the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) instrument on the Very Large Telescope (VLT; Pinilla et al. 2015), indicating a nearly face-on geometry. Complementary observations are indicative of a misaligned inner disk with respect to the outer disk. Its variable light curve is that of an irregular dipper (Ansdell et al. 2020), infrared scattered-light observations show the presence of two shadows with variable morphology on timescales possibly shorter than a day (Pinilla et al. 2018), and ALMA ¹²CO (J = 3–2) line observations show deviations from Keplerian rotation in the cavity (Mayama et al. 2018). The position of the scattered-light shadows are suggestive of a large misalignment (∼70–90°). Measurements of the projected rotational velocity (v sin i) indicate that the star is aligned with the inner disk, and thus misaligned with the outer disk (Sicilia-Aguilar et al. 2020).

In this work, we present new ALMA observations of J1604 and focus on the kinematics of the ¹²CO (J = 2–1) line. Section 2 presents the observations and Sects. 3 and 4 our methodology and results, respectively. Section 5 provides a discussion of the results and Sect. 6, our conclusions.

2. Observations, calibration, and imaging

We present new ALMA Band 6 observations (2018.1.01255.S; PI: Benisty) with five executions spread over two years, obtained on 2019 April 4, 2019 July 30 and 31, 2021 April 29, and 2021 September 27. The spectral setup was designed for continuum detection but includes the ¹²CO J = 2–1 line. The data were combined with archival data from program 2015.1.00964.S (PI Oberg; see Table A.1). The data calibration and imaging were performed following the procedure of Andrews et al. (2018), with CASA v.5.6.1 (McMullin et al. 2007), and is detailed in Appendix A. The synthesized beam of the ¹²CO line and dust continuum images is 0.18″ × 0.15″ (102°) and 0.060″ × 0.039″ (−78°), respectively. The rms in a line-free channel was measured to be 1.1 mJy beam⁻¹ (4.3 K) for CO and 10 mJy beam⁻¹ for the dust continuum. Figure 1 shows the dust continuum map (left), which displays a cavity and a bright dust ring peaking at R∼ 0.56″ (∼81 au), and the ¹²CO peak brightness temperature map (center), which indicates a smaller cavity in gas, with a peak at R ∼ 0.39″ (∼56 au). A selection of channel maps can be found in Fig. A.1.

Fig. 1. ALMA observations of J1604. Panel a: 231 GHz dust continuum. The solid black contours are drawn at [5, 15, 25, 35, 45]σ, and the image is plotted with a power-law scaling of γ = 0.6. Panel b: 12CO peak brightness temperature map computed from I0 using the Planck law with black solid contours drawn at [5, 10, 20, 40, 60, 65, 70] σ. Pixels below 5σ are masked. Panel c: peak intensity residuals after subtracting an azimuthally averaged radial profile from the data. We have adjusted the color scale such that residuals smaller than 1σ are white. The beam sizes are shown in the lower-left corner, and the position of the star is marked by a green cross. In panels b and c, we have overlaid the continuum contours in white and black, respectively.

3. Methodology

3.1. Channel maps model

To model the disk line intensity and kinematics, we used the discminer package of Izquierdo et al. (2021). The code uses parametric prescriptions for the line peak intensity, line width, rotational velocity, and disk emission height to produce channel maps; it uses emcee (Foreman-Mackey et al. 2013) to maximize a χ² log-likelihood function of the difference between the model and input intensity for each pixel in a channel map. To prescribe the model intensity, we used a generalized bell kernel function of the disk cylindrical coordinates (R, z):

$$\begin{array}{c}\hfill I_{\mathrm{m}}(R,z;v_{\mathrm{ch}})=I_{\mathrm{p}}(R,z){(1+{\left|\frac{v_{\mathrm{ch}}-v}{L_{\mathrm{w}}(R,z)}\right|}^{2L_{\mathrm{s}}})}^{-1},\end{array}$$(1)

where I_p is the peak intensity, L_w is half the line width at half power, hereafter “the line width,” and L_s is the line slope. The v_ch is the channel velocity at which the intensity is computed and v the observed Keplerian line-of-sight velocity. As the disk is nearly face-on, the code is unable to infer an emission height, and we therefore assumed a flat emission surface. We additionally fixed the inclination, i, of the disk to the one inferred from the dust continuum (i = 6.0°; Dong et al. 2017) to break the degeneracy of M_⋆ ⋅ sin i. The fitting procedure and the Markov chain Monte Carlo (MCMC) search are explained in detail in Appendix B, where the functional form of each model parameter together with its best-fit parameters are summarized in Table B.1. We compare selected channel maps to the best-fit model using these parameters in Fig. A.1.

3.2. Moment maps

The moment maps were computed with bettermoments (Teague & Foreman-Mackey 2018). Since the ¹²CO line emission is optically thick, we fitted the following line profile to both the data and model channel maps:

$$\begin{array}{c}\hfill I(v)=I_0·\frac{1-exp(-\tau \left(v\right))}{1-exp(-{\tau }_0)}\text{with}\tau ={\tau }_{0}exp\left(\frac{-{(v-v_0)}^2}{\mathrm{\Delta }{V}^2}\right),\end{array}$$(2)

where I₀ is the peak intensity of the line and the optical depth, τ (v), varies like a Gaussian; v₀ is the line centroid, τ₀ the peak optical depth, and ΔV the width of the line – where the full width at half maximum (FWHM) = $2\sqrt{\ln 2}\mathrm{\Delta }V$ – as used in Teague et al. (2022). In Fig. 1b we show I₀ for ¹²CO in units of brightness temperature. The corresponding v₀ maps for the data and model are displayed in panels a and b of Fig. 2, respectively. The moment maps for ΔV and τ₀, as well the error of the line centroid fitting, δv₀, can be found in Fig. A.4.

Fig. 2. Line-of-sight velocity maps for v0 data (a) and the discminer model, vmod (b). Panel c: velocity residual map after subtracting vmod from v0. The dust continuum is overlaid in solid contours with equal levels, as in Fig. 1. The innermost region was masked during the fit by one beam size in radius, shown as the gray shaded ellipse. The insets in panels a and c zoom into the innermost region of the disk to highlight the non-Keplerian velocities. Contours are drawn at vsys = (4.62 ± 0.60) km s−1 in steps of 0.1 km s−1 and from −60 to 60 m s−1 in steps of 10 m s−1, respectively. All maps show the synthesized beam for CO (black) and the continuum (white) in the lower-left corner. The beams are masked where the CO peak intensity falls below the 5σ level for panel a and where R > Rout for the rest.

4. Results

4.1. Dust and gas radial and azimuthal brightness profiles

Figure 1 shows the 1.3 mm dust continuum together with the peak brightness temperature map, I₀, of the ¹²CO (J = 2 − 1) line emission. Both dust and gas tracers show a cavity, and the ¹²CO (J = 2 − 1) line emission extends inward of the dust continuum, as expected if the continuum ring results from dust trapping (e.g., Facchini et al. 2018b, see our Fig. A.2). We note that the ¹²CO cavity appears asymmetric with respect to the position of the star and that we observe a gap-like feature in the ¹²CO peak intensity map at R ∼ 1.2″, apparent as a plateau of I₀ ≈ 31 mJy beam⁻¹ stretching over ΔR ≈ 0.1″ (Fig. A.2). Interestingly, the disk shows significant azimuthal intensity variations (34% at R = 0.56″; 19% at 0.39″ for continuum and gas, respectively; see also Fig. A.3). Figure 1c shows the residuals obtained after subtracting an azimuthally averaged radial profile from the ¹²CO peak brightness temperature map. Azimuthal variations are clearly apparent within the dust cavity, with residual values of > 10σ. The fainter regions, distributed along the east-west direction, are broadly aligned with fainter regions seen in the continuum (see the contours of Figs. 1b and A.3) and with the shadows reported in scattered light (Pinilla et al. 2018; Kurtovic et al., in prep.).

4.2. Kinematical features

4.2.1. Localized velocity residuals

The centroid residual map in Fig. 2c shows a prominent, localized non-Keplerian velocity feature of δv ≈ 60 m s⁻¹, between ∼0.35″ and 0.55″ (i.e., 50–80 au), which is at the edge of the dust continuum ring and oriented at PA ≈ (270 ± 15)°. To assess its significance, we followed the variance peak method from Izquierdo et al. (2021). First, the centroid velocity residuals were folded and subtracted along the disk minor axis to remove axisymmetric features. Second, a 2D scan was performed to search for peak velocity residuals and obtain their locations in the folded map. Using these detected points, a K-means clustering algorithm then searched for coherent velocity perturbations within predefined radial and azimuthal bins (MacQueen 1967; Pedregosa et al. 2011). We considered seven radial and ten azimuthal bins, which corresponds to a width of roughly one beam size, to identify clusters. The algorithm then subdivided the input residual points such that the center of each cluster is the closest to all points in the cluster, by iteratively minimizing the sum of squared distances from the data points to the center of the cluster. This led to irregularly spaced bin boundaries since the cluster centers are near the densest accumulations of points.

In Fig. 3 we show the folded velocity residual map together with the detected peak velocity residuals (gray points). The location of the detected peak velocity residuals in azimuth and radius, within identified clusters, can be found in Fig. B.1. Clusters with high significance (those with peak velocity residuals larger than three times the variance in other clusters) are located within one radial and azimuthal bin shown in Fig. 3. Taking the centers of the selected clusters allowed us to identify a localized perturbation at 0.28″ ± 0.07″ (R = 41 ± 10 au) and PA = 280° ±2°. The reported errors are the standard deviation of the peak residual point (R, ϕ) locations within the selected clusters. The detection yields a cluster significance of 5.4σ_ϕ in azimuth and 5.3σ_R in radius, where σ represents the standard deviation of background cluster variances with a mean of σ_ϕ = 0.034 km² s⁻² and σ_R = 0.018 km² s⁻² (see black crosses in Fig. B.1). We note that a localized signature is robustly detected regardless of the amount of clusters defined, which we tested using 7–12 azimuthal or 5–9 radial clusters. We reported the clusters associated with the highest significance. Additionally, we note that there are other detections with lower significance at 0.65″ (94 au). This means that the radial extent of the prominent perturbation is roughly 0.40″ (58 au) and the global peak of the folded velocity residuals is at 0.39″ (56 au; as can be seen in the middle panel of Fig. B.1). This analysis confirms the presence of a significant localized non-Keplerian feature, as identified visually in Fig. 2c, within the continuum ring.

Fig. 3. Folded velocity residuals (left) and detected clusters of peak velocities (right) in the disk reference frame. The green shaded annulus and wedge in the right plot mark the significant clusters in radius and azimuth, respectively. The position of the localized velocity perturbation inferred from these clusters is marked with a magenta point with error bars. The gray region (one beam size in radius) indicates the masked area, and the black dashed lines the FWHM of the dust ring.

4.2.2. Spiral feature

Figure 2c also shows an extended arc-like positive residual feature, beyond the dust continuum emission and covering nearly 300° in azimuth, that is more evident in the polar deprojection of the velocity residual map (Fig. 4). To assess if this feature is a coherent structure, we used the FilFinder package (Koch & Rosolowsky 2015) implemented in discminer between 0.30″ and 1.25″ (43–180 au; see Fig. A.5). As indicated by the coherent filaments, the strong localized positive velocity residual discussed in Sect. 4.2.1 seems to be the starting point of a spiral tracing outward up to roughly 1.1″ (159 au). In Fig. 4 we overlay an Archimedean (linear) spiral, prescribed by r_spiral = a + b ϕ_spiral, using {a, b}={0.48, 0.12}. Computing the pitch angles tan(β) = − (dr/dϕ)/r, we obtain values ranging from 14° to 6° over the spiral extent.

Fig. 4. Polar projection of the velocity residual map. The solid black line shows a linear spiral trace. The gray region indicates the masked area and the black dashed lines, the FWHM of the dust ring. The y-axis extends farther than 360° to enhance the visibility of the spiral.

4.2.3. A possible warp in the ¹²CO cavity

An additional feature clearly evident from the velocity maps is the highly perturbed inner disk regions. As seen in the inset of the v₀ line centroid map in Fig. 2a, the iso-velocity lines show strong bending in the inner region (∼3 beam sizes in diameter), indicative of non-Keplerian velocities. It likely traces a warp or a misaligned inner disk, as reported in Mayama et al. (2018), Pinilla et al. (2018), Sicilia-Aguilar et al. (2020), and Ansdell et al. (2020) to explain the scattered-light shadows and variable photometry of J1604. Higher angular resolution deep gas observations are, however, needed to assess its morphology and kinematics.

4.3. Deprojected velocity components

To understand the contributions from v_ϕ, v_r, v_z, we produced three centroid residual maps for each velocity component, after deprojection (Eq. (C.1)), assuming that “all” the velocities are azimuthal, radial, or vertical (see Fig. C.2; Teague et al. 2022). The localized residual feature at the edge of the dust ring appears to trace variations in the vertical motions, v_z, or in the rotational motions, v_ϕ, or a combination of both. Radial perturbations can be ruled out since the feature is located close to the disk redshifted major axis (PA_disk = 258°), where v_r, proj ≈ 0. Assuming purely rotational velocities, this corresponds to perturbations as high as δv_ϕ ≈ 600 m s⁻¹ (∼0.4 × v_kepler) due to the low disk inclination.

As seen in Fig. 2c, the spiral-like velocity residual feature does not change sign around the disk major or minor axes, which would occur for rotational, v_ϕ, or radial, v_r, velocity perturbations, respectively (see Eq. (C.1)). We are likely seeing vertical perturbations, which we are most sensitive to in a nearly face-on disk.

Figure C.1 shows the deprojected and azimuthally averaged radial profiles of each velocity component determined with eddy. For v_ϕ, we observe super-Keplerian rotation from R ∼ 0.35″–0.70″ (51–101 au), peaking at 0.45″ (65 au), right beyond the dust continuum. The rotational velocities then sharply drop to sub-Keplerian in the inner disk regions. However, we stress that the azimuthally averaged velocities at the radial location of the strong localized perturbation (R ∼ 0.3 − 0.6″,43–87 au) are likely affected by the feature. We tentatively observe radial inflow inward of the CO intensity peak but with very large uncertainties on v_r. Finally, we mostly detect downward vertical motion of the disk within R ∼ 1.25″ (181 au).

5. Discussion

5.1. Origin of v₀ residuals

In this Letter we report the detection of two main non-Keplerian features, in addition to highly perturbed gas velocities in the gas cavity: (1) a localized positive residual near the edge of the dust ring and (2) an extended spiral-like feature, possibly starting from (1). A variety of velocity residual features were detected in other systems, with a diverse range of inclinations (e.g.,Wölfer et al. 2023). In the case of TW Hya, a similarly face-on disk, the detected perturbations are ∼40 m s⁻¹ (Teague et al. 2022), which is lower than what is derived for (1), and in our case can account for 40% of the local Keplerian velocity assuming that the perturbation is purely due to rotational velocities. This is also larger in the magnitude of deviation than the Doppler flip reported in the HD 100546 transition disk (Casassus et al. 2022). These velocity residuals are often interpreted as tracing planet–disk interactions from massive companions (Pinte et al. 2022) that potentially carve out gaps. It is thus worth noting that our inferred planet location (R = 41 ± 10 au) is close to the gap location in ¹³CO (J = 2 − 1) at 37 au reported in van der Marel et al. (2021). Comparison with simulations (Rabago & Zhu 2021; Izquierdo et al. 2022) or semi-analytical prescriptions (Bollati et al. 2021) allowed us to estimate a possible planet mass from the velocity deviations. To this end, we considered Eq. (14) from Yun et al. (2019), which relates the change in rotational velocity, δv_ϕ, to the planet mass, M_p, through 2D hydrodynamic simulations.

Since δ_V is the difference between the super- and sub-Keplerian peak, we considered the peak of the super-Keplerian motion, (v_ϕ/v_kep)/v_kep (see Fig. C.1), as a lower limit for the dimensionless amplitude of the perturbed rotational velocity (${\delta }_{\text{V}}^{\text{min}}$ = 0.06), and its double (${\delta }_{\text{V}}^{\text{max}}=0.12$) as an upper limit. Assuming (H/R)_p = 0.1 at the planet location (R = 41 au), we estimate the planet mass to roughly be in the range M_p ≈ (1.6 − 2.9)M_jup (α/10⁻³)^0.5.

The extended spiral-like feature appears to be related to the significant localized velocity residual. Due to its low pitch angles and the consistent positive velocity residuals, we speculate that the spiral is caused by buoyancy resonances excited by a massive planet located within the dust ring. Indeed, in contrast with Lindblad spirals, buoyancy spirals are shown to exhibit a tightly wound morphology with predominantly vertical motions (Bae et al. 2021). Such a spiral has also been suggested in TW Hya (Teague et al. 2019b), whose radial extent is similar to the one reported here. Interestingly, Wölfer et al. (2023) report a tentative arc feature in J1604, at R ∼ 1.0″ and ranging from PA ≈ 160 − 200°, probed by the kinematics of the ¹²CO (J = 3 − 2) line emission, that partly coincides with the spiral-like feature that we detect. Additional observations in optically thin tracers would be very useful to assess its nature.

5.2. Kinematic perturbations due to shadows

The localized residual feature seems to roughly align with the shadows detected in scattered light. Comparing Figs. 1c and 2c, positive and negative v₀ residuals broadly align with cold and hot regions in the brightness temperature of ¹²CO, respectively (see also Fig. A.6). In particular, the orientation of the significant localized velocity perturbation coincides with the western shadow. Such a shadow can cool down the disk material and possibly induce a local drop in pressure support and therefore impact the gas velocity. In this section we estimate whether the detected velocity perturbations could be caused by azimuthal variations in temperature. We related the azimuthal change in temperature, ΔT_ϕ, to variations in rotational velocity, Δv_ϕ, by solving for the Navier–Stokes equation in cylindrical coordinates. Following the derivation in Appendix D, we obtain

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{v}_{\varphi }}{v_{\mathrm{kep}}}\approx {\left(\frac{H}{R}\right)}^2\frac{\mathrm{\Delta }{T}_{\varphi }}{T},\end{array}$$(3)

with H/R the disk aspect ratio, T the ¹²CO brightness temperature tracing the gas temperature (as ¹²CO is optically thick), and v_kep the Keplerian velocity. Evaluating this equation for a large H/R = 0.2 along a radially averaged annulus centered at R = (0.39 ± 0.07)″, where we experience the strongest azimuthal changes in the ¹²CO brightness temperature of up to ΔT ≈ 15 K over T = 60 K, we estimate the change in azimuthal velocity to be a mere δv_ϕ/v_kep ≈ 1%. Hence, the shadows cannot be solely responsible for the localized velocity residual feature. In addition, as the shadows are nearly symmetric, we would expect a similar feature opposite to its east.

Montesinos et al. (2016) investigated the development of spirals due to pressure gradients caused by temperature differences between obscured and illuminated regions. In their simulations, symmetric shadows always form two-armed spirals. However, they only develop for massive (0.25 M_⋆) and/or strongly illuminated disks (100 L_⊙), which does not seem to be the case for J1604 (∼0.02 M_⊙ and 0.7 L_⊙, Manara et al. 2020), where we also only observe one spiral. Additionally, we would also expect such spiral features to appear in the brightness temperature residuals (Fig. 1c). It is therefore unlikely that shadows are responsible for the extended spiral-like velocity residual feature.

5.3. A warped or misaligned inner disk?

The bending of the iso-velocity curves that we observe in the inset of Fig. 2a is reminiscent of a warped or (broken) misaligned inner disk (Juhász & Facchini 2017; Facchini et al. 2018a). However, as the inner disk is unresolved in our observations, the warp morphology cannot be derived. We note that radial inflows are also degenerate in appearance with warps, as shown by Rosenfeld et al. (2014), and that our observations do not allow us to be conclusive regarding the origin of the disturbed kinematics in the innermost disk. We attempted to infer if the position angle or inclination vary with radius by fitting the innermost disk only (R ≤ 0.5″) with eddy and considering a fixed stellar mass, but we did not find any significant variations relative to our best-fit values. We therefore constrain the warp to be confined within one beam size (∼0.18″, i.e., 26 au) from the center. We note that we obtain a 5% higher dynamical mass of the system when fitting for M_⋆ while masking the innermost beam size in radius, an effect predicted by hydrodynamical simulations of warps (Young et al. 2022).

While our observations cannot provide a full picture of the system due to a limited angular resolution, the very low mass accretion rate and near-infrared excess (Sicilia-Aguilar et al. 2020), as well as the gas cavity in ¹²CO with non-Keplerian velocities, suggest the presence of an additional, very massive (possibly stellar) companion within the inner ∼0.25″ (∼35 au). Such a companion would need to be on an inclined (nearly polar) orbit to misalign the inner disk (Zhu 2019), which would then lead to the shadows (Nealon et al. 2019) and variable extinction events seen in the light curves (Ansdell et al. 2020; Sicilia-Aguilar et al. 2020). It would not, however, explain the strong localized velocity residual feature that we report, which we speculate traces a planetary-mass object located at the edge of the dust continuum. Detailed modeling of the system is thus needed to assess the need for an additional companion. An interesting comparison is the CS Cha spectro-binary system (separation of ∼7 au), which shows a similar dust continuum and gas emission at a similarly low inclination but no departure from Keplerian rotation in the ¹²CO kinematics (Kurtovic et al. 2022).

6. Conclusions

In this Letter we present new ALMA observations of the 1.3 mm dust continuum and the ¹²CO (J = 2 − 1) line emission from the transition disk around RXJ1604.3–2130 A. The dust continuum shows a large cavity enclosing a smaller ¹²CO cavity. Azimuthal brightness variations in the ¹²CO line and dust continuum are broadly aligned with shadows detected in scattered light (Pinilla et al. 2018). Using the discminer package (Izquierdo et al. 2021), we modeled the channel-by-channel line emission and calculated the line-of-sight velocity maps. We report the detection of a coherent, localized non-Keplerian feature at R = (41 ± 10) au and PA = 280° ±2°, that is, within the continuum ring. While broadly aligned with the scattered-light shadows, the localized non-Keplerian feature cannot be due to changes in temperature. Instead, we interpret the kinematical perturbation as tracing the presence of a massive companion of M_p ≈ (1.6 − 2.9) M_jup. We also detect a tightly wound spiral that extends over 300° in azimuth, possibly connected to the localized feature and caused by buoyancy resonances driven by planet–disk interactions. Bending of the iso-velocity contours within the gas cavity indicates a highly perturbed inner region, possibly related to the presence of a misaligned inner disk. However, as the putative planet at ∼41 au cannot explain the gas cavity, the low accretion rate, or the misaligned inner disk, we speculate that another massive companion, likely on an inclined orbit, shapes the inner ∼0.25″(∼35 au).

Acknowledgments

We would like to thank the anonymous referee for the constructive feedback, as well as Clement Baruteau, Kees Dullemond, Guillaume Laibe and Andrew Winter for helpful discussions. This Letter makes use of the following ALMA data: ADS/JAO.ALMA#2017.A.01255.S and ADS/JAO.ALMA#2015.1.00964.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. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (PROTOPLANETS, grant agreement No. 101002188). L.P. gratefully acknowledges support by the ANID BASAL projects ACE210002 and FB210003, and by ANID, – Millennium Science Initiative Program – NCN19_171. Software: CARTA (Comrie et al. 2021), CASA (McMullin et al. 2007), Discminer (Izquierdo et al. 2021), Eddy (Teague 2019a), FilFinder (Koch & Rosolowsky 2015), GoFish (Teague 2019b), Matplotlib (Hunter 2007), Numpy (van der Walt et al. 2011), Scipy (Virtanen et al. 2020).

References

Appendix A: Observations, calibration, and imaging

To self-calibrate our observations, we proceeded as follows. We first flagged the channels containing the line to produce a continuum data set. We centered the individual execution blocks (EBs) by fitting the continuum visibilities with a ring model, allowing for a different center and amplitude. This enabled us to recover the phase shift and amplitude re-scaling and apply them to the EBs before combining them. To determine a good initial model for the self-calibration, we used multi-scale cleaning with the tclean task using a threshold of two times the rms noise level of the image. Using the tasks gaincal and applycal, we corrected for phase offsets between spectral windows, and between polarizations, considering a solution interval of the scan length (solint=inf). Executions obtained in 2019 were concatenated and self-calibrated together, and a similar procedure was applied to those obtained in 2021. In addition to the first round of self-calibration, two additional iterations of phase self-calibration were done, with solution intervals of 300s and 180s for the 2019 data; only one round was done for the 2021 data, with a solution interval of 360s. For both data sets, a round of amplitude self-calibration was applied with solint=inf. The solutions were then applied to the gas data. While these two epochs will be analyzed separately for the continuum in a forthcoming paper (Kurtovic et al.), to analyze the gas data, we concatenated them after checking that the data do not show significant variations between the two epochs. We imaged the resulting visibilities with the tclean task using the multi-scale CLEAN algorithm with scales of 0, 1, 3, and 6 times the beam FWHM and an elliptic CLEAN mask encompassing the disk emission. The ¹²CO (2-1) molecular line observations are imaged with a robust value of 1.0 and a channel width of 0.1 km s⁻¹ and masked by the 4.0 σ threshold. The data were tapered to 0.″1, and we used the “JvM correction” (Jorsater & van Moorsel 1995; Czekala et al. 2021).

In Fig. A.7 we show a comparison of velocity residual maps for additional cubes, imaged using different imaging parameters, to assess the robustness of our detections. We compare residual maps for the same cube as used in the main text, but without JvM correction, and for a cube imaged with a different tapering (0.15″ instead of 0.10″). Best-fit Keplerian models were subtracted from each of the cubes. As is evident from the comparison of the residual maps, the detection of the reported non-Keplerian features are robust irrespective of the imaging parameters. However, the detailed morphology of the velocity residual peaks changes with imaging parameters and, as a consequence, the value of the inferred planet location from the discminer analysis. We find that the best-fit discminer models to the non-JvM-corrected cubes are similar within 3% with respect to the best-fit parameters listed in Appendix B, with the exception of the line slope and some of the peak intensity parameters, which vary by up to 15%. While estimating the systematics due to imaging parameters is beyond the scope of this Letter, Fig. A.7 provides evidence that the detection of non-Keplerian features is robust. We measure the rms in a line-free channel to be 2.6 mJy beam⁻¹ and 2.9 mJy beam⁻¹ for the non-JvM-corrected cubes with a taper of 0.10 and 0.15, respectively.

Fig. A.1. Gallery of selected channel maps. Panels show the 12CO data channel maps (top row) and best-fit model channel maps (middle row), together with intensity residuals in kelvins for each channel (bottom row). In the bottom row the colorbar has been adjusted such that residuals smaller than 1σ are white. The beam size is depicted in the lower-left corner of each channel. For reference, the best-fit systemic velocity was found to be vsys = 4.62 km s−1 and the channel spacing is 100 m s−1.

Fig. A.2. Radial profile of the surface brightness for different tracers. Profiles are normalized to the peak of the emission for the 231 GHz continuum and the CO peak flux (both for the data and discminer model) and for the SPHERE scattered-light observation. Shaded regions show the standard deviation of each azimuthal average. The lines in the lower-right corner show the major beam size (resolution) for each profile in the corresponding color.

Fig. A.3. Azimuthal profiles of the surface brightness, normalized to the peak of the emission. Profiles were extracted at an annulus with a width of approximately one corresponding beam size centered at 0.56″ and 0.39″ for the 231 GHz continuum and the CO peak flux, which both show significant azimuthal intensity variations, of 34% and 19%, respectively. Shaded regions show the standard deviation of each radial average.

Fig. A.4. Additional moment maps of the centroid fitting. Panels show the line width, ΔV (left), the peak optical depth, τ0 (center), and the error of the centroid fitting, δv0 (right). Note that for τ0 < 1 one can assume the line profile to be well presented by a Gaussian, while for τ0 > 5 the line profile has a saturated core, i.e., a very flat top (see Eq. 2). The beam size is depicted in the lower-left corner, and only regions where I0 > 5σ with σ = 1.1 mJy beam−1 are shown.

Fig. A.5. Polar map of the velocity residuals. Same as Fig. 4, but now overlaid with filamentary structures found by FilFinder. The overplotted red and blue lines are the medial axes of the filamentary structures found by the algorithm. To trace the apparent spiral in the residuals, we restricted the algorithm to search for filaments in the radial locations r = [0.3, 1.25]″. For the filamentary detection, we assumed a smoothing size of one synthesized beam size and a minimum size of 500 pixels for a filament to be considered.

Fig. A.6. Polar contour map of the centroid residuals. The upper panel shows residuals inside the cavity, the middle panel residuals between the cavity and the outer edge of dust ring, and the lower panel the residuals in the outer disk. The radial spacing between each contour is ∼1.8 au, and the opacity of the lines increases with radius. We would like to emphasize that the bump between PA≈ − (110 − 70)° in the middle panel includes most of the residual points we detect with the peak variance method of discminer, which can be readily seen in the left panel of Fig. B.1.

Fig. A.7. Comparison of the velocity residual maps for different imaging parameters. The left panel corresponds to the cube used in the main text, while the middle panel corresponds to the same cube without JvM correction. The right panel shows the residual maps for a non-JvM-corrected cube with a different taper (0.15″). In all panels the dust continuum is overlaid in solid contours with equal levels, as in Fig. 1. Best-fit Keplerian models were subtracted from each of the cubes. The detection of the non-Keplerian features is quite robust, irrespective of the imaging procedure.

Appendix B: Model best-fit parameters

For the initial emcee run, we used literature values for the position angle and stellar mass (PA = 260°, M_⋆ = 1.24 M_⊙; Dong et al. 2017; Manara et al. 2020, respectively). The initial values of the other parameters were found by comparing the overall morphology between the data and a prototype model. We performed the MCMC search with 150 walkers that evolved for 2000 steps for an initial burn-in stage. We proceeded in two steps. First, we masked the disk region inward of the dust continuum and only fitted the outer disk (R > 90 au) to get a robust estimate of the stellar mass and avoid the code being confusing by the strongly non-Keplerian velocity features in the inner regions. In this run, we interestingly find a strong offset from the disk center in the x direction of −8.0 au. In a second step, we fixed the stellar mass and now fitted for the whole disk, masking an inner region corresponding to one major beam size (26 au) in radius, where effects of beam smearing are strongest. We ran 150 walkers for 20000 steps until reaching convergence, which presented as a nearly normal distribution of the walkers. The variance and median of the parameter walkers remain effectively unchanged after ∼7000 steps. The best-fit parameters are the median of the posterior distributions, and given errors are the 16 and 84 percentiles in the last 5000 steps of the 20000-step run, summarized in Table B.1.

Fig. B.1. Location of the folded peak velocity residuals. The detected points are shown in azimuth (left) and radius (middle) and were obtained using the peak variance method. Colors indicate the seven different radial clusters specified, where blue peak residual points are within a detected significant radial cluster. The black crosses are the velocity variances of the clusters plotted at the (R, ϕ) location of each cluster center. The centers of the accepted clusters (those with peak velocity residuals larger than three times the variance in other clusters) in radius and azimuth are marked with vertical black lines in both panels. The right-hand plot shows the normal distribution of the peak residual points in a histogram. Note that outliers of the distribution are related to the localized perturbation. The maximum value of all peak folded centroid residuals is at 0.39″ (57 au), its mean value is 39 m/s, and 1σv = 20 m/s (not to be mistaken with the cluster variances).

Appendix C: Decomposition and deprojection of velocity components

To determine the rotation curves of each velocity component, we used the code eddy (Teague 2019a). We followed the method presented in Teague et al. (2018b), which uses a Gaussian process to determine the azimuthal, v_ϕ, and radial, v_r, velocity components along a given annulus. To this end, we divided the disk into concentric annuli with a radial width of one-fourth of the synthesized beam (∼0.05″), which ranges from 0.18″ to 1.85″, and extracted the velocities over 20 iterations to minimize their standard deviation. To obtain the vertical velocity component, v_z, we used the measured azimuthally averaged profiles of v_ϕ and v_r and extended them to produce 2D maps, considering the projection of the following components along the line of sight:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill v_{\varphi ,\mathrm{proj}}& =v_{\varphi }cos(\varphi )sin(|i|),\hfill \\ \hfill v_{r,\mathrm{proj}}& =v_{r}sin(\varphi )sin(i),\hfill \\ \hfill v_{z,\mathrm{proj}}& =-v_{z}cos(i),\hfill \end{array}\end{array}$$(C.1)

where ϕ is the polar angle in the disk frame (such that ϕ= 0 corresponds to the redshifted major axis) and i the inclination of the disk. In the case of J1604, the disk rotates clockwise, which corresponds to a negative inclination in the above definition. We then subtracted these maps together with the systemic velocity, v_sys, from the line-of-sight velocity, v₀, map (Fig. 2,a) to obtain a map of the vertical velocity component, v_z,  proj:

$$\begin{array}{c}\hfill v_{z,\mathrm{proj}}=v_0-v_{\mathrm{sys}}-v_{\varphi ,\mathrm{proj}}-v_{r,\mathrm{proj}}.\end{array}$$(C.2)

The radial profile of v_z was obtained by deprojecting and azimuthally averaging its 2D velocity map. The radial profiles of the deprojected velocity components can be found in Fig. C.1.

Fig. C.1. Azimuthally averaged and deprojected azimuthal, radial, and vertical velocity components. The radial width of each annulus is one-fourth of the synthesized beam size. The error bars are given by the standard deviation for each velocity component averaged over the 20 iterations used.

Fig. C.2. J1604 deprojected velocity components. It is assumed that all velocities are azimuthal (left column), radial (central column), or vertical (right column). For the azimuthal and radial components, wedges along the minor and major axis have been masked as the observations are insensitive to these components (see Eq. C.1). In each panel the synthesized beam is shown in the lower-left corner.

Appendix D: Derivation of Δv_ϕ as a function of azimuthal temperature variations, ΔT

To relate the change in the brightness temperature, ΔT, of ¹²CO to variations in the rotational velocity, v_ϕ, we solved the Navier-Stokes equation in cylindrical coordinates in the ϕ direction:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\rho }_g\frac{v_{\varphi }}{R}\frac{\partial v_{\varphi }}{\partial \varphi }=& \frac{1}{R}\frac{\partial p}{\partial \varphi }+\mu \left(\frac{\partial }{\partial R}\left(\frac{1}{R}\frac{\partial }{\partial R}(R{v}_{\varphi })\right)\frac{1}{R^2}\frac{{\partial }^2{v}_{\varphi }}{\partial {\varphi }^2}\right),\hfill \end{array}\end{array}$$(D.1)

where R is the cylindrical radius, ρ_g the gas density, and μ the mean molecular weight. In a first step, we assumed that the radial variations within the chosen annulus are negligible and inserted the disk gas pressure in the vertically isothermal assumption: $p={\rho }_g{c}_s^2={\rho }_g\frac{k_{\mathrm{B}}T}{\mu m_p}$, where k_B is the Boltzmann constant and m_p the proton mass. In a second step, we further assumed the gas density to be constant along the annulus ρ_g = const. and re-arranged the equation:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\rho }_g\frac{v_{\varphi }}{R}\frac{\partial v_{\varphi }}{\partial \varphi }=& \frac{1}{R}\frac{\partial }{\partial \varphi }\left({\rho }_g\frac{k_{\mathrm{B}}T}{\mu m_p}\right)+\frac{\mu }{R^2}\frac{{\partial }^2{v}_{\varphi }}{\partial {\varphi }^2}\hfill \\ \hfill \frac{\partial v_{\varphi }}{\partial \varphi }-\frac{\mu }{v_{\varphi }R{\rho }_g}\frac{{\partial }^2{v}_{\varphi }}{\partial {\varphi }^2}=& \frac{k_{\mathrm{B}}}{v_{\varphi }\mu m_p}\frac{\partial T}{\partial \varphi }\hfill \\ \hfill \frac{\partial v_{\varphi }}{\partial \varphi }-\frac{0.4\mu }{v_{\varphi }\sqrt{2\pi }{\mathrm{\Sigma }}_g}\frac{{\partial }^2{v}_{\varphi }}{\partial {\varphi }^2}=& \frac{k_{\mathrm{B}}}{v_{\varphi }\mu m_p}\frac{\partial T}{\partial \varphi }\hfill \end{array}.\end{array}$$(D.2)

In the last step, we inserted the gas midplane density, ${\rho }_g={\mathrm{\Sigma }}_g/(\sqrt{2\pi }H)={\mathrm{\Sigma }}_g/(\sqrt{2\pi }0.2R),$ assuming a disk aspect ratio of H/R = 0.2. Assuming that v_ϕ ≈ v_kep and Σ_g ≈ 1g/cm² (see Fig. 3 of Dong et al. 2017) at the location of the annulus at R ∼ 0.4″ (58 au), we can assess the order of magnitude of the second term on the left-hand side of the equation, which is only on the order of 10⁻⁶. Therefore, we neglected the second-order derivative and identified the sound speed, c_s:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{\Delta }{v}_{\varphi }\approx & \frac{k_{\mathrm{B}}}{v_{\mathrm{kep}}\mu m_p}\frac{T}{T}\mathrm{\Delta }{T}_{\varphi }\approx \frac{c_s^2}{v_{\mathrm{kep}}}\frac{\mathrm{\Delta }{T}_{\varphi }}{T}|÷v_{\mathrm{kep}}\hfill \\ \hfill \frac{\mathrm{\Delta }{v}_{\varphi }}{v_{\mathrm{kep}}}\approx & {\left(\frac{c_s}{v_{\mathrm{kep}}}\right)}^2\frac{\mathrm{\Delta }{T}_{\varphi }}{T}\approx {\left(\frac{H}{R}\right)}^2\frac{\mathrm{\Delta }{T}_{\varphi }}{T}\hfill \end{array}.\end{array}$$(D.3)

Equation D.3 now connects the fractional azimuthal temperature variation, ΔT_ϕ/T, to the rotational velocity deviation relative to Keplerian, Δv_ϕ/v_kep.

All Tables

All Figures

Fig. 1. ALMA observations of J1604. Panel a: 231 GHz dust continuum. The solid black contours are drawn at [5, 15, 25, 35, 45]σ, and the image is plotted with a power-law scaling of γ = 0.6. Panel b: 12CO peak brightness temperature map computed from I0 using the Planck law with black solid contours drawn at [5, 10, 20, 40, 60, 65, 70] σ. Pixels below 5σ are masked. Panel c: peak intensity residuals after subtracting an azimuthally averaged radial profile from the data. We have adjusted the color scale such that residuals smaller than 1σ are white. The beam sizes are shown in the lower-left corner, and the position of the star is marked by a green cross. In panels b and c, we have overlaid the continuum contours in white and black, respectively.

In the text

Fig. 2. Line-of-sight velocity maps for v0 data (a) and the discminer model, vmod (b). Panel c: velocity residual map after subtracting vmod from v0. The dust continuum is overlaid in solid contours with equal levels, as in Fig. 1. The innermost region was masked during the fit by one beam size in radius, shown as the gray shaded ellipse. The insets in panels a and c zoom into the innermost region of the disk to highlight the non-Keplerian velocities. Contours are drawn at vsys = (4.62 ± 0.60) km s−1 in steps of 0.1 km s−1 and from −60 to 60 m s−1 in steps of 10 m s−1, respectively. All maps show the synthesized beam for CO (black) and the continuum (white) in the lower-left corner. The beams are masked where the CO peak intensity falls below the 5σ level for panel a and where R > Rout for the rest.

In the text

Fig. 3. Folded velocity residuals (left) and detected clusters of peak velocities (right) in the disk reference frame. The green shaded annulus and wedge in the right plot mark the significant clusters in radius and azimuth, respectively. The position of the localized velocity perturbation inferred from these clusters is marked with a magenta point with error bars. The gray region (one beam size in radius) indicates the masked area, and the black dashed lines the FWHM of the dust ring.

In the text

	Fig. 4. Polar projection of the velocity residual map. The solid black line shows a linear spiral trace. The gray region indicates the masked area and the black dashed lines, the FWHM of the dust ring. The y-axis extends farther than 360° to enhance the visibility of the spiral.
In the text

Fig. A.1. Gallery of selected channel maps. Panels show the 12CO data channel maps (top row) and best-fit model channel maps (middle row), together with intensity residuals in kelvins for each channel (bottom row). In the bottom row the colorbar has been adjusted such that residuals smaller than 1σ are white. The beam size is depicted in the lower-left corner of each channel. For reference, the best-fit systemic velocity was found to be vsys = 4.62 km s−1 and the channel spacing is 100 m s−1.

In the text

Fig. A.2. Radial profile of the surface brightness for different tracers. Profiles are normalized to the peak of the emission for the 231 GHz continuum and the CO peak flux (both for the data and discminer model) and for the SPHERE scattered-light observation. Shaded regions show the standard deviation of each azimuthal average. The lines in the lower-right corner show the major beam size (resolution) for each profile in the corresponding color.

In the text

Fig. A.3. Azimuthal profiles of the surface brightness, normalized to the peak of the emission. Profiles were extracted at an annulus with a width of approximately one corresponding beam size centered at 0.56″ and 0.39″ for the 231 GHz continuum and the CO peak flux, which both show significant azimuthal intensity variations, of 34% and 19%, respectively. Shaded regions show the standard deviation of each radial average.

In the text

Fig. A.4. Additional moment maps of the centroid fitting. Panels show the line width, ΔV (left), the peak optical depth, τ0 (center), and the error of the centroid fitting, δv0 (right). Note that for τ0 < 1 one can assume the line profile to be well presented by a Gaussian, while for τ0 > 5 the line profile has a saturated core, i.e., a very flat top (see Eq. 2). The beam size is depicted in the lower-left corner, and only regions where I0 > 5σ with σ = 1.1 mJy beam−1 are shown.

In the text

Fig. A.5. Polar map of the velocity residuals. Same as Fig. 4, but now overlaid with filamentary structures found by FilFinder. The overplotted red and blue lines are the medial axes of the filamentary structures found by the algorithm. To trace the apparent spiral in the residuals, we restricted the algorithm to search for filaments in the radial locations r = [0.3, 1.25]″. For the filamentary detection, we assumed a smoothing size of one synthesized beam size and a minimum size of 500 pixels for a filament to be considered.

In the text

Fig. A.6. Polar contour map of the centroid residuals. The upper panel shows residuals inside the cavity, the middle panel residuals between the cavity and the outer edge of dust ring, and the lower panel the residuals in the outer disk. The radial spacing between each contour is ∼1.8 au, and the opacity of the lines increases with radius. We would like to emphasize that the bump between PA≈ − (110 − 70)° in the middle panel includes most of the residual points we detect with the peak variance method of discminer, which can be readily seen in the left panel of Fig. B.1.

In the text

Fig. A.7. Comparison of the velocity residual maps for different imaging parameters. The left panel corresponds to the cube used in the main text, while the middle panel corresponds to the same cube without JvM correction. The right panel shows the residual maps for a non-JvM-corrected cube with a different taper (0.15″). In all panels the dust continuum is overlaid in solid contours with equal levels, as in Fig. 1. Best-fit Keplerian models were subtracted from each of the cubes. The detection of the non-Keplerian features is quite robust, irrespective of the imaging procedure.

In the text

Fig. B.1. Location of the folded peak velocity residuals. The detected points are shown in azimuth (left) and radius (middle) and were obtained using the peak variance method. Colors indicate the seven different radial clusters specified, where blue peak residual points are within a detected significant radial cluster. The black crosses are the velocity variances of the clusters plotted at the (R, ϕ) location of each cluster center. The centers of the accepted clusters (those with peak velocity residuals larger than three times the variance in other clusters) in radius and azimuth are marked with vertical black lines in both panels. The right-hand plot shows the normal distribution of the peak residual points in a histogram. Note that outliers of the distribution are related to the localized perturbation. The maximum value of all peak folded centroid residuals is at 0.39″ (57 au), its mean value is 39 m/s, and 1σv = 20 m/s (not to be mistaken with the cluster variances).

In the text

	Fig. C.1. Azimuthally averaged and deprojected azimuthal, radial, and vertical velocity components. The radial width of each annulus is one-fourth of the synthesized beam size. The error bars are given by the standard deviation for each velocity component averaged over the 20 iterations used.
In the text

Fig. C.2. J1604 deprojected velocity components. It is assumed that all velocities are azimuthal (left column), radial (central column), or vertical (right column). For the azimuthal and radial components, wedges along the minor and major axis have been masked as the observations are insensitive to these components (see Eq. C.1). In each panel the synthesized beam is shown in the lower-left corner.

In the text