© ESO 2021

1 Introduction

In order to detect young planets, it is imperative to understand the footprints they leave on protoplanetary discs, their place of formation. Recent observations in concert with theoretical efforts suggest that dust substructures, mostly rings and gaps but also cavities, spirals, and asymmetric features, are possibly ubiquitous in protoplanetary discs (ALMA Partnership 2015; Isella et al. 2016; Pérez et al. 2016; Long et al. 2018; Andrews et al. 2018). In multiple cases, embedded planets may play a key role in shaping some of these substructures observed in infrared and (sub)millimetre wavelengths (see e.g. Benisty et al. 2015; Dipierro et al. 2015; Pinilla et al. 2018; Zhang et al. 2018; Ubeira Gabellini et al. 2019; Facchini et al. 2020). However, planet–disc interactions are far from being the only driving mechanism behind dust signatures in young discs. Magnetic, hydrodynamic, and gravitational instabilities can also lead to dust substructure (Armitage 2011; Andrews 2020), meaning that looking at the dust emission alone is generally insufficient to unambiguously claim the presence of planets. On top of that, thermal and accretion emission from young planets is hard to detect through direct imaging, whose range of action is currently narrowed to massive planets and low-dust-extinction scenarios (Testi et al. 2015; Sanchis et al. 2020). To date, PDS 70 is the only system in which forming planets have been convincingly detected by direct imaging (PDS 70b and PDS 70c, Keppler et al. 2018; Haffert et al. 2019).

Luckily, not only dust but also gas stores valuable information that can help disentangle the physical and chemical processes at work (and often coupled) in planet-forming discs (Bruderer et al. 2012; Henning & Semenov 2013; Dutrey et al. 2014). The presence of deep gas cavities, smaller than those in dust at (sub)mm wavelengths, is a clear diagnostic of the disc interaction with planet and stellar companion(s) (Bruderer et al. 2014; Perez et al. 2015a; van der Marel et al. 2015, 2016). On the other hand, the rich molecular gas disc inventory has been meticulously examined over recent years to reveal density and temperature structure, as well as elemental abundances in a number of objects (Piétu et al. 2007; Rosenfeld et al. 2013; Williams & Best 2014; Miotello et al. 2016; Dutrey et al. 2017; Pinte et al. 2018a; Dullemond et al. 2020; Rosotti et al. 2020; Teague et al. 2020; Facchini et al. 2021). These are all crucial pieces in the vast puzzle of planet formation (Benz et al. 2014; Johansen & Lambrechts 2017; Öberg & Bergin 2021). Moreover, spectral lines from molecules provide a useful window onto gas velocities, and can be ‘mined’ with appropriate modelling to understand the mechanisms driving the disc dynamics (see e.g. Rosenfeld et al. 2013; Price et al. 2018; Teague et al. 2019; Wölfer et al. 2021; Paneque-Carreno et al. 2021), and in consequence, to better constrain the presence of planets.

It is wellknown from previous theoretical works that hydrodynamical and gravitational interactions between forming planets and gas discs produce particular signatures in molecular line observations (Perez et al. 2015b; Pérez et al. 2018). Line Doppler shifts due to localised deviations from Keplerian rotation are expected in circumplanetary material and along spiral wakes launched from the planet. Furthermore, embedded planets carve density gapswhich induce pressure gradients in the gas disc, producing azimuthally extended non-Keplerian velocities at the edges of the gaps (see e.g. Disk Dynamics Collaboration 2020). In this context, and supported by the high angular and spectral resolution offered by the Atacama Large (sub)Millimeter Array (ALMA), many are increasingly turning their attention to spotting the kinematical clues left by planets in protoplanetary discs.

More specifically, Pinte et al. (2018b, 2019) claimed kinematical detections of giant planets embedded in the discs around HD 163296 and HD 97048 by empirical comparison of planet-disc hydrodynamic simulations capable of producing kink-like features similar to those observed in CO intensity channel maps. A handful of new kinks in ¹²CO were reported later by Pinte et al. (2020) but most of them await confirmation because of limitations in signal-to-noise ratio (S/N) and spectral resolution. The kink velocity can be linked to the driving planet mass through simple relationships as recently theorised by Rabago & Zhu (2021) and Bollati et al. (2021), yet both of these latter studies omit radiative transfer effects. On the other hand, Teague et al. (2018a) proposed the presence of two other giant planets in shorter orbits around HD 163296 by looking at azimuthally symmetric gas velocity deviations from Keplerian rotation due to pressure gradients driven by planet-induced gaps. However, this method is limited by the fact that gaps opened by other mechanisms would drive similar kinematical signatures in the disc (see e.g. Rabago & Zhu 2021). Using rotation curves on the disc of HD 100546, Casassus & Pérez (2019) detected a localised ‘Doppler flip’ in ¹²CO reminiscent of the velocity perturbations expected along spiral wakes induced by a planet. Nevertheless, although pivotal, these studies are still source-specific and some lack statistical significance. In particular, the presence of kinks is currently assessed by visual inspection of channel maps, which works fine in some clear cases, but does not allow for an estimate of the significance of the detection, implying that the method can be misleading in less apparent cases.

In this paper, we introduce a statistical framework to overcome these limitations and to robustly detect localised perturbations due to unseen planets using molecular line observations. The technique is also applied to systematically extract observable velocity perturbations driven by different planet masses at a number of azimuths in a synthetic disc. The outline of the work is as follows. Section 2 describes the planet–disc interaction simulations and synthetic observations used throughout the work. Section 3 introduces the DISCMINER package and how we use it to fit Keplerian intensity channel maps on the synthetic observations. Section 4 presents the extraction and analysis of non-Keplerian velocities, as well as the statistical framework designed to infer the location of planets as a function of planet mass and azimuth in the disc. We expand on the advantages of our method and discuss the observable signatures of planets in Sect. 5, and summarise the main results in Sect. 6.

2 Hydrodynamic simulations and radiative transfer

To simulate velocity perturbations triggered by planet–disc interactions, we use the multifluid open source code FARGO3D (Benítez-Llambay & Masset 2016) to solve the hydrodynamic evolution of a gas disc, which responds to thermodynamic variables via the Navier-Stokes and continuity equations, and to the gravitational potential of point-like sources.

Assuming viscous accretion (Lynden-Bell & Pringle 1974; Hartmann et al. 1998), we initialise the simulations with a gas disc surface density profile of the form Σ_gas(R) = ${\Sigma }_{\text{c}}{(R/R_{\text{c}})}^{-\gamma }\exp [-{(R/R_{\text{c}})}^{2-\gamma }]$, with R_c = 100 au as the characteristic radius, Σ_c = 3.0 g cm⁻² the density normalisation at R_c, and γ = 1.0 the viscous power-law exponent. We adopt a uniform kinematic viscosity of α = 10⁻³, and a 1 M_⊙ point-like source at the centre of the mesh. The mesh spans from − π to π in ϕ, 15–700 au in R, and −240 to 240 au in z. The number of cells is 2048, 1250, and 149, respectively, evenly spaced in ϕ and z, and logarithmically spaced in R. The inner and outer radii are assumed to be R_in = 15 au and R_out = 400 au, yielding a gas mass of 0.01 M_⊙.

