Kinematic evidence for an embedded planet in the IM Lupi disc

We test the hypothesis that an embedded giant planet in the IM Lupi protostellar disc can produce velocity kinks seen in CO line observations as well as the spiral arms seen in scattered light and continuum emission. We inject planets into 3D hydrodynamics simulations of IM Lupi, generating synthetic observations using Monte Carlo radiative transfer. We find that an embedded planet of 2-3 times the mass of Jupiter can reproduce non-Keplerian velocity perturbations, or `kinks', in the 12CO J=2-1 channel maps. Such a planet can also explain the spiral arms seen in 1.25mm dust continuum emission and 1.6 micron scattered light images. We show that the wake of the planet can be traced in the observed peak velocity map, which appears to closely follow the morphology expected from our simulations and from analytic models of planet-disc interaction.


INTRODUCTION
IM Lupi, a young star located 155.8±0.5 pc away (Klioner et al. 2018), hosts a large and spectacular protoplanetary disc (Panić et al. 2009;Cleeves et al. 2016;Avenhaus et al. 2018;Pinte et al. 2018a). Observations in the 12 CO J=2-1 spectral line from the Disk Substructures at High Angular Resolution Project (DSHARP) Andrews et al. 2018) showed localised deviations from Keplerian velocity in multiple velocity channels. More recently, the Molecules with ALMA at Planet-forming Scales Survey (MAPS) (Oberg et al. 2021) have confirmed these structures are present within the disc and are not observational artefacts. Pinte et al. (2020) predicted that a massive planet located at 117 au from the host star could be the cause of these deviations, referred to as kinks (see definition in Calcino et al. 2022). Detecting this planet using traditional observational methods, such as radial velocities, transits or direct imaging, is currently impossible as its orbit is too large and the planet is hidden within an optically thick disc of gas and dust (e.g. Pinte et al. 2018a).
Scattered light images of IM Lupi (Avenhaus et al. 2018) provide further evidence for embedded planets. This spectacular image revealed spiral structures in the upper layer of the disc traced by scattered light from sub-micron grains. Additionally, observations from DSHARP in the 1.25 mm dust continuum show spiral arms and gaps in the midplane dust disc .
IM Lupi is important because it is one of only a few discs that show large-scale spiral arms in mm-continuum emission. The other well-studied example is Elias 2-27 where there are two main hypotheses (Meru et al. 2017): an embedded planet (Meru et al. 2017) or gravitational instability (Forgan et al. 2018;Hall et al. 2018;Veronesi et al. 2021;Paneque-Carreño et al. 2021). The main prediction has been that gravitational instability will produce flocculent spiral structure and hence 'kinks everywhere' (Hall et al. 2020), whereas planets should produce more localised deviations from Keplerian motion (Pinte et al. 2020). Gravitational instability also requires a massive disc, while embedded planets do not require this.
Our aim in this Letter is to investigate whether an embedded planet in the disc can explain the observed substructures in IM Lupi. We model discs using hydrodynamical simulations, creating synthetic observations using Monte Carlo radiative transfer. Comparing our models to observations enables us to constrain the mass and location of the planet. We also predict kinematic perturbations from the planet wake in the peak velocity map, which we confirm are present in the observational data. Our paper is organised as follows: We describe our methods and initial conditions in Section 2, present our findings in Section 3, discuss the implications and limitations in Section 4 and conclude in Section 5 2. METHODS

Initial Conditions
We performed smoothed particle hydrodynamics (SPH) simulations of planet-disc interaction using Phantom . We modelled the disc with 10 7 SPH particles, set up initially to follow a tapered power-law surface density profile given by We assumed p = 0.48 (Pinte et al. 2018a). We used R c = 150 au instead of the Pinte et al. (2018a) value of R c = 284 au for computational convenience. We assumed a vertically isothermal equation of state P = c 2 s (R)ρ with c s ∝ R −q , q = 0.31 and the sound speed normalised to give an aspect ratio H/R = 0.129 at a radius of 100 au (Cleeves et al. 2016;Pinte et al. 2018a). We adopted a stellar mass of 1.12 M ) with the star modelled as a sink particle with an accretion radius of 1 au. We setup the disc initially between r in = 30 au and r out = 970 au (Panić et al. 2009), as we were not concerned with the inner disc structure. The vertically averaged ratio of smoothing length to disc scale height, < h > /H, varies between ≈ 1/20 at 200 au and 1/5 at 30 au, with around 10 resolution lengths per scale height at the final planet location and ≈ 5 at the height of the CO and scattered-light emitting layer in IM Lupi (Law et al. 2021). We adopted a total gas mass of 0.1 M , as determined by Cleeves et al. (2016). We also performed simulations with a disc mass of 0.01 M , but found that the higher disc mass better reproduces the scattered light image because the sub-micron sized grains remain well coupled in the top layers of the disc. We included dust in the simulations using the multigrain one-fluid algorithm Ballabio et al. 2018;Hutchison et al. 2018;Price et al. 2018). We modeled 11 grain sizes spanning a min = 1.0µm to a max = 2300µm on a logarithmically-spaced grid with a power law distribution with a slope of −3.5. We assumed a gas to dust ratio of 57 in order to give a total dust mass of 1.7 × 10 −3 M , as determined by Pinte et al. (2018a).

Embedded Planets
Pinte et al. (2020) predicted a planet orbiting IM Lupi with a semi-major axis of 117 au, based on the location of the velocity kink in the 12 CO channel in the DSHARP data. Since our planets migrate over the course of the simulation, we placed the planet initially at 150 au, which results in a planet at approximately the correct radius. The initial fast migration (10 au/kyr) occurs mainly because it takes several orbits for the planet to open a gap. After this the migration steadies to ≈ 2.3 au/kyr.
A consequence of the relatively high disc mass is that our planets, if modelled as 'vacuum cleaner' sink particles would accrete a large amount of mass (up to ten times their initial mass in the first fifty orbits). Since this is largely an artefact of the boundary condition employed (Ayliffe & Bate 2010), we instead modelled our planet as a non-accreting sink particle with a softened gravitational potential, similar to the procedure in Szulágyi et al. (2016). We assumed a cubic B-spline softening kernel (Price & Monaghan 2007) which becomes Keplerian at a radius 2h soft , where h soft is the softening length (see e.g. Price et al. 2018). This allows the sink particle mass for the planet to remain constant throughout the simulation, while producing the correct interactions with the surrounding gas and dust. We use a constant softening length of h soft = 1.75 au for the planet. To model the effects that the mass of the planet has on the disc, we performed four simulations, each with a single planet of mass 2, 3, 5 and 7 Jupiter-masses, respectively, plus one simulation with no planet. We kept all other initial conditions the same between each simulation.
We produced channel maps with a channel spacing of 0.05 km/s, from -2.4 km/s to +2.4 km/s. For the channel maps we assumed turbulent broadening of v turb = 0.08 km/s and f UV = 0.09 (Pinte et al. 2018a). We adjusted the azimuthal location of the planet based on Pinte et al. (2020), giving an azimuthal angle of 30 • (clockwise from North). To avoid CO being too bright in the outer disc we adopted a slightly lower (fixed) 12 COto-H 2 ratio of 1 × 10 −5 compared to 5 × 10 −5 assumed previously (Pinte et al. 2018a). This is likely due to our reduced value of R c .
We assumed a wavelength of 1.25 mm to produce the dust continuum image ) and a wavelength of 1.6 microns to produce the scattered light im-age (Avenhaus et al. 2018). For the scattered light comparison, we computed the polarized intensity (PI) image by combining the Q and U Stokes parameters (Pinte et al. 2006). We assumed spherical grains (Mie scattering) to compute the polarised emission, as described in Pinte et al. (2006)  the notional Stokes numbers from the hydrodynamical simulations (down) by a factor of 10. This mainly suggests that grains are not spherical (see Section 4). Figure 1 shows the observed 12 CO J=2-1 channel maps of the IM Lupi disc from the MAPS project (Oberg et al. 2021). A typical Keplerian butterfly pattern is visible. Distinct CO-emitting surfaces from the top and bottom of the disc are also evident (see Pinte et al. 2018a). On top of this butterfly pattern we indicate (with arrows) at 16 locations where 'kinks' or other perturbations away Keplerian motion are visible in the data. This seemingly argues against a planetary origin since previous papers suggested that an embedded planet should produce only localised velocity perturbations (Pinte et al. 2018b(Pinte et al. , 2020. Figure 2 shows the synthetic channels from our simulation with a 3M J planet. In all locations where a non-Keplerian signature was indicated in the data (arrows), we find a counterpart in the model which is significant in residuals from a no planet model (Appendix A). Essentially a kink occurs every time the planet wake crosses the channel, as predicted analytically by Bollati et al. (2021). The close match between the predicted and actual kinematic signatures in the channel maps, while not proof, shows that the hypothesis of a single embedded planet can explain the widespread non-Keplerian signatures in the disc without recourse to large scale instabilities. Our model replicates the upper and lower disc surface shape, extending out to the same approximate radius. Most interesting is that we reproduce the inner spirals in the SPHERE data. In our model these are caused by the propagation of the planet wake in the upper layers of the disc. The simulated spirals from the planet extend to the edge of the emitting surface, as in the observations. The dark line along the upper-left diagonal from the centre of the model image is caused by the phase function of the polarisation (Pinte et al. 2006). A drop in polarised intensity is also visible along the same diagonal in the observations.

Tracing the planet wake in the outer disc
A key prediction from our model is that the planetary wake should produce coherent velocity perturbations tracing the wake in the outer disc. That is, the non-Keplerian features seen in Figure 1 should trace the specific morphology of a planet wake, as opposed to being a series of flocculent spiral arms from gravitational instability. Calcino et al. (2022) found that the planet wake in the HD163296 disc was best traced in the peak velocity map. Figure 4 compares the observed peak velocity map computed from the CO cube for IM Lupi supplied by the MAPS team (Law et al. 2021; left) to same predicted from our 2M J model (right). To make the comparison clear, we overlaid both plots with a dotted line showing the analytic spiral wake (Rafikov 2002) best matching our simulation model, projected to the upper disc surface using the 12 CO emitting surface height calculated by Law et al. (2021) (following Pinte et al. 2018a. In both the model and the observations, we observe velocity perturbations that trace the planet wake through the outer disc. The distorted contours in the observational data appear to follow the planet wake both inside and outside the planet location; particularly in the second spiral. Particularly intriguing is the 'N-wave' structure  apparent in the distorted contours around the dashed line in the observational map, which are predicted by both our simulations (right panel) and by semi-analytic models of planet wake propagation (Rafikov 2002;Bollati et al. 2021).
Patches of red and blue seen in the bottom right of the data are simply where the lower surface of the disc becomes visible in the peak velocity map; this is reproduced in our model. Figure 5 shows the continuum emission (left column) and selected CO channels (right columns) from each of the four simulations to data from DSHARP and MAPS Czekala et al. 2021) (top row).

Continuum emission and planet mass estimate
The kink observed by Pinte et al. (2020) is located in the ∆v = −1.60km/s channel, which is also seen in our simulations. Some additional perturbations (or 'secondary kinks') can be seen on both the upper and lower surfaces, as in Figure 1. These correspond to velocity perturbations created by the planet wake.
Comparing the DSHARP continuum top panel) to the four simulations (bottom four panels), on all simulation models we observe an m = 2 spiral in the disc, increasing in amplitude with planet mass. For larger planet masses the simulations also show the formation of a gap in the outer ring, carved by the planet. However for planet masses ≥ 5M Jup the amplitude of the scattered light spirals are too high compared to the data (Figure 3), so we favour a planet mass of 2-3 M Jup .

DISCUSSION
In this Letter we explored whether an embedded planet can explain the spiral arms, non-Keplerian motions and scattered light emission in the IM Lupi disc. With caveats, we find that a single planet of several times the mass of Jupiter orbiting at ≈ 100-120 au can simultaneously explain all of these features. Our best evidence for a planet is in the kinematics. We confirm that the tentative kink found by Pinte et al. (2020) in the new data from MAPS (Oberg et al. 2021). However this is not the only kink in the data. In particular, a series of 'secondary kinks' (Bollati et al. 2021) appear to trace the planet wake in the observations in a manner predicted by our simulations. These are best seen in the observed peak velocity map (Figure 4), as recently demonstrated for the planet in the HD163296 disc (Pinte et al. 2018b;Calcino et al. 2022).
To reproduce the disc morphology in scattered light (with or without a planet) we needed the grains responsible for scattered light emission to be well-coupled to the gas, and for the disc to be optically thick at 1.6µm. In our hydrodynamic model we simulated a range of grain 'sizes' -assuming spherical grains and an intrinsic grain density of 3 g/cm 3 -from 1µm to 2.3 mm, corresponding to Stokes numbers of 10 −4 and 0.2 at 100 au, respectively (e.g. using Eq. 3 in Dipierro et al. 2015). In the radiative transfer calculation we needed an additional factor of ten scaling (down) in the Stokes number to avoid settling of the grains responsible for the scattered light emission out of the disc atmosphere. This is not surprising as Cleeves et al. (2016) already inferred that micron-sized grains have the same scale height as the gas in IM Lupi (h = 12au at r = 100au) and that the scale height of mm-emitting grains is ≈ 4× smaller. This likely indicates that grains are fluffy aggregates rather than compact spheres, as already suggested for IM Lupi by Pinte et al. (2008) and by ALMA polarization observations (Hull et al. 2018). Elias 2-27 is the other main disc that shows m = 2 spiral arms in continuum emission . The hypotheses are a planet, or gravitational instability (Meru et al. 2017;Forgan et al. 2018). Gravitational instability can lead to formation of multiple spiral arms and possibly fragmentation to form gas giants (Kimura & Tsuribe 2012;Kratter & Lodato 2016) and hence may also explain the main phenomena seen in IM Lupi. But one would need to posit that the inner disc is unstable (to produce spiral arms) while the outer disc, being relatively featureless in continuum, is not. This is opposite to the expected scenario where the outer disc is more unstable (Kratter & Lodato 2016). The global condition for instability is M disc /M * H/R (Toomre 1964). With M disc = 0.1 M the disc-to-star mass ratio is 0.1/1.12 = 0.09 which is smaller than the inferred H/R ≈ 0.12 at 100 au and hence suggests that the disc is stable. Sierra et al. (2021) estimated a Toomre Q ≈2 at r 50 au in IM Lupi, suggesting the disc is gravitationally stable. We neglected self-gravity in our models but the above discussion suggests that it could be important in IM Lupi, even if not ultimately responsible for the spiral arms.
Planets in protoplanetary discs migrate (Ward 1986). In Phantom, sink particles are free to interact with gas and dust, and hence migrate. We found that the 5M J planet migrated at ≈ 2 au/kyr after gap opening. Additionally, the migration depends on the planet mass, with more massive planets migrating faster. This can be seen in Figure 5 by comparing the grey circles (indicating the 117 au estimate) and the simulated location. Migration -which can be inward or outward depending on disc properties -makes it difficult to predict the radial location of the planet and brings a coincidence problem as to why any planet would be observed at a particular location. Nevertheless, we have shown that a planet in the right location can explain all the main substructures observed in IM Lupi. 5. CONCLUSIONS 1. A single embedded 2-3 M Jup planet orbiting at ≈110 au in the IM Lupi circumstellar disc can simultaneously explain the 16 different localised deviations from Keplerian motion ('kinks') seen in the 12 CO channel maps.
2. The same planet can explain the spiral arms seen in the upper layers of the disc in scattered light and the m = 2 spiral arms seen in continuum emission.
3. We predict, and confirm, that the wake from the planet should be visible in the observational peak velocity map from the 12 CO emission line. The perturbations seen in this map ( Fig. 4; left panel) closely match the prediction for the planetgenerated spiral arm We required a relatively massive disc (≈ 0.1M ) as well as a scaled Stokes number for the dust grains to remain well-coupled in both the upper layers and midplane of the disc. Such a disc mass suggests that gravitational instability may also be possible (but with Sierra et al. 2021 finding Q min ∼ 2). To assess the significance of the non-Keplerian features in the channel maps presented in Figures 1 and 2, Figure 6 shows the absolute brightness temperature residuals in the channel maps when the model with no planet is subtracted from the model shown in Figure 2. Arrows reproduced from Figures 1-2 all point to correspondingly significant features (∆T b > 10 K; corresponding to a signal-to-noise ratio of 7) in the residual maps. We found similar results when subtracting an azimuthally averaged model.