A kinematically detected planet candidate in a transition disk

,


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, respec- The ALMA 12 CO molecular line emission cube is only available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/670/L1 tively (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(Izquierdo et al. , 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−2130A (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 scatteredlight observations show the presence of two shadows with variable morphology on timescales possibly shorter than a day (Pinilla et al. 2018), and ALMA 12 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 sini) 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 12 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.

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 12 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 12 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 −1 (4.3 K) for CO and 10 mJy beam −1 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 12 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.

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 χ 2 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): 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

Moment maps
The moment maps were computed with bettermoments (Teague & Foreman-Mackey 2018).Since the 12 CO line emission is optically thick, we fitted the following line profile to both the data and model channel maps: where I 0 is the peak intensity of the line and the optical depth, τ (v), varies like a Gaussian; v 0 is the line centroid, τ 0 the peak optical depth, and ∆V the width of the line -where the full width at half maximum (FWHM) = 2 √ ln2∆V -as used in Teague et al. (2022).In Fig. 1b we show I 0 for 12 CO in units of brightness temperature.The corresponding v 0 maps for the data and model are displayed in panels a and b of Fig. 2, respectively.The moment maps for ∆V and τ 0 , as well the error of the line centroid fitting, δv 0 , can be found in

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 0 , of the 12 CO (J = 2−1) line emission.Both dust and gas tracers show a cavity, and the L1, page 2 of 13  [5,15,25,35,45]σ, and the image is plotted with a power-law scaling of γ = 0.6.Panel b:12 CO peak brightness temperature map computed from I 0 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.

Localized velocity residuals
The centroid residual map in Fig. 2c shows a prominent, localized non-Keplerian velocity feature of δv ≈ 60 m s −1 , 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 L1, page 3 of 13 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 2 s −2 and σ R = 0.018 km 2 s −2 (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.

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

A possible warp in the 12 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 0 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.

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 −1 (∼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 L1, page 4 of 13 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).

Origin of v 0 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 −1 (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 13 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 (δ min V =0.06), and its double (δ max V = 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)Mjup (α/10 −3 ) 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 12 CO (J = 3−2) line emission, that partly coincides with the spirallike feature that we detect.Additional observations in optically thin tracers would be very useful to assess its nature.

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 0 residuals broadly align with cold and hot regions in the brightness temperature of 12 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 with H/R the disk aspect ratio, T the 12 CO brightness temperature tracing the gas temperature (as 12 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 12 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.

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 12 CO with non-Keplerian velocities, suggest the presence of an additional, very massive (possibly stellar) companion within the inner ∼0.25 (∼35 au).
L1, page 5 of 13 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 12 CO kinematics (Kurtovic et al. 2022).

Conclusions
In this Letter we present new ALMA observations of the 1.3 mm dust continuum and the 12 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 12 CO cavity.Azimuthal brightness variations in the 12 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).

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 selfcalibration, 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 selfcalibration 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 12 CO (2-1) molecular line observations are imaged with a robust value of 1.0 and a channel width of 0.1 km s −1 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 −1 and 2.9 mJy beam −1 for the non-JvM-corrected cubes with a taper of 0.10 and 0.15, respectively.
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 = ρ g c 2 s = ρ g k B T µ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: (D.2) In the last step, we inserted the gas midplane density, ρ g = Σ g /( √ 2πH) = Σ g /( √ 2π 0.2R), assuming a disk aspect ratio of H/R = 0.2.Assuming that v φ ≈ v kep and Σ g ≈ 1g/cm 2 (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 −6 .Therefore, we neglected the second-order derivative and identified the sound speed, c s : L1, page 13 of 13 Fig. A.4.

Fig. 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: 12 CO peak brightness temperature map computed from I 0 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.

Fig. 2 .
Fig. 2. Line-of-sight velocity maps for v 0 data (a) and the discminer model, v mod (b).Panel c: velocity residual map after subtracting v mod from v 0 .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 v sys = (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 > R out for the rest.

Fig. 3 .
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.
Figure2calso 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 see Fig. A.5).As indicated by the coherent filaments, the strong localized positive

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

Fig. A. 1 .
Fig. A.1.Gallery of selected channel maps.Panels show the 12 CO 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 v sys = 4.62 km s −1 and the channel spacing is 100 m s −1 .

Fig. A. 4 .
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, δv 0 (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 I 0 > 5σ with σ = 1.1 mJy beam −1 are shown.

Fig. A. 5 .
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. 7 .
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.

Fig. C. 2 .
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.
connects the fractional azimuthal temperature variation, ∆T φ /T , to the rotational velocity deviation relative to Keplerian, ∆v φ /v kep .

Table A .
1. Summary of the ALMA Band 6 observations of J1604 presented in this paper.