We ran three independent 2D simulations varying the mass of the embedded planet (0.3, 1.0, and 3.0 M_Jup) at a fixed radial location R_p = 100 au, and let them evolve for 1000 orbits to reach a steady state. Even though the mesh of our simulations is not adaptive, its spatial resolution guarantees that the Hill radius of the planet is resolved (e.g. with ~22 cells for the 1.0 M_Jup planet). Also, the planets are not treated as sink particles, so they are not accreting material. Their gravitational potential is smoothed with 0.6 of the disc scale height, which is the closest 2D approximation of 3D gravity (Müller et al. 2012). On the other hand, we do not find fully developed vortices in any of the simulations at steady state. This is a consequence of the kinematic viscosity being high enough to efficiently quench vortices (see e.g. Fu et al. 2014; Hammer et al. 2017).

To extend the simulations to three dimensions, we assume a cylindrical velocity field. Because of the planet gravity, it is not straightforward to calculate the rotation velocity at different heights in the disc based on the midplane velocities and the central force of the star only, let alone the radial and vertical velocity components. Full 3D simulations are required to self-consistently tackle this. On the other hand, we consider hydrostatic equilibrium to compute the gas volume density along the vertical direction, ${\rho }_{\text{gas}}(R,z)=\left({\Sigma }_{\text{gas}}(R)/\sqrt{2\pi }H\right)$ $\exp [-0.5{(z/H)}^2]$, where $H(R)=H_0{(R/100\text{au})}^{\psi }$ is the scale height of the disc, with a normalisation H₀ = 6.5 au, and a flaring index ψ = 1.25. Figure 1 shows edge-on and face-on slices of the resulting gas number density and intrinsic deviations from Keplerian rotation for the 1.0 M_Jup simulation.

Based on the simulated gas densities and velocities, we perform radiative transfer calculations to obtain synthetic emission maps of ¹²CO J = 2−1 assuming a parametrised CO abundance (see below). Velocity binning, line-of-sight projection and optical depth effects are therefore considered in the analysis. These are unavoidable factors that must be taken into account if one is to study any observable, especially when it is expected to be faint and confined to small scales as is generally the case for planet-driven perturbations (Pérez et al. 2018). In real data, additional sources of uncertainty such as noise in the signal, beam smearing, and non-linear artefacts from image reconstruction algorithms are also important (see e.g. Disk Dynamics Collaboration 2020) and will be explored in future releases of the Disc Miner series, which are more focused on observations.

We use the RT modules of the SF3DMODELS package (Izquierdo et al. 2018) to bridge our FARGO3D simulation snapshots with the POLARIS radiative transfer code (Reissl et al. 2016). We first run POLARIS to compute the three-dimensional dust temperature, which is done by propagating photon packages semi-randomly from a T Tauri star at the centre of the disc, with a luminosity of 1.9 L_⊙ and a photosphereat an effective temperature of 4000 K. We assume a standard Mathis radiation field at a galactocentric distance of 10 kpc as anexternal source of radiation (Mathis et al. 1983). The disc is located 100 pc away from the observer and inclined at − 45° with respect to the plane of the sky, with the north half being the side closest to the observer. For simplicity, the position angle of the disc was fixed at 0°. The dust model consists only of silicate grains with an intrinsic density of 3.5 g cm⁻³ (with optical properties from Weingartner & Draine 2001). We use a standard ISM gas-to-dust mass ratio of 100, and a grain size distribution n(a) = a^−3.5, between a_min = 5 nm and a_max = 1 mm. We consider a direct conversion of temperatures T_gas = T_dust, which is a good assumption at scale-heights z∕R < 0.3 where our analysis takes place (see e.g. Kamp & Dullemond 2004; Jonkheid et al. 2004; Miotello et al. 2014). Second, we use SF3DMODELS again to read the output temperatures and add simplified chemical processes such as CO freeze-out on dust grains, where T_dust < 20 K, and photodissociation, at scale heights z∕R > 0.3. We adopt a ¹²CO abundance of 5 × 10⁻⁵ in the gas phase, and 5 × 10⁻¹¹ in freeze-outand photodissociated regions. Lastly, SF3DMODELS provides POLARIS with the gas velocity, temperature, and the processed CO abundance distribution, which this time are used to compute the ray-traced intensity channel maps of the simulations.

In order to consider projection effects on the planet-driven perturbations, we ran the ray-tracing 13 times per planet, each time rotating the input gas properties along the vertical axis of the disc from –90° to 90°, in steps of 15°, so that the planets can be observed at 13 different azimuths in the posterior analysis (see Fig. 1). As such, we end up with 39 position-position-velocity (ppv) cubes, each with 101 channel maps ranging from −5 to 5 km s⁻¹, taking a channel width of 0.1 km s⁻¹. We adopt a pixel size of 3.1 au which translates to 31 mas at 100 pc. A selection of these channel maps is presented in Fig. 2, where we also compare some of the kinematical features triggered by the three planets in our sample. There are many kinks evident to the eye in the channel maps, especially around massive planets. However, we note that the gap alone can also produce kink-like features, suggesting that empirical methods that rely on visual inspection are prone to false positive detection of planets. Conversely, as explained in Sect. 4, our statistical method does not get confused by these apparent features.

Fig. 1 Left panel: edge-on gas number density of one of our planet–disc simulation snapshots (1.0 MJup, ϕp = 0°). The white contour encloses the CO freeze-out region (T < 20 K). The black lines correspond to z∕R = 0.1, 0.2, 0.3, 0.4 scale heights, with the solid line being the threshold adopted for photodissociation. Middle panel: face-on view of the gas number density in the midplane of the disc. Right panel: azimuthal (main panel) and radial (zoom-in) deviations from Keplerian velocity. The solid and dashed contours in the zoomed-in panel are ±60 per cent peak azimuthal and peak radial perturbations, respectively, and illustrate that both components do not necessarily overlap and in turn contribute independently to the observed peak velocity perturbations. Also shown is the planet Hill sphere, with a radius of 6.8 au. The green circle indicates the current position of the planet, and the grey circles are additional planet azimuths explored in this work.

3 Fitting channel maps

3.1 The Discminer package

In this section we present the basic functionality of the DISCMINER code, carefully designed to capture structural and kinematical features from circumstellar discs using molecular line emission. The open source version of the DISCMINER will be released with the second paper of the Disc Miner series.

As in previous fitting methods that dig into the kinematics of discs (e.g. Teague et al. 2018b; Casassus & Pérez 2019), our package is based on parametric descriptions of the physical and geometrical properties that make up a simplified disc of gas¹. However, the DISCMINER provides for the first time a simultaneous representation of line profiles and velocities by fitting individual intensity channel maps rather than projected velocity maps (see Fig. 3). To do so, the line intensity is modelled considering the following attributes,

• (i)

  We assume that the disc emission comes from two thin (upper and lower) surfaces. Their vertical location is described by two parametric prescriptions of the form z = f(R).

• (ii)

  A Keplerian velocity field υ_k = f(R, z;M_⋆, υ_sys), which is used to determine the shift in velocity space of the emission from any given pixel.

• (iii)

  A kernel to shape the line intensity profile as a function of the disc coordinates. The kernel combines peak intensity (I_p), line-width (L_w) and line-slope (L_s) information, which are also parametric prescriptions of the form A = f(R, z). See further details in Sect. 3.2.

• (iv)

  Inclination and position angle of the disc. These are used for geometrical projection of intensities and line-of-sight velocities onthe plane of the sky.

The model is then coupled with a Markov chain Monte Carlo (MCMC) random sampler, EMCEE (Foreman-Mackey et al. 2013), which efficiently walks over a vast range of parameters to determine the subset of them that best reproduces the projected line intensity of the input disc, often encoded in a three-dimensional ppv cube. We configure the MCMC sampler to maximise a χ² log-likelihood defined as, $${\chi }^2=-0.5{\sum }_j^{n_{\text{ch}}}{\sum }_i^{n_{\text{pix}}}{w}_i^{-2}{\left[I_{\text{m}}(r_i,{\upsilon }_j)-I_{\text{d}}(r_i,{\upsilon }_j)\right]}^2,$$(1)

where the index i runs over the pixel location in (x, y)_sky coordinates, and the index j runs over the input velocity channels. The χ² kernel is the difference between the model intensity, I_m, and the input intensity, I_d, whose uncertainty is encompassed in the weighting factor w which is computed using residual intensity from line-free velocity channels.

Fig. 2 Selected 12CO J = 2−1 synthetic channel maps for the three simulation snapshots (0.3, 1.0 and 3.0 MJup, from top to bottom), with the planets at Rp = 100 au radius, and ϕp = 45° azimuth, marked as green circles. The disc is inclined at − 45° with respect to the plane of the sky. The thin grey line is the projected circular orbit of the planet. The solid lines are centroid velocity contours extracted from the simulations at the velocity channel indicated on the top header; the dashed lines show the same but for the Keplerian best-fit model. Small arrows indicate kink-like features identified by visual inspection. In the bottom row a zoom-in around planets is shown for better comparison of centroid velocities and their deviation from Keplerian rotation as a function of planet mass.

Fig. 3 Summary of the main attributes making up the line emission of a disc in the DISCMINER. Left panel: projected intensity of the disc for a channel centred on υch = − 1.0 km s−1. Right panel:line intensity profile extracted from the marker on the left. The grey annotations indicate the role of each attribute listed in Table 1.

3.2 Line intensity profile

In this work, we assume a generalised bell function for the line profile kernel, $$I_{\text{m}}(R,z;{\upsilon }_{\text{ch}})=I_{\text{p}}{\left(1+{\left|\frac{{\upsilon }_{\text{ch}}-{\upsilon }_{{\text{k}}^{\text{l}.\text{o}.\text{s}}}}{L_{\text{w}}}\right|}^{2L_s}\right)}^{-1},$$(2)

where I_p is the peak intensity, υ_ch is the channel velocity and ${\upsilon }_{{\text{k}}^{\text{l}.\text{o}.\text{s}}}$ is the Keplerian line-of-sight velocity. The choice of this kernel is motivated by the fact that at the cost of only one additional parameter (L_s), it performs better than a Gaussian function at reproducing optically thick lines (such as those from ¹²CO transitions) which are flat at the top and decay rapidly towards the wings. The line width (L_w) of the generalised bell function is the half width of the profile at half power. The additional (dimensionless) parameter, the line slope (L_s), controls how steep the signal drops at the wings and in turn also determines the spectral extent of the plateau at the top of the profile. It is easy to deduce from the first derivative of the bell function that the slope of a line tangent to either points of the (normalised) profile, at half power, is the ratio 2L_s ∕L_w.

For simplicity, we parametrise the peak intensity (I_p) and line width (L_w) as power laws of the disc cylindrical coordinates (R, z), which reproduce the overall line profiles of the synthetic observations reasonably well (see Sect. 4). The line slope (L_s) is allowed to vary but it is assumed uniform everywhere. We note that the model attributes perform well at describing the observed emission on the upper and lower surfaces of the input disc only, and hence any extrapolation to other scale heights should be done with caution.

3.3 Best-fit model of the simulation

To reveal deviations from Keplerian rotation, we use the DISCMINER to fit Keplerian channel maps on the simulated synthetic observations. In particular, we obtain three best-fit models, one for each planet mass (0.3, 1.0 and 3.0 M_Jup) at 0° azimuth. To achieve this, we first use the prototyping tool² of the DISCMINER to find a reasonable set of seeding parameters for EMCEE. We then let 256 walkers evolve over 1000 steps in a first trial, and allow them to fully converge to the final parameters over another 1500 steps. As our simulations are noise-free, the weighting factor w in Eq. (1) is equal to one everywhere. Each run takes about six hours on a 48-core machine with a clock rate of 2.3 GHz per core.

Table 1 presents a summary of the functional forms considered for the model attributes and the best-fit parameters obtained for the three snapshots. Figure 4 shows the best-fit attributes computed for the 0.3 M_Jup snapshot, deprojected on the upper and lower emitting surfaces of the disc, and the vertical location of both the surfaces (see parameter walkers in Appendix B). We notice slight variations in the best-fit parameters of the other two (1.0 and 3.0 M_Jup) snapshots due to the combined effect of the planet mass and the depth of the gap it carves. More specifically, the stellar mass (M_⋆) retrieved by the model increases with planet mass (M_⋆ = 1.013, 1.025M_⊙, respectively), while the emission surfaces are shallower because of the increasingly deeper gaps (e.g. z₀ = 15.59, 14.35 au for the 1.0 and 3.0 M_Jup upper surfaces). Nevertheless, these variations have negligible impact on the detection and quantification of planet perturbations.

Taken together, the parameters retrieved by the model converge notably well to the features we knew beforehand from the input simulations. The mass of the star (1 M_⊙), the disc inclination (− 45°), and the height of the lower surface tracing the back side of the CO freeze-out region are all closely reproduced even when a planet is present. This means that the DISCMINER is also well suited to future studies of the three-dimensional structure of discs using molecular line emission.

Fig. 4 Best-fit attributes obtained by the DISCMINER for the 0.3 MJup snapshot. Left panel: freeze-out and gas phase regions for 12CO. The coloured circles indicate the height of the model emitting surfaces. Right panels: model peak brightness temperature and line width as a function of radius on both emitting surfaces.

4 Finding the kinematical footprints of planets

4.1 Intrinsic deviations from Keplerian rotation

Before analysing observables, it is worth looking at the pristine simulation velocities to understand the posterior effect of radiative transfer and disc structure on the retrieved planet-driven perturbations. In Fig. 5, we present peak line-of-sight velocity deviations from Keplerian rotation (in magnitude) for all three planet masses, varying planet azimuths from − 90° to 90°. We use the same disc inclination (i = −45°) adopted for the synthetic observations, and assume a line-of-sight parallel to the vertical axis. There are two overall features that stand out here. First, it is clear that the projected peak deviations are different for positive and negative planet azimuths. This is because the outer and inner spiral wakes launched from the planet, are geometrically different, and aside from direction, perturbations along these spirals generally differ in magnitude as well. Hence one or the other dominates depending on the planet mass (inner[outer] spiral perturbations are higher for the 0.3[1.0, 3.0] M_Jup planet), and on the projected direction of the fluctuations at each planet azimuth. Second, in all cases there is a sharp turnover azimuth where the radial component of the perturbation starts dominating over the azimuthal one. This is explained by the fact that, unlike azimuthal perturbations, the highest radial perturbations do not necessarily occur along spiral wakes but on circumplanetary material (see Fig. 1, right panel), implying that the radial and azimuthal perturbations generally contribute independently to the projected perturbation.

In the following sections, we perform a similar analysis on synthetic observations of the simulations and find that radiative transfer anddisc structure can strongly influence the observed velocity fluctuations as compared to the actual gas velocities.

4.2 Residual maps

In order to extract observables and quantify any line profile differences between the simulation and the smooth Keplerian (best-fit) model, we fit a Gaussian profile to each pixel on both the simulation and model intensity cubes for each planet mass and azimuth (see Appendix A). The mean, amplitude, and standard deviation of each Gaussian profile are assumed to be the line centroid, peak intensity, and line width of the corresponding pixel. We then subtract the best-fit model line profile properties from those of the simulations and produce residual maps as illustrated in Fig. 6. There, we conveniently present a full scan of residuals from one of our snapshots as a function of azimuth along constant radii contours whose colours represent their closeness to the actual radial distance of the planet (R_p = 100 au).

A number of features of these residual maps are worthy of discussion. First, the three types of residuals are all mostly uniform along the green and red contours which correspond to the outermost parts of the disc, meaning that the simulation and the model line profiles are close to identical in regions away from the planet. Second, residuals along the blue contours are rather high and are mainly associated with the gap carved by the planet (see Appendix A), whose contribution is not considered by the smooth Keplerian model. However, in all cases these residuals are symmetric along the projected minor axis (hereafter the ‘symmetry axis’) of the disc, ϕ = ±90°. Unsurprisingly, this symmetry is notably disrupted by the contribution of the embedded planet, whose strongest effect is localised both in radius and azimuth. In Sect. 4.3, we exploit these (a)symmetries to isolate kinematical perturbations driven by planets as a function of their mass and azimuth in the disc, and distinguish them from other perturbations.

Fig. 5 Peak line-of-sight deviations from Keplerian rotation versus positive and negative planet azimuths, for all three planet masses, extracted from the native simulations (without radiative transfer effects). The diamonds highlight turnover points where the radial velocity perturbations become important and eventually dominant over the projected azimuthal perturbations.

4.3 Line centroid folding

Taking advantage of the fact that the strongest kinematical perturbations driven by planets are spatially and spectrally localised, we propose a line centroid folding method to remove any symmetric contribution to the velocity field arising from the natural rotation of the disc and from the large-scale contribution of the gap. The method consists of subtracting line centroid velocities from one half of the disc from those of the other half, exactly as if the disc was folded along its symmetry axis, ϕ = ±90°.

We illustrate the outcome of folding line centroids in Fig. 7, for all three planet masses, namely 0.3, 1.0 and 3.0 M_Jup, and three azimuthal locations, 0°, 45° and 90°. Unlike in the raw residual maps, the highest velocity residuals are this time closely related to the embedded planet thanks to the localised nature of the perturbation, whereas most of the contribution from the gap is cancelled out. Of particular interest, the magnitude of the residuals, as well as their spatial extent, appear to be tightly linked to the mass and azimuthal location of the planet. Quantitatively, typical peak centroid residuals range within 40−70 m s⁻¹, 70−170 m s⁻¹, and 130−450 m s⁻¹ for the 0.3, 1.0, and 3.0 M_Jup planets, respectively. At ϕ_p = 90° azimuth, all planets trigger the lowest velocity residuals, and from ϕ_p = 30° to 60°, the highest. The angular dependence of residuals can be readily understood by noting that at ϕ_p = 90° and neighbouring angles, most of the azimuthal component of the planet perturbation is cancelled out because it is orthogonal to the line-of-sight. The equivalent occurs around ϕ_p = 0° azimuth for the radial component of the perturbation. Nevertheless, interpreting the high-velocity fluctuations observed at intermediate angles is less straightforward. The observed fluctuations at these angles do not only depend on the intrinsic magnitude of the perturbation, but also on the scale height at which the perturbation is measured and hence on the structure of the disc itself. This effect is further discussed in Sect. 5.3.

Now, we discuss whether or not it is possible to determine not only the presence, but also the actual location of the planets from the information gathered immediately above. From Fig. 7, this appears to be the case, especially for those planets at azimuths smaller than90°, where the perturbation is obvious to the eye. However, line centroid residuals from less massive planets, such as the 0.3 M_Jup included in our analysis, may be particularly challenging to detect just by visual inspection because the localised signature, intrinsic to the planet, is small and can easily get confused with other asymmetric, though more extended perturbations. Also, planets of any mass near the symmetry axis of the disc, ϕ = 90°, are equally challenging because their signature is no longer seen as localised to the naked eye. In the following section, we provide two statistical methods that conveniently examine line centroid residuals and do not rely on visual inspection to robustly infer the location of embedded planets in discs.

Fig. 6 Line centroid (top), peak intensity (middle), and line width (bottom) residuals for the (1.0 MJup, ϕp = 45°) snapshot. The azimuthal scans on the left run along constant radii contours in disc coordinates; their colours represent their closeness to the radial location of the planet, with blue being the closest. The solid black contour runs along the projected distance of the planet (Rp = 100 au) and the dashed black line shows the azimuth of the planet.

Fig. 7 Folded (magnitude of) line centroid residuals as a function of azimuth, varying the analysis radius (in colours), for different planet masses (0.3, 1.0 and 3.0 MJup, from left to right) and azimuths (0°, 45°, and 90° from top to bottom). The solid black contour runs along the projected distance of the planet (Rp = 100 au), and the dashed black line corresponds to the azimuth of the planet.

4.4 Statistical methods to detect planets

In order to detect an embedded planet through kinematical analysis, one needs to be able to spatially isolate the velocity fluctuations that the planet induces on the gas disc, but also to provide a robust measurement of their magnitude. Here we present two statistical methods to detect and quantify planet-driven kinematical perturbations using synthetic observations of simulations with different planet masses and azimuths (see Sect. 2). Both methods make use of the folded centroid residuals presented in Sect. 4.3.

From inspection of Fig. 7, a first natural approach would be to determine the location of the global maximum centroid residual, which is apparently well correlated with the azimuthal and radial location of the planet. Following this idea, we extract the peak velocity residual for each radius in the disc as a function of the azimuth where it occurs. From the resulting distribution of velocity residuals, we compute a 3σ threshold below which residuals are discarded from the analysis. The inferred location of the planet is assumed to be the median azimuth and the median radius of the leftover peak residuals, as illustrated in Fig. 8. This procedure is what we call the Global Peak detection method. In Appendix B.2, we show the dependence of peak fluctuations as functions of planet mass by fitting δυ = aM^b powerlaws, for each planet azimuth, with a and b as free parameters. We note that there are substantial differences between the observed fluctuations at negative and positive planet azimuths. As explained later in Sect. 5.3, the reason is that the disc vertical structure and the gap – combined with projection effects – are important, and their contribution to the observed velocity fluctuations depends on the location of the planet.

Typically, after comparing the inferred and the actual location of planets, the Global Peak method is accurate within ± 3° in azimuth and ± 8 au in radius, but its significance, which in all cases is between 3 and 7σ, might be a weak spot in noisy scenarios. This leads us to introduce a rather different method, based on the assumption that the strongest velocity fluctuations driven by planets should be simultaneously localised and coherent, which in other words implies that the velocity field, as well as its observables, should vary smoothly as one approaches the planet. If this is true, in addition to large velocity residuals there should also be a higher density of peak velocity residuals around the planet compared to undisturbed regions of the disc. This is pictured in Fig. 9, where we find clusters of peak residuals, independently for each spatial coordinate (R, ϕ), using a K-means clustering algorithm (MacQueen 1967; Lloyd 1982). The K-means algorithm subdivides the input residuals into a predefined number of clusters in such a way that the centre of each cluster is the closest centre to all the residuals in the cluster. In other words, the input data are iteratively partitioned into Voronoi cells until convergence is reached, which in this case means until the sum of squared distances from the data to the centre of their clusters is minimised. We note that the minimisation distance used by the iterative procedure to identify clusters in the top panel of Fig. 9 is defined in 2D (azimuth, velocity residual) space, but in practice the azimuthal distance dominates over the distance in velocity residuals so that the clustering is similar to a binning in azimuth (the same applies for the radial clustering). However, unlike simple manual binning, the bin boundaries are in this case irregular, and the bin centres tend to be near the densest accumulation of points, which is ideal for localising coherent velocities. We adopt ten clusters such that the azimuthal and radial extent of the localised perturbation from the planets (δϕ_p ≈ 20° − 50°, δR ≈ 50−100 au, at a radius of 100 au) is always within one to three clusters. Next, the cluster with the highest velocity variance is attributed to the planet-driven perturbation as long as it meets the requirement of being above the standard 3σ threshold, which refers to the variances of the background clusters. In this case, the retrieved azimuthal location of the planet is the azimuth of the centre of the cluster with the peak variance. The same recipe is followed to infer the radial location of the planet³. This method, the Variance Peak, is equally accurate to the Global Peak method when it comes to determining the location of planets. However, it strikingly boosts the significance of the detection thanks to the coherent nature of velocities around planets.

The absolute value of peak residuals in the Global Peak method, and the peak variances in the Variance Peak method are presented inFig. 10. Unsurprisingly, regardless of the planet azimuth, peak fluctuations and peak variances increase steadily with the mass of the planet. Also, because of projection effects, the highest peaks are reached at intermediate angles (30°, 45°, 60°). Furthermore, all of the detected radial locations are approximately within R = 100 ± 10 au except for the ϕ_p = 90° snapshots where there is a higher dispersion of peak residuals, making it more challenging for the technique to find the radial location of the planet. Interestingly, we note that for the same ϕ = 90° snapshots, the azimuthal location of all planets is reasonably well determined, as opposed to the first impression left by the bottom row of Fig. 7, where velocity residuals were anything but localised to the naked eye.

However, of greater interest is the statistical significance of both methods. Even though both approaches work well at inferring the location of planets for all masses and azimuths, the Variance Peak is far more robust than the Global Peak method. As illustrated in Fig. 11, the significance of the Variance Peak detections are in almost all cases larger (or even much larger) than 10σ, whereas the Global Peak detections are always between 3 and 7σ. Again, this is thanks to the coherence of the velocity field around planets leading to a significant accumulation of peak residuals to which the Variance Peak method is fairly sensitive. The Global Peak method, on the other hand, is not sensitive to the bulk of velocity residuals around planets, but is limited to only the highest of them.

We note that massive planets such as the 3.0 M_Jup planet in our sample do not always yield the highest detection significance. The reason behind this is mainly that the velocity perturbation from massive planets is so spatially extended that some of the velocity residuals associated with the planet are actually contributing to the background residuals, which in turn reduces the significance of the detected global and variance peaks. We softened this effect in the Variance Peak method by allowing more than one cluster to be considered as part of the planet perturbation so that the background is better constrained prior to calculation of the significance of the detection. In any case, the fact that the least massive planets in our sample (0.3, 1.0 M_Jup) are in almost all cases robustly detected is encouraging, given that such low masses are almost undetectable by empirical methods.

Fig. 8 Left panel: peak residuals and their azimuthal location for the (1.0 MJup, ϕp = 45°) snapshot. Right panel: peak residuals rearranged as a function of their radial distance. The panel on the right shows the peak residual distribution (green line) and 1, 2, and 3σ significance thresholds (dotted lines). The solid black lines are the azimuthal and radial location of the global peak, while the dashed black lines show the actual location of the planet.

5 Discussion

5.1 Comparison with other methods

The fitting methods encompassed in the DISCMINER have proven to be effective at describing the large-scale structure of discs but also to be accurate when detecting localised signatures on intensity and velocity driven by embedded planets. In particular, our fitting technique has a number of advantages over previous modelling efforts. First, it can describe intensity, line width, and velocity field simultaneously. Second, it is able to reproduce the height and line properties of both the upper and lower emitting surfaces of the disc independently, which is essential to properly assessing the flux and kinematics of discs. Third, it consistently gains velocity accuracy by fitting multiple channel maps at once, allowing for subspectral measurements of small-scale perturbations even on cubes with standard spectral resolutions (e.g. Δ_chan = 100 m s⁻¹).

Modelling all these ingredients together is key to understanding the actual contribution from embedded planets to the velocity perturbations, whose observables are secondary products that depend on the underlying intensity of the disc. To illustrate this, in Appendix A we manually removed the planet from one of our simulations and imposed the gas velocities to be either fully Keplerian or non-Keplerian but azimuthally symmetric. We then noticed that the wavy behaviour of the velocity residuals displayed in Fig. 6 is mostly driven by the gap, or more specifically, by the fact that the differences between simulation and model line intensities on the upper and lower emitting surfaces of the disc are more prominent at the location of the gap. This discrepancy causes the resulting line centroids to be red- or blueshifted in relation to one other depending on their projected location only. One way to work around this is by comparing velocities at peak intensities (namely on the upper surface of the disc only), but at the cost of a substantial loss of velocity accuracy as this method is closely dependent on the channel width of the data. Additionally, as explained in Sect. 5.3, the observed velocity perturbations may be stronger and more clearly projected on the lower surface of the disc depending on the planet location, in which case they would go unnoticed by peak intensity methods. Our line centroid folding procedure gets rid of the unwanted contribution of the gap while keeping both (sub)spectral resolution and information from the lower surface of the disc. Equivalently, as illustrated in Appendix A, this procedure allows us to distinguish between azimuthally localised planet-driven fluctuations and axisymmetric non-planet-induced velocity perturbations, with the second scenario effectively leading to non-detections.

On the other hand, previous techniques that rely on visual inspection of kink-like features (e.g. Pinte et al. 2018b, 2020) are unable to find low-mass planets (<1 M_Jup) or planets of any mass near the projected minor axis of the disc, because kinks are camouflaged with the background velocities in both of those contexts. Our statistical method can instead capture localised, coherent velocities, and therefore detect planets even when kinks are not visually manifest. This also implies that even though random noise or large-scale fluctuations may be comparable in magnitude to the localised planet-driven velocity fluctuations, the latter are more densely assembled and therefore more easily detected by our technique.

Fig. 9 Top row, left panel: KDE contours derived from K-means clusters of peak velocity residuals for the (1.0 MJup, ϕp = 45°) snapshot. The contours enclose 33 and 95 per cent of the peak residuals. The white crosses are the cluster centres retrieved by the K-means algorithm. Right panel: spectral variance of peak residuals for each cluster on the left. In the panel attached on the right we show 1, 2, and 3σ significance thresholds. The yellow line highlights the azimuthal location of the peak variance and the dashed line the actual location of the planet. The green cross is the peak variance significant enough to be attributed to the planet perturbation. Bottom row: 2D visualisation of the detection using the Variance Peak method for three planets (0.3, 1.0 and 3.0 MJup), all at ϕp = 45° azimuth and Rp = 100 au radius. The green circles and crosses are the actual and inferred location of the planets, respectively. The highlighted sectors are the azimuthal and radial clusters with the highest spectral variance extracted by the method. The boundaries of the sectors mark the maximum spatial coverage of the clusters. The circles show the location of the observed peak velocity residuals, whose magnitude is indicated by their size. The red circles are residuals within the azimuthal or radial peak variance clusters. The solid (δυϕ) and dashed (υR) contours trace spiral wakes and radial planet perturbations, and correspond to 60 per cent peak velocity fluctuations extracted from the simulation velocities.

Fig. 10 Peak variance and peak velocity fluctuations extracted with the Variance and Global peak methods as a function of the detected planet azimuth for all three planet masses. The detected radial distance is shown in the bottom panels. Empty circles indicate the actual locations of planets.

Fig. 11 Significance of planet detections using the Variance and Global peak methods as a function of planet azimuth for all three planet masses. The red line highlights the 3σ threshold considered to validate a detection. Empty circles indicate the actual location of planets.

5.2 Robustness of the spectral resolution of the method

We assessed the robustness of the spectral resolution of our method by performing synthetic observations with the same setup as in Sect. 2 but using half the original channel width this time (i.e. Δ_chan = 50 m s⁻¹). The detected planets and the significance of the measurements remained almost the same except for a single snapshot (0.3 M_Jup, ϕ_p = 75°), where the detection was indeed closer to the actual location of the planet. This suggests that our analysis is already very close to measuring velocity fluctuations without suffering from channelisation effects.

5.3 Magnitude of the observed planet-driven perturbations

Our findings suggest that the observable velocity fluctuations driven by planets depend on the intrinsic magnitude of the perturbations, given by Fig. 5, but also on the disc vertical structure and the depth of the gap carved by the planet.

Figure 12 demonstrates the impact of disc structure on the observed velocity fluctuations for a 3.0 M_Jup planet at ϕ_p = +30° and ϕ_p = −30° azimuths. Despite the fact that the intrinsic peak fluctuations are almost the same for both angles, the observed peak fluctuations differ significantly. To understand this, we first highlight the fact that the emission height of the perturbation determines the background velocity (at the other side of the disc) against which it will be compared to compute the magnitude of the perturbation. If the test background velocity falls within the gap, it will trace a different Keplerian velocity from the one it would trace at the same emitting surface of the perturbation. More specifically, perturbations projected outside the orbit of the planet will be overestimated, whereas those projected inside the orbit will be underestimated, by a factor that depends on the depth of the gap. For a similar reason, the observed peak fluctuations are generally dominated by one of the emitting surfaces of the disc, which at the same time determines how much the line profile centroid can be shifted. In particular, the brighter, upper surface of the disc shifts line centroids more than the fainter, lower surface. This effect is especially prominent for massive planets, which carve deeper gaps, and for intermediate planet azimuths, where the projected perturbation is more extended and the observed line profiles are clearly shaped by both emitting surfaces. No less important is the role of the disc inclination, which determines how much of thelower surface can be observed and how far apart the emission from both surfaces is in the velocity space.

With this example as background information, it is worth discussing some of the trends obtained for the observed peak velocity fluctuations. As illustrated in Fig. 10, planets at intermediate angles are better detected and yield the highest velocity residuals, δυ = 70, 170, 450 m s⁻¹ for 0.3, 1.0, and 3.0 M_Jup, respectively.The retrieved velocity fluctuations are enhanced in those cases by the combined contribution of the gap and projection effects on the upper and lower emitting surfaces of the disc, as explained above. Overall, the retrievable velocity fluctuations from planets seem periodic as a function of azimuth, and the amplitude of such a pattern correlates well with the mass of the planet. It is also interesting to analyse the limiting planet azimuths ϕ_p = 0°, 90° because they provide constraints on the orthogonal components of the perturbation by cancelling out one of them by line-of-sight projection. For the ϕ_p = 0° case, we obtain peak residuals δυ_ϕ = 40, 120, 220 m s⁻¹ for each planet mass, where the subscript ϕ indicates that the observed perturbation is nearly fully azimuthal. Similarly, for the ϕ_p = 90° snapshot we find peak residuals δυ_r = 50, 75, 130 m s⁻¹. These limiting cases better match the trend of the intrinsic peak deviations from Keplerian presented in Fig. 5, where the azimuthal perturbations are often stronger than the radial ones. However, we note that the observed peak perturbations tend to be lower than the intrinsic ones. This is because most of the background reference gas at the opposite side of the perturbation stands on the non-Keplerian edges of the gap, which in turn softens the magnitude of the observed perturbation around the planet.

Fig. 12 Illustrating how line centroids are shifted towards faster velocities around a 3.0 MJup planet, at two azimuths ϕp = ±30°. The crosses on the left indicate the location of the line profiles with the same colours on the right. The blue cross is the location of the observed peak perturbation, and the pink cross is the mirror pixel at the other side of the disc. We note that depending on the planet azimuth, the peak perturbation is projected either on the upper or on the lower emitting surface, which contribute differently to the observed line profiles. This contrast leads to variations in the observed line centroids (vertical black lines on the right) and the retrieved velocity fluctuations (in red).

5.4 Caveats

We warn the reader that we do not claim to present a full characterisation of planet-driven perturbations, nor do we include an exhaustive list of all variables that real discs encompass. Instead, the goal of the paper is to provide a new methodology for detecting embedded planets through careful inspection of line emission profiles under certain disc conditions. Nevertheless, the technique can potentially be extrapolated to other disc scenarios. As such, to keep degeneracy at its lowest level, we did not include a full three-dimensional treatment of the gas velocities. We assumed simple cylindrical rotation implying that the central force from the star, per unit mass, is simply f_R = GM_⋆∕R². This means that there is no differential rotation along the vertical coordinate of the disc, which leads to errors of between ~ 40−60 m s⁻¹ in the rotation velocities as measured on the model emitting surfaces. However, because of the axisymmetric nature of this effect, it does not have any impact on the detection of planet-driven perturbations with our method. On the other hand, the vertical gravitational pull from the planet and hydrodynamic meridional flows driven by the carved gap are not considered either. Also, some observational biases such as noise and beam smearing are excluded from the analysis. Future releases of the Disc Miner series more focused on observations will incorporate these effects.

6 Conclusions

We introduce a novel statistical technique to detect kinematical perturbations driven by embedded planets in discs. The method is sensitive to localised deviations from Keplerian rotation by examining line centroid differences using intensity channel maps, which allows us to locate and quantify velocity fluctuations around planets with high accuracy, all while preserving line width and intensity information. Our approach is powered by the DISCMINER package, which aims to model channel maps by simultaneously fitting the intensity, rotation velocity, and height of the upper and lower emitting surfaces of the disc. The package was originally developed for kinematical analyses, but it is also well suited to studying the three-dimensional structure of discs.

We tested this new method on synthetic observations of the ¹²CO J = 2−1 line from simulations of planet–disc interactions to explore variations in the disc kinematics as a function of the planet mass (0.3, 1.0 and 3.0 M_Jup) and planet azimuth (from − 90° to 90° in steps of 15°) for a disc inclination of − 45°. As expected, the observed velocity fluctuations increase with planet mass. In all cases, the highest deviations from Keplerian rotation are found at intermediate azimuths, ϕ_p = ±(30°, 45°, 60°), which are strongly enhanced by the influence of the gap and the vertical structure of the disc combined with projection effects. The lowest velocity fluctuations are obtained for planets on or near the projected main axes of the disc (i.e. along ϕ = 0° or ϕ = ±90°). The method can detect all planets at all azimuths, despite the fact that some of them do not exhibit clear kinks in the channel maps, such as the 0.3 M_Jup planet, or planets of any mass near the projected minor axis of the disc, ϕ = ±90°. Likewise, our approach does not get confused by apparent kinks triggered by gaps in the gas disc or, more generally, by any axisymmetric velocity perturbation field, regardless of its origin.

Our technique takes advantage of the coherence of velocity fluctuations around planets to boost the detection of planet-driven kinematical perturbations in gas discs. We find this to be a substantial improvement to the previous methods for determining whether or not localised velocity fluctuations are unambiguously manufactured by embedded planets.

Acknowledgements

The authors would like to thank the anonymous referee for their constructive comments that helped improve the manuscript. This work was partly supported by the Italian Ministero dell Istruzione, Università e Ricerca through the grant Progetti Premiali 2012 - iALMA (CUP C52I13000140001), by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) - Ref no. FOR 2634/1 TE 1024/1-1, and by the DFG cluster of excellence Origins (www.origins-cluster.de). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 823823 (DUSTBUSTERS) and from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). S.F. acknowledges an ESO fellowship. G.R. acknowledges support from the Netherlands Organisation for Scientific Research (NWO, program number 016.Veni.192.233) and from an STFC Ernest Rutherford Fellowship (grant number ST/T003855/1). This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure.

Appendix A Impact of the gap and emitting surfaces on the observed gas velocities

Our analysis indicates that gaps and emitting surfaces must be taken into account to fully understand kinematical observables. In particular, the gap and the uneven intensity contribution of the upper and lower surfaces of the disc can trigger artificial deviations from Keplerian rotation, which might appear similar to the perturbations on the gap edges (see Fig. 1, right panel) but should not be confused with these latter. For instance, from Fig. A.1 we notice that even after forcing the gas motions to be Keplerian, our line centroid method overestimates(underestimates) velocities on the gap at negative(positive) azimuths of the disc. As illustrated in Fig. A.2, this is due to the drop of intensity at the gap which abruptly alters the upper-to-lower surface intensity ratio, shifting the observed line centroids to non-Keplerian velocities. Another more direct observational consequence is that gaps alone are also able to produce kinks (see e.g. Fig. 2, middle rows) which may lead empirical methods to false positive inference of planets.

Fig. A.1 Line centroid (top), peak intensity (middle), and line width (bottom) residuals for the 1.0 MJup snapshot. The planet was removed and the gas velocity set to be fully Keplerian in order to analyse the contribution of the gas gap alone. The azimuthal scans on the left run along constant radii contours, whose colours represent their closeness to the gap, centred at R = 100 au. Although the velocities are Keplerian, the high (symmetric) centroid velocity residuals remain due to intensity differences between simulation and model at the location of the gap (see Fig. A.2). By simple comparison with Fig. 6, it is easy to identify the impact of the planet on the intensity and velocity residuals.

Fig. A.2 Comparison between simulation (dashed black) and model intensity profiles (dashed grey) extracted from different azimuths in the gap. Asymmetric line profiles are due to the contribution of both upper and lower emitting surfaces as indicated in the top right panel, whereas symmetric profiles are shaped bythe upper surface only. The solid black and grey lines are the Gaussian fits of the profiles and the vertical lines are the associated line centroids. Middle panel: centroid velocity residuals, in km s−1, zoomed-in on the central region of the Keplerian disc of Fig. A.1. The crosses are all at R = 100 au, and indicate the exact location of the intensity profiles in the surrounding panels. The gap was initially carved by a 1.0 MJup planet; thenthe planet was removed and the gas velocity set to be fully Keplerian for this analysis. The model, on the other hand, does not contain a gap, and so it systematically overestimates peak intensities there. Even though the disc rotation is Keplerian, there are high centroid velocity residuals due to differences between the simulation and model intensities. This effect appears all over the gap where the line profiles are asymmetric due to the contribution of the lower surface of the disc.

A proper model of the gap intensity would help reproduce the Keplerian pattern of the gap and in turn would facilitate the extraction of velocity fluctuations within it. Unfortunately, even a simple prescription of a gap implies at least three more model parameters (gap location, width, and depth), making our main goal of detecting planet-driven perturbations unnecessarily complex. Alternatively, one could simply study gas velocities at peak intensities (namely on the upper surface of the disc) only, instead of using full intensity profiles as our method does. However, this approach affects the accuracy of the retrieved velocities and omits any information coming from the lower surface of the disc, which we proved to be crucial for extracting velocity perturbations of planets at certain azimuths (see Sect. 5.3 and Fig. 12). Instead, we propose a simpler solution. The contribution of the gap can be readily excluded from the analysis by exploiting its symmetry around the projected minor axis, which we do in Sect. 4.3 by subtracting line centroids on the east from those on the west side of the disc.

In addition, we conducted an experiment to study the impact of non-planet-induced zonal flows on our detection analysis. To do this, instead of forcing the disc velocity field to be fully Keplerian, we smoothed it out by taking azimuthally averaged velocities across all pixels. This step eliminates any localised planet perturbation while keeping high-velocity fluctuations at the gap edges. We then applied the same analysis detailed in Sect. 4. Such a scenario effectively leads to non-detection of planets, with both azimuthal and radial peak velocity variances barely reaching a level of 1σ (see Fig. A.3). Again, this is because our line centroid folding procedure cancels any contribution from an azimuthally symmetric velocity field, regardless of its origin.

Fig. A.3 Same as Figs. 7 and 9 but for azimuthally symmetric velocity perturbations. Unlike planet-driven perturbations, such a scenario leads to non-localised velocity fluctuations and hence to non-detections according toour Variance Peak method.

Appendix B Supporting figures

Fig. B.1 Converging parameter walkers obtained with the DISCMINER for the 0.3 MJup snapshot, using 256 walkers and 1500 steps. This execution is preceded by an initial 1000-step run which is useful to find the seeding parameters before convergence. The dashed red line highlights the last quarter of walkers whose median corresponds to the reported best-fit parameters (in blue).

Fig. B.2 Peak velocity fluctuations against planet mass for different planet azimuths (colouredlines, left: negative, right: positive azimuths). A power-law fit of the form δυ = aMb is shown in colour according to each planet azimuth.

References

All Tables

All Figures

Fig. 1 Left panel: edge-on gas number density of one of our planet–disc simulation snapshots (1.0 MJup, ϕp = 0°). The white contour encloses the CO freeze-out region (T < 20 K). The black lines correspond to z∕R = 0.1, 0.2, 0.3, 0.4 scale heights, with the solid line being the threshold adopted for photodissociation. Middle panel: face-on view of the gas number density in the midplane of the disc. Right panel: azimuthal (main panel) and radial (zoom-in) deviations from Keplerian velocity. The solid and dashed contours in the zoomed-in panel are ±60 per cent peak azimuthal and peak radial perturbations, respectively, and illustrate that both components do not necessarily overlap and in turn contribute independently to the observed peak velocity perturbations. Also shown is the planet Hill sphere, with a radius of 6.8 au. The green circle indicates the current position of the planet, and the grey circles are additional planet azimuths explored in this work.

In the text

Fig. 2 Selected 12CO J = 2−1 synthetic channel maps for the three simulation snapshots (0.3, 1.0 and 3.0 MJup, from top to bottom), with the planets at Rp = 100 au radius, and ϕp = 45° azimuth, marked as green circles. The disc is inclined at − 45° with respect to the plane of the sky. The thin grey line is the projected circular orbit of the planet. The solid lines are centroid velocity contours extracted from the simulations at the velocity channel indicated on the top header; the dashed lines show the same but for the Keplerian best-fit model. Small arrows indicate kink-like features identified by visual inspection. In the bottom row a zoom-in around planets is shown for better comparison of centroid velocities and their deviation from Keplerian rotation as a function of planet mass.

In the text

	Fig. 3 Summary of the main attributes making up the line emission of a disc in the DISCMINER. Left panel: projected intensity of the disc for a channel centred on υch = − 1.0 km s−1. Right panel:line intensity profile extracted from the marker on the left. The grey annotations indicate the role of each attribute listed in Table 1.
In the text

	Fig. 4 Best-fit attributes obtained by the DISCMINER for the 0.3 MJup snapshot. Left panel: freeze-out and gas phase regions for 12CO. The coloured circles indicate the height of the model emitting surfaces. Right panels: model peak brightness temperature and line width as a function of radius on both emitting surfaces.
In the text

	Fig. 5 Peak line-of-sight deviations from Keplerian rotation versus positive and negative planet azimuths, for all three planet masses, extracted from the native simulations (without radiative transfer effects). The diamonds highlight turnover points where the radial velocity perturbations become important and eventually dominant over the projected azimuthal perturbations.
In the text

Fig. 6 Line centroid (top), peak intensity (middle), and line width (bottom) residuals for the (1.0 MJup, ϕp = 45°) snapshot. The azimuthal scans on the left run along constant radii contours in disc coordinates; their colours represent their closeness to the radial location of the planet, with blue being the closest. The solid black contour runs along the projected distance of the planet (Rp = 100 au) and the dashed black line shows the azimuth of the planet.

In the text

	Fig. 7 Folded (magnitude of) line centroid residuals as a function of azimuth, varying the analysis radius (in colours), for different planet masses (0.3, 1.0 and 3.0 MJup, from left to right) and azimuths (0°, 45°, and 90° from top to bottom). The solid black contour runs along the projected distance of the planet (Rp = 100 au), and the dashed black line corresponds to the azimuth of the planet.
In the text

Fig. 8 Left panel: peak residuals and their azimuthal location for the (1.0 MJup, ϕp = 45°) snapshot. Right panel: peak residuals rearranged as a function of their radial distance. The panel on the right shows the peak residual distribution (green line) and 1, 2, and 3σ significance thresholds (dotted lines). The solid black lines are the azimuthal and radial location of the global peak, while the dashed black lines show the actual location of the planet.

In the text

Fig. 9 Top row, left panel: KDE contours derived from K-means clusters of peak velocity residuals for the (1.0 MJup, ϕp = 45°) snapshot. The contours enclose 33 and 95 per cent of the peak residuals. The white crosses are the cluster centres retrieved by the K-means algorithm. Right panel: spectral variance of peak residuals for each cluster on the left. In the panel attached on the right we show 1, 2, and 3σ significance thresholds. The yellow line highlights the azimuthal location of the peak variance and the dashed line the actual location of the planet. The green cross is the peak variance significant enough to be attributed to the planet perturbation. Bottom row: 2D visualisation of the detection using the Variance Peak method for three planets (0.3, 1.0 and 3.0 MJup), all at ϕp = 45° azimuth and Rp = 100 au radius. The green circles and crosses are the actual and inferred location of the planets, respectively. The highlighted sectors are the azimuthal and radial clusters with the highest spectral variance extracted by the method. The boundaries of the sectors mark the maximum spatial coverage of the clusters. The circles show the location of the observed peak velocity residuals, whose magnitude is indicated by their size. The red circles are residuals within the azimuthal or radial peak variance clusters. The solid (δυϕ) and dashed (υR) contours trace spiral wakes and radial planet perturbations, and correspond to 60 per cent peak velocity fluctuations extracted from the simulation velocities.

In the text

	Fig. 10 Peak variance and peak velocity fluctuations extracted with the Variance and Global peak methods as a function of the detected planet azimuth for all three planet masses. The detected radial distance is shown in the bottom panels. Empty circles indicate the actual locations of planets.
In the text

	Fig. 11 Significance of planet detections using the Variance and Global peak methods as a function of planet azimuth for all three planet masses. The red line highlights the 3σ threshold considered to validate a detection. Empty circles indicate the actual location of planets.
In the text

Fig. 12 Illustrating how line centroids are shifted towards faster velocities around a 3.0 MJup planet, at two azimuths ϕp = ±30°. The crosses on the left indicate the location of the line profiles with the same colours on the right. The blue cross is the location of the observed peak perturbation, and the pink cross is the mirror pixel at the other side of the disc. We note that depending on the planet azimuth, the peak perturbation is projected either on the upper or on the lower emitting surface, which contribute differently to the observed line profiles. This contrast leads to variations in the observed line centroids (vertical black lines on the right) and the retrieved velocity fluctuations (in red).

In the text

Fig. A.1 Line centroid (top), peak intensity (middle), and line width (bottom) residuals for the 1.0 MJup snapshot. The planet was removed and the gas velocity set to be fully Keplerian in order to analyse the contribution of the gas gap alone. The azimuthal scans on the left run along constant radii contours, whose colours represent their closeness to the gap, centred at R = 100 au. Although the velocities are Keplerian, the high (symmetric) centroid velocity residuals remain due to intensity differences between simulation and model at the location of the gap (see Fig. A.2). By simple comparison with Fig. 6, it is easy to identify the impact of the planet on the intensity and velocity residuals.

In the text

Fig. A.2 Comparison between simulation (dashed black) and model intensity profiles (dashed grey) extracted from different azimuths in the gap. Asymmetric line profiles are due to the contribution of both upper and lower emitting surfaces as indicated in the top right panel, whereas symmetric profiles are shaped bythe upper surface only. The solid black and grey lines are the Gaussian fits of the profiles and the vertical lines are the associated line centroids. Middle panel: centroid velocity residuals, in km s−1, zoomed-in on the central region of the Keplerian disc of Fig. A.1. The crosses are all at R = 100 au, and indicate the exact location of the intensity profiles in the surrounding panels. The gap was initially carved by a 1.0 MJup planet; thenthe planet was removed and the gas velocity set to be fully Keplerian for this analysis. The model, on the other hand, does not contain a gap, and so it systematically overestimates peak intensities there. Even though the disc rotation is Keplerian, there are high centroid velocity residuals due to differences between the simulation and model intensities. This effect appears all over the gap where the line profiles are asymmetric due to the contribution of the lower surface of the disc.

In the text

	Fig. A.3 Same as Figs. 7 and 9 but for azimuthally symmetric velocity perturbations. Unlike planet-driven perturbations, such a scenario leads to non-localised velocity fluctuations and hence to non-detections according toour Variance Peak method.
In the text

	Fig. B.1 Converging parameter walkers obtained with the DISCMINER for the 0.3 MJup snapshot, using 256 walkers and 1500 steps. This execution is preceded by an initial 1000-step run which is useful to find the seeding parameters before convergence. The dashed red line highlights the last quarter of walkers whose median corresponds to the reported best-fit parameters (in blue).
In the text

	Fig. B.2 Peak velocity fluctuations against planet mass for different planet azimuths (colouredlines, left: negative, right: positive azimuths). A power-law fit of the form δυ = aMb is shown in colour according to each planet azimuth.
In the text