Planet migration, resonant locking and accretion streams in PDS70: Comparing models and data
Claudia Toci, Giuseppe Lodato, Valentin Christiaens, Davide Fedele, Christophe Pinte, Daniel J. Price, Leonardo Testi
MMNRAS , 1–14 (2020) Preprint 9 October 2020 Compiled using MNRAS L A TEX style file v3.0
Planet migration, resonant locking and accretion streams in PDS 70:Comparing models and data
Claudia Toci, , ★ Giuseppe Lodato, Valentin Christiaens, Davide Fedele, Christophe Pinte, Daniel J. Price, Leonardo Testi INAF OA Brera, via Brera 28, 20121, Milano MI, Italy Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano MI, Italy School of Physics and Astronomy, Monash University, VIC 3800, Australia INAF OA Arcetri, Largo Enrico Fermi 5, 50125 Firenze FI European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei Munchen, Germany
Accepted 2020 September 19. Received 2020 September 18; in original form 2020 August 10
ABSTRACT
The disc surrounding PDS 70, with two directly imaged embedded giant planets, is an ideal laboratory to study planet-discinteraction. We present three-dimensional smoothed particle hydrodynamics simulations of the system. In our simulations,planets, which are free to migrate and accrete mass, end up in a locked resonant configuration that is dynamically stable. Weshow that features observed at infrared (scattered light) and millimetre (thermal continuum) wavelengths are naturally explainedby the accretion stream onto the outer planet, without requiring a circumplanetary disc around planet c. We post-processed ournear-infrared synthetic images in order to account for observational biases known to affect high-contrast images. Our successfulreproduction of the observations indicates that planet-disc dynamical interactions alone are sufficient to explain the observationsof PDS 70.
Key words: planet-disc interaction – hydrodynamics — stars: individual (PDS 70)
Thanks to the exceptional capabilities of new observing facilities suchas the Atacama Large Millimeter Array (ALMA) and the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) instru-ment at the Very Large Telescope, high angular resolution observa-tions of protoplanetary discs are now available. These high-resolutionobservations have dramatically changed our view of disc evolution(Testi et al. 2015): small-scale substructures such as rings, gaps, cav-ities and spiral arms appear to be common (e.g. Benisty et al. 2017;Andrews et al. 2018). Interpreting these as signatures of planet-discinteractions (e.g., Dipierro et al. 2015; Pinilla et al. 2012; Fedeleet al. 2018; Dipierro et al. 2018; Pinte et al. 2018b; Zhang et al.2018; Lodato et al. 2019) suggests that giant planets may be alreadyfully formed at an age of only 1 Myr or less.This explanation is however not unique — other mechanisms havebeen invoked as the cause of the observed features, such as self-induced reconnection in magnetised disc-winds (Suriano et al. 2017),dust trapping (Gonzalez et al. 2015) and vortices (Barge et al. 2017).Planet signatures may appear in other disc diagnostics, such as in thekinematics of gas as traced by CO lines (Pinte et al. 2018b; Teagueet al. 2018; Pinte et al. 2019; Casassus & Pérez 2019). In some casesplanets detected with kinematical methods lie within observed dustgaps, strengthening the idea of planet-carved gaps.PDS 70 is a ∼ ± ★ E-mail: [email protected] the Scorpius-Centaurus association (Pecaut & Mamajek 2016). Thisyoung star is surrounded by a protoplanetary disc, first inferred frominfrared excess in the spectral energy distribution (Gregorio-Hetem& Hetem 2002) and lately directly imaged with different techniquesin the near infrared and sub-mm wavelengths (Hashimoto et al. 2012;Dong et al. 2012; Hashimoto et al. 2015; Long et al. 2018; Mülleret al. 2018; Keppler et al. 2019; Christiaens et al. 2019a; Isella et al.2019; Mesa et al. 2019). The disc has a wide gap (from ∼ −
50 au)clear of dust in the near infrared observations, that becomes larger( ∼
74 au) when observed at sub-mm wavelengths both in the dust(Keppler et al. 2018) and in the CO emission (Keppler et al. 2019).This suggests that dust grains of different sizes are spatially segre-gated by the radial pressure gradient of the disc, an hypothesis testednumerically in Bae et al. (2019).Two planets, PDS 70b and c, have been directly imaged in the disc(Keppler et al. 2018; Müller et al. 2018; Christiaens et al. 2019a; Haf-fert et al. 2019) at 195 and 234 milliarcsecond (mas), correspondingto deprojected distances of ∼
20 and 35 au, respectively. The lattersuggest that the planets might be in mean-motion resonance ( ∼ J ) for PDS 70b (Müller et al. 2018; Christiaens et al. 2019b; Haffertet al. 2019) and 2-10 M J for PDS 70c (Haffert et al. 2019; Bae et al.2019; Isella et al. 2019). Based on the excess emission of 𝐻𝛼 at thelocation of the protoplanet candidates, Haffert et al. (2019) estimateda mass accretion rate of ∼ × − M J yr − for the inner planetPDS 70b and ∼ − M J yr − for the outer planet PDS 70c, while © a r X i v : . [ a s t r o - ph . E P ] O c t C. Toci al. the stellar accretion rate is ∼ . × − M J yr − (Haffert et al.2019). However, these measurements are uncertain.The two planets can be responsible for the dynamical clearing ofthe inner disc and for the dust trap (Pinilla et al. 2012). Observa-tions in the sub-mm (Isella et al. 2019) and IR (Christiaens et al.2019b) wavelengths suggest the presence of a circumplanetary discsurrounding each planet, a possibility that has been theoreticallypredicted (see e.g., Lubow et al. 1999) and tested using numericalsimulations (e.g., Ayliffe & Bate 2009, Szulágyi et al. 2017 for a nu-merical study on the formation mechanisms of circumplanetary discsaround giant planets or Szulágyi et al. 2018 for numerical estimatesof the detection range of these structures using ALMA). Moreover, CO-integrated intensity maps (Keppler et al. 2019) shows evidenceof an inner decreasing of the flux density, located at ∼ ∼ 𝜇 m scat-tered light images, close to a broader and brighter feature (Keppleret al. 2018; Müller et al. 2018). The nature of these features is stilluncertain.Bae et al. (2019) performed two dimensional simulations of thissource. They showed that the presence of the inner planet alone(PDS 70b) fails to reproduce the observed peak in the radial profileof the sub-mm observation of (Keppler et al. 2019), as well as thespatially unresolved continuum emission close to PDS 70b of Isellaet al. (2019). They initialised their simulations with two grown giantplanets placed close to the 2:1 mean motion resonance, testing dif-ferent values for the mass of PDS 70c while fixing the mass of PDS70b to 5 M J , in order to find the best configuration to reproduce theobservations. They focussed on dust trapping, size segregation andfiltration. Bae et al. (2019) also discussed whether or not the observedunresolved continuum flux in the vicinity of the two planets (Isellaet al. 2019) could be explained with circumplanetary (CPD) discs,with their study providing an order of magnitude estimate of theirdust mass, ∼ . × − and 4 . × − Earth masses respectively.In this work we present new 3D smoothed particle hydrodynamics(SPH) simulations of the PDS 70 disc assuming two embedded giantplanets. We model both dust and gas, aiming to understand the phys-ical origin of features observed in mm and scattered light images. Incontrast to Bae et al. (2019), we start our simulations with lower massplanets and let them grow and migrate to eventually reach a meanmotion resonance. This condition better represents a real dynamicalscenario for planets embedded in the disc. In order to analyse howsensitive the final configuration is to the initial position and mass ofthe planets, we present the results of three simulations with differentinitial conditions.The paper is organised as follows: in Sec. 2 we describe the simu-lation initial conditions and radiative transfer models used to obtainsynthetic observations. A comparison between our models and obser-vational data and the dynamics of the planets are shown in Sec. 3. Wediscuss the implication of our results in Sec. 4 and give conclusionsin Sec. 5.
We perform a set of three dust and gas hydrodynamics simulationsin 3D using phantom, a smoothed particle hydrodynamics (SPH)code (Price et al. 2018). SPH (see Monaghan 2005 or Price 2012 for reviews) is particularly suitable for the study of accretion flowsand asymmetric interactions between planets and the disc, becausethere is no preferred geometry. Angular momentum is conserved tomachine precision, allowing for accurate orbital dynamics. phantomhas already been used widely for studies of dust and gas in discs (e.g.Dipierro et al. 2015; Ragusa et al. 2017; Cuello et al. 2018; Mentiplayet al. 2019; van der Plas et al. 2019; Veronesi et al. 2019; Toci et al.2020). In this work, we assume a single grain size of 1 .
00 mm,employing the one fluid algorithm (Laibe & Price 2014; Price &Laibe 2015; Ballabio et al. 2018) based on the terminal velocityapproximation (Youdin & Goodman 2005). In all our simulations,we include the dust back-reaction and the gravity between the gasand dust discs and the planets.
We present a small sample of three simulations where we fix thedisc parameters and we vary the initial conditions of two embeddedgiant planets. The system is modelled with a central star of massM ★ = . (cid:12) and the two giant planets embedded in a dust and gasdisc. The planets are allowed to migrate and accrete mass.In all the simulations our time unit is the period of the outer planet(Planet c) at its initial position in Sim 2 and 3, 𝑇 orb ( 𝑇 orb ∼
336 yr).All of our calculations have been evolved for at least 500 outer planetorbits, about ∼ We set up an initial surface density profile for the disc as in Lodato& Price 2010, representing the disc with 10 SPH particles. Thevalue of the initial inner and outer radii are and 𝑅 in = 𝑅 in = 𝑅 out =
120 au. The initial surface density profile for thegas is Σ ( 𝑅 ) = Σ (cid:18) 𝑅𝑅 (cid:19) − 𝑝 exp (cid:104) −( 𝑅 / 𝑅 ) ( − 𝑞 ) (cid:105) (cid:32) − √︂ 𝑅 in 𝑅 (cid:33) , (1)where 𝑅 =
100 au, the density normalisation Σ = .
07 g cm − is setto have an initial total gas mass of ∼ − M (cid:12) and we choose 𝑝 = 𝑐 s = 𝑐 s , (cid:18) 𝑅𝑅 (cid:19) − 𝑞 , (2)where 𝑞 = . 𝐻𝑅 = (cid:18) 𝐻 𝑅 (cid:19) (cid:18) 𝑅𝑅 (cid:19) ( / − 𝑞 ) , (3)where we set 𝐻 / 𝑅 = .
05 at 𝑅 = 𝑅 . In all our simulations we setthe value of 𝛼 SPH viscosity as in Lodato & Price 2010, in order tohave a value for the effective Shakura & Sunyaev 1973 viscosity of 𝛼 SS = . dust = − M (cid:12) . In all our dust simulation the chosensize of the dust is 𝑎 = 𝜌 d = − . For these values the initial midplane Stokesnumber, defined as 𝑆 𝑡 = 𝜌 d 𝑎 / Σ gas , is smaller than unity. This impliesthat dust and the gas are initially coupled in the disc. This justifiesthe use of the one fluid phantom mode (Laibe & Price 2014). MNRAS000
05 at 𝑅 = 𝑅 . In all our simulations we setthe value of 𝛼 SPH viscosity as in Lodato & Price 2010, in order tohave a value for the effective Shakura & Sunyaev 1973 viscosity of 𝛼 SS = . dust = − M (cid:12) . In all our dust simulation the chosensize of the dust is 𝑎 = 𝜌 d = − . For these values the initial midplane Stokesnumber, defined as 𝑆 𝑡 = 𝜌 d 𝑎 / Σ gas , is smaller than unity. This impliesthat dust and the gas are initially coupled in the disc. This justifiesthe use of the one fluid phantom mode (Laibe & Price 2014). MNRAS000 , 1–14 (2020) ynamical model of PDS 70 Table 1.
Initial conditions for the three simulations, called sim 1, sim 2, sim 3.Initial conditions for the masses are in Jupiter mass (M J ), the units for theinitial planet separations are au while the sink radii 𝑅 acc are measured infractions 𝑓 of the Hill radii (see sec 2.3 for the definition of Hill radius).SimulationsSim 1 Sim 2 Sim 3 𝑀 b ( 𝑀 J ) 2.5 3 3 𝑅 b (au) 20 15 20 𝑓 b = 𝑅 acc , b / 𝑅 H 𝑀 c ( 𝑀 J ) 0.5 0.5 0.5 𝑅 c (au) 55 45 45 𝑓 c = 𝑅 acc , c / 𝑅 H 𝑅 c / 𝑅 b The planets and the central star are included in the code as sinkparticles, characterized by an initial mass, position and sink radius.In phantom, sinks are free to migrate and are able to accrete gasand dust, interacting with dust and gas particle through gravitationalforces (for an overview of sink particles in SPH see Bate et al. 1995and Price et al. 2018). All the material that enters the sink radius isconsidered as instantaneously accreted on the sink . We chose thesink radius of each planet to be a fraction 𝑓 of the Hill radius 𝑅 𝐻 ,defined as 𝑅 H = (cid:18) 𝑀 p 𝑀 ★ (cid:19) / 𝑅 p , (4)where 𝑀 p and 𝑅 p are the planet mass and radial distance from the star.The values of 𝑓 for the two planets are different in the simulations,and are listed in Table 1. The value of the sink radius may affect theplanetary accretion rate and thus their final mass, so we decided tocheck the dependence of our results on this parameter.The initial conditions for the mass and the position of the planetsare significantly different from previous models (Bae et al. 2019).The fact that dust and gas can accrete onto the sinks represents acondition that mimics a real accretion scenario, and the fact thatplanets are free to migrate due to planet-disc interaction, planet-planet interaction and viscous disc evolution allow us to set up oursimulation without assuming that the planets are already close to themean motion resonance as inferred from observations. Our aim isindeed to test if planets that are initially massive enough to open agap in the gas surface density according to Crida et al. (2006) but areinitially smaller in masses and further out in position with respect tothe observational results can evolve due to accretion and migrationand obtain a final configuration with masses, position and motionresonances comparable with the observational results. The initialconditions for the planets can be found in Table 1. In all simulationsthe initial mass of the inner planet (M b ) is heavier than the one of theouter planet (M c ) and closer to the observed value, while the mass ofthe outer planet is far smaller than the observed value. This choice ismotivated by our previous study on discs with two embedded giantplanets (Toci et al. 2020): the initial mass of the outer planet has tobe small enough to allow it to migrate inward, while the higher massof the inner planet prevents this phenomenon. We also expect Planetc to accrete more mass than Planet b because of the larger gas anddust surface density in its neighbourhood. Technically instant accretion occurs at 80% of the sink radius by default,particles are accreted also at 80-100% the sink radius in special occasionssuch as if they are bound or infalling
To compare our simulations with observations, we post-processedour simulation using the Monte Carlo radiative transfer code mcfost(Pinte et al. 2006, 2009). mcfost maps the dust and gas densitiesfrom the SPH particles to the radiative transfer grid, using a Voronoitesselation to match each cell with an SPH particle, avoiding the needto interpolate.The dust distribution in each cell is obtained by assuming dustgrains smaller than 1 𝜇 m follow the gas, and by linear interpolationin log of the grain size between 1 𝜇 m and 1 mm. We assume a power-law grain size dust distribution 𝑑𝑛 / 𝑑𝑎 ∝ 𝑎 − . (integrating over thewhole disc), spanning the value in a range from 𝑎 𝑚𝑖𝑛 = . 𝜇 m and 𝑎 𝑚𝑎𝑥 = 𝜇 m. Grains larger than 1 mm are assumed to followthe distribution of the 1 mm grains. The dust optical properties arecomputed assuming spherical and homogeneous dust grains (Mietheory), with chemical composition of 60% astronomical silicatesand 15% amorphous carbons assumed, as the diana standard dustcomposition (Woitke et al. 2016), and a porosity of 10%. The totalgas mass is directly taken from our SPH simulation and the total dustmass is derived from the gas to ratio of the SPH simulation, which is100. The effective local dust to gas ratio for the various simulations isdiscussed in Sec. 3. We consider the emission of the central source asa black body with an effective temperature T e = 4090 K and a radiusof R ★ = 1.57 Solar radii, as derived from stellar evolution models fora 0.8 M (cid:12) at 5Myr (Siess et al. 2000).In all our simulations, we fixed the values of the source distance 𝑑 =
113 pc (Bailer-Jones et al. 2018), disc inclination 𝑖 = − ◦ ( following Müller et al. (2018) to have an observed clockwiseorbital motion of the planets 𝑖 > 90 ◦ ) and position angle PA=159 ◦ .We use 1 . × photon packets in the Monte Carlo simulation tocompute the temperature and specific intensities of our models. Wethen compute the source function by ray-tracing to produce the outputimages at the desired wavelengths. We produced 855 𝜇 m continuumimages and both 2.11 𝜇 m and 2.25 𝜇 m scattered light images ofthe t = 500 T orb ( ∼ CO moment 0 map, to test if our models canreproduce the flux value and the observed depth of the cavity in CO ofKeppler et al. (2019) observations. To compute the temperature andthe synthetic line maps, following Pinte et al. (2018b) and Pinte et al.(2019), we assumed local thermal equilibrium (LTE), 𝑇 gas = 𝑇 dust and a constant CO abundance of 10 − . We added a local spatiallyunresolved turbulent velocity of 50 m s − and we also consideredphoto-dissociation of molecules where the UV radiation is high.Precisely, as can be found in Pinte et al. (2018a), the condition forthe photo-dissociation of the CO is log ( 𝜒 / 𝑛 ) > −
6, where 𝑛 is thenumber density of hydrogen atoms and 𝜒 is the intensity of the UVradiation field. The intensity of the UV radiation field is evaluated asWoitke et al. (2009): 𝜒 = 𝐹 Draine ∫ 𝑛𝑚 . 𝑛𝑚 𝑢 𝜆 𝑑𝜆, (5)where 𝐹 Draine = . × − 𝑚 − 𝑠 − and where 𝑢 𝜆 is the energydensity of the radiation field computed by MCFOST in each point ofthe model. Finally, we convolved the sub-mm images with a Gaussianbeam to match the beam size of the observations (74 ×
57 mas forthe 855 𝜇 m image from Keppler et al. 2019. The IR synthetic imagesare further processed as explained in Section 2.5. MNRAS , 1–14 (2020)
C. Toci al.
The observation techniques and post-processing algorithms used toobtain high-contrast IR images of discs induce geometric biasesto the incoming signals. In particular, angular differential imaging(ADI; Marois et al. 2006) is known to induce significant distortionto extended disc signals (e.g. Milli et al. 2012; Christiaens et al.2019a). For a meaningful comparison of our synthetic IR imageto observations of PDS 70, we have both reduced the 2018-02-24 SPHERE/IRDIS dataset presented in Müller et al. (2018, ESOprogram 1100.C-0481), and post-processed our IR predictions in asimilar way to the observations.We calibrated the IRDIS data with a custom-made python pipelinewhich performs dark subtraction, flat-fielding, bad pixel correction,sky subtraction, star centering based on satellite spots, and gaussianfit of the unsaturated PSF in order to infer stellar flux and FWHM.The pipeline uses routines of the Vortex Image Processing package(VIP; Gomez Gonzalez et al. 2017), an open-source compilation ofpython routines used in high-contrast imaging. We then processedthe calibrated IRDIS cube using the median-ADI algorithm (Maroiset al. 2006) as in Müller et al. (2018).In order to mimic biases induced by high-contrast observations andimage processing, we applied the following effects to the MCFOSTimages of PDS 70 predicted at 2.11 and 2.25 𝜇 m (i.e. at the centralwavelength of the IRDIS K1 and K2 filters, respectively):(i) Injection of the flux of protoplanets PDS 70 b and c at thelocation of the corresponding sink particles, using the contrast withrespect to the star measured in Müller et al. (2018); Mesa et al.(2019);(ii) Convolution with the non-coronagraphic normalized pointspread function observed in the IRDIS 2018-02-24 dataset, whichhas a FWHM of ∼ ∼ Figure 1 shows the gas and dust (1 mm) surface density distributionfor the t = 500 T orb ( ∼ -1 -0.5 0log (cid:1) (cid:2) g [g/cm ] dz [cm]100 au (cid:1) (cid:2) d [g/cm ] dz [cm]100 au Figure 1.
Top: logarithmic plot of the gas surface density map of thet = 500 T orb ( ∼ orb ( ∼ is trapped in a ring due to the combined effect of dust trapping inthe pressure maximum (Pinilla et al. 2012) and radial drift. This isconsistent with the results of Bae et al. (2019), where only smaller( < 𝜇𝑚 ) dust grains can filtrate inside the cavity. As visible in Fig.1,an accretion stream connects Planet c with the gas rim of the disc,highlighting that the planet is actively interacting with the disc.In our simulations, a circumplanetary disc is created around allthe planets and quickly (in about 50-100 T orb ) accreted. In Figure 1,shown after 500 T orb , only gas around the planets can be spotted.Figure 2 shows how the synthetic dust grain population is obtainedfrom the simulation during the post-processing with MCFOST. InFig. 2 (top) the surface density profiles of gas (light blue) and dust(blue) are plotted as a function of the radius for the above snapshot, inFig. 2 (middle) we show the density profiles of the two synthetic dustpopulations (large and small grains) and of the total dust distribution.Finally, Fig. 2 (bottom) displays the local value of the dust-to-gas ratioany given radial point. MNRAS000
Top: logarithmic plot of the gas surface density map of thet = 500 T orb ( ∼ orb ( ∼ is trapped in a ring due to the combined effect of dust trapping inthe pressure maximum (Pinilla et al. 2012) and radial drift. This isconsistent with the results of Bae et al. (2019), where only smaller( < 𝜇𝑚 ) dust grains can filtrate inside the cavity. As visible in Fig.1,an accretion stream connects Planet c with the gas rim of the disc,highlighting that the planet is actively interacting with the disc.In our simulations, a circumplanetary disc is created around allthe planets and quickly (in about 50-100 T orb ) accreted. In Figure 1,shown after 500 T orb , only gas around the planets can be spotted.Figure 2 shows how the synthetic dust grain population is obtainedfrom the simulation during the post-processing with MCFOST. InFig. 2 (top) the surface density profiles of gas (light blue) and dust(blue) are plotted as a function of the radius for the above snapshot, inFig. 2 (middle) we show the density profiles of the two synthetic dustpopulations (large and small grains) and of the total dust distribution.Finally, Fig. 2 (bottom) displays the local value of the dust-to-gas ratioany given radial point. MNRAS000 , 1–14 (2020) ynamical model of PDS 70 [ g c m ] dust [g cm ] gas [g cm ] d u s t c o m p o s i t i o n dust , large dust , small dust , tot d u s t / g a s r a t i o Figure 2.
Top: surface density profiles of gas (light blue) and dust (blue)plotted as a function of the radius for the t = 500 T orb ( ∼ Figure 3 shows the surface density profile of gas (light blue) anddust (blue) plotted as a function of the radius for the t = 500 T orb snapshots (Panels a). For each simulation, we also show the simulateddust continuum image at 855 𝜇 m convolved with a 74 ×
57 masGaussian beam (Panels b), and the average between the 2.11 𝜇 mand 2.25 𝜇 m scattered light image, either obtained after the steps(i) to (ix) detailed in Section 2.5 (Panels c), or only convolved witha PSF of 0.05 mas (Panels d). The latter aim to show the authenticmorphology of scattered light signals before ADI processing. Table 2 Table 2.
Masses and positions of the two planets, gas and dust disc massesfor our set of simulations at t=500 T orb . Planets’ masses are in Jupiter mass(M J ), distances are in au and the disc masses are in Solar masses.SimulationsSim 1 Sim 2 Sim 3 𝑀 b ( 𝑀 J ) 4.4 5. 5.3 𝑅 b (au) 19.4 14.2 19.1 𝑀 c (M J ) 4.1 4.3 3.9 𝑅 c (au) 42.6 34.9 36.1 𝑀 gas (M (cid:12) ) 2.5 × − × − × − 𝑀 dust (M (cid:12) ) 5.4 × − × − × − lists the masses and positions of the planets as well as the disc massat the time of the snapshots.All our models reproduce the observed morphology of the dustand gas rings, as well as both the small features that are found inthe images (the dust skewness seen from Isella et al. 2019 and Kep-pler et al. 2019 and most of the features seen in the cavity in the2.11/2.15 𝜇 m SPHERE images from Müller et al. 2018).Focusing on the surface density profiles (Panels a in Fig. 3), theposition of the inner Planet b (plotted as an orange vertical lines)does not appear to impact the shape of the surface density profile:all the simulations show a large and deep cavity in the inner — ∼
50 au — part of the disc (the surface density in the inner part is atleast two orders of magnitude lower than the outer part of the disc),regardless of the final positions of Planet b, that migrates very littlein all cases. On the contrary, the final position of the outer planet(Planet c, shown as a yellow vertical line) affects both the width ofthe gas cavity (larger for larger values of the initial Planet c distance,as in Sim 1) and the position of the dust ring (that appears at ∼ ∼ 𝜇 m image of Sim 1 has a fainter andlarger ring compared to the images of Sim 2 and 3 (see panels b inFigure 3). However, a spur is always present, in the same position,in all cases. This is in agreement with what observed with ALMA(Isella et al. 2019). This feature corresponds exactly to the locationof Planet c, and is due to the interaction between the planet and theouter part of the disc. Indeed, the planet is actively accreting materialfrom the dust and the gas ring that may form a circumplanetary dustring. This feature has already been reproduced and analysed in (Baeet al. 2019).Panels c and d of Figure 3 show the synthetic 2.11/2.15 𝜇 m imagesof our simulations with and without median-ADI post-processing,respectively. Contrary to panels d, panels c were obtained with theflux of planets b and c injected in the synthetic ADI cube at thelocation of the sink particles. This is for a fairer comparison to theobservations, given that the bright flux of the planets also affectsthe median-ADI algorithm. Sim 1 can be seen to produce a largercavity and a broader ring compared to the other two cases, as alreadydiscussed. However, we stress that we naturally explain most of theobserved bright features in the scattered light image without anyspecific assumption. In order to interpret the origin of these features,labelled i to iv in Figure 3, we also processed a snapshot of thesimulation corresponding to planet c located at the opposite point ofits orbit (+180 degrees true anomaly) compared to the observations.This is shown in Figure 4 for Sim 3.The dust spur (seen in the 855 𝜇 m image, Panels b) is now presenton the other side of the disc, confirming that this is not a numericalfeature but is due to the presence of dust surrounding Planet c.However, some of the intra-cavity emission features in the 2.11/2.15 MNRAS , 1–14 (2020)
C. Toci al.
Keppler et al. 2019 Muller et al. 2018
Sim1 Sim 2 Sim 3 c i ii b iii iv
Figure 3.
Comparison between our three simulations (bottom three rows) and observations from Keppler et al. (2019) and Müller et al. (2018) (top row). Panelsa show the surface density profile of gas (light blue) and dust (blue) plotted as a function of radius for the t = 500 T orb snapshot for the three simulations (1 to3 from top to bottom), in units g cm − and au respectively. Instantaneous planets positions are also plotted as red and orange vertical lines. Panels b show thecorresponding simulated continuum image at 855 𝜇 m, with fluxes in Jy beam − . The image has been smoothed with a 74 ×
53 mas Gaussian beam to match thebeam size of (Keppler et al. 2019). Panels c and d show the average between the 2.11/2.15 𝜇 m scattered light images with and without median-ADI processing,respectively. In panels c the flux of the planets were injected in the ADI cube before processing, while panels d do not include emission/heating from the planets.The flux density is in Jy/arcsec − and the image has been convolved with a psf of 0.05 mas. 𝜇 m image (panels c) are still present at roughly the same location.Features i to iv are interpreted as follow:(i) Feature i is an artefact of median-ADI processing of the brightforward-scattered edge of the cavity. Similar artefacts were obtainedin the forward modelling tests shown in Milli et al. (2012).(ii) Feature ii likely consists of two contributions: scattered light signal from small intra-cavity dust (as suggested from the blob at thesame location in the bottom panel of Figure 4) and from the spiralaccretion stream extending inward from c. The latter is confirmedby the residual blob seen at the location of feature ii in the imageobtained when subtracting the scattered light snapshot where planet MNRAS000
53 mas Gaussian beam to match thebeam size of (Keppler et al. 2019). Panels c and d show the average between the 2.11/2.15 𝜇 m scattered light images with and without median-ADI processing,respectively. In panels c the flux of the planets were injected in the ADI cube before processing, while panels d do not include emission/heating from the planets.The flux density is in Jy/arcsec − and the image has been convolved with a psf of 0.05 mas. 𝜇 m image (panels c) are still present at roughly the same location.Features i to iv are interpreted as follow:(i) Feature i is an artefact of median-ADI processing of the brightforward-scattered edge of the cavity. Similar artefacts were obtainedin the forward modelling tests shown in Milli et al. (2012).(ii) Feature ii likely consists of two contributions: scattered light signal from small intra-cavity dust (as suggested from the blob at thesame location in the bottom panel of Figure 4) and from the spiralaccretion stream extending inward from c. The latter is confirmedby the residual blob seen at the location of feature ii in the imageobtained when subtracting the scattered light snapshot where planet MNRAS000 , 1–14 (2020) ynamical model of PDS 70 RA ["] D e c [ " ] Continuum 855 m
RA ["] D e c [ " ] Scattered light 2.11/2.25 m
RA ["] D e c [ " ] Classical ADI K K F l u x d e n s i t y [ m J y . b e a m ] Figure 4.
Simulated continuum image at 855 𝜇 m (top) and the average be-tween 2.11/2.15 𝜇 m scattered light images (bottom) for the t = 500 T orb snapshot of sim 3 when planet c is located 180deg forward in true anomaly,compared to the observed one. While the forward-scattered intra-cavity dustsignal remains in the scattered light image, the spiral accretion stream, in-cluding the bright blob surrounding c, disappears. c is 180deg away from its observed location from the snapshot whereit is at its observed location (Appendix A).(iii) Feature iii could trace the back-scattered emission from smalldust grains in the cavity, as a faint signal is also seen on that side ofthe cavity in our predicted images.(iv) Feature iv is not reproduced in our processed RT images. It is unclear whether this point-like feature discussed in Mesa et al.(2019) could trace the edge of a denser inner disc – not present inour simulations – or whether it is a residual speckle.A more detailed discussion of post-processing artefacts includ-ing a comparison between polarimetric differential imaging (PDI),median-ADI and principal component analysis (PCA) ADI observa-tions on the one hand, and our predicted processed images in theother is provided in Appendix B.In all our images the dusty feature observed around the inner planet(Planet b) in Isella et al. 2019 is absent. This does not mean that adusty circumplanetary disc should not be present around this planet.Indeed, in our simulations dust material surrounding this planet formsa disc that is quickly accreted on the sink due to numerical effects.For a more complete discussion on this feature see Bae et al. 2019. Figures 5 and 6 illustrate the dynamics of the planets for the threesimulations. In particular, Figure 5 shows the position of Planet b(left), Planet c (middle) and the ratio between the two orbital periods(right), to check whether or not they settle into a mean motion reso-nance, as a function of the normalized time. In the first ∼
50 orbits thesystems relaxes out of the initial conditions: Planet c quickly migrates ∼
10 au inward (middle panel), while the inner planet (Planet b) re-mains closer to its initial position, due to its larger initial mass (leftpanel). The migration of the outer planet ends when the orbits of thetwo objects lock in a period resonance (right panel) that is differentin the various simulations. This value appears to be dependent on theinitial position of the outer planet (Planet c): a larger initial distancebetween the star and Planet c (Sim 1) leads to an equilibrium con-figuration when Planet c more distant compared with the other cases(Sim 2 and 3). No significant dependence on the initial position ofthe inner planet (Planet b) is found (see the behaviour of Sim 2 andSim 3). The value of the ratio between the orbital periods dependsalso on the position of the inner planet, Planet b. In order to matchthe 2:1 expected resonance for this system as well as the ALMA ob-served 855 continuum flux, Planet b should be located at ∼ ∼ ∼
100 T orb
Planet b has already accreted almost all ofthe available mass, reaching the value of 4.5 - 5 M J , depending on theinitial conditions, while the outer planet continues to accrete massfrom the external (R >
50 au) reservoir of gas, reaching ∼ J massat the end of the simulation, with a significant mass accretion rate.This suggests that at the end of the evolutionary path of the disc, themass of the outer planet could be even larger. In all our models theorbital eccentricity of the planets remains relatively small, less thanabout 0.1 for Planet c (with an exact value dependent on the planetmass) and ∼ ∼
100 orbits but it is quickly accreted about ∼ MNRAS , 1–14 (2020)
C. Toci al. orb a b [ a u ] Sim 1Sim 2Sim 3 orb a c [ a u ] Sim 1Sim 2Sim 3 t / T orb T c / T b Sim 1Sim 2Sim 3
Figure 5.
Planet positions for the simulations: Semi-major axis of Planet b (left), Planet c (middle) and the ratio between the two orbital periods (right) plottedas a function of the normalized time for the three simulations (see legend). The period ratio indicates that planets reach resonances like 5:2 and 7:2. orb m b [ M J Sim 1Sim 2Sim 3 orb m c [ M J ] Sim 1Sim 2Sim 3 t / T orb m b / m c Sim 1Sim 2Sim 3
Figure 6.
Masses of Planet b (left) and Planet c (middle) as a function of the normalized time, and their relative ratio (right) for the three simulations.
Figure 7.
CO J=3-2 integrated intensity (moment 0, continuum subtracted) for all the simulations for the t = 500 T orb snapshot and, as a reference, the observedimage of Keppler et al. (2019) (left). The position of the star and the planets is shown as cyan dots. The image is convolved with a Gaussian 76 ×
61 mas beamwith PA of 63 degrees as in Keppler et al. (2019) observations.
In order to test whether our models are compatible with CO observa-tions of Keppler et al. (2019) we also produced synthetic CO momentzero maps. Figure 7 shows the synthetic emission in CO J = 3 - 2of our three models (last three panels) compared to the observations(left panel). The two planets open a ∼
40 au wide cavity with a de-crease in surface density of at least two orders of magnitude in thegas (see Fig.3), with an emptier inner region (R < 25 au). This resultsin the CO moment zero map in an inner gap with no emission (all the gas is photo-dissociated) and the presence of an asymmetric ringat ∼ MNRAS000
40 au wide cavity with a de-crease in surface density of at least two orders of magnitude in thegas (see Fig.3), with an emptier inner region (R < 25 au). This resultsin the CO moment zero map in an inner gap with no emission (all the gas is photo-dissociated) and the presence of an asymmetric ringat ∼ MNRAS000 , 1–14 (2020) ynamical model of PDS 70 F l u x [ m J y / B e a m k m / s ] Azimuthally averaged flux CO Sim 1Planet bPlanet ccontinuum F l u x [ m J y / B e a m k m / s ] Azimuthally averaged flux CO Sim 2Planet bPlanet ccontinuum F l u x [ m J y / B e a m k m / s ] Azimuthally averaged flux CO Sim 3Planet bPlanet ccontinuum
Figure 8.
Azimuthally averaged CO moment zero maps (continuum subtracted) for for all the simulations (left: Sim 1, middle: Sim 2, right: Sim 3) for thet = 500 T orb snapshot. The position of Planet b and c is shown as the red and orange vertical line respectively, while the grey region shows the extent of thecontinuum ring. neighbourhood of the planets’ orbits. We also find a change in theslope of the profiles located at 0.6 arcsec. Comparing our imageswith Keppler et al. (2019) observations, we observe that the the ob-served CO images peaks at about 0.4 arcsec while the models showthe peak at about 0.3 arcsec. This difference is probably due to ahigher gas mass inside the cavities in our models with respect to thesource, accentuated by the absence of an inner ring that can shieldthe material from the radiation of the star. Higher resolution models,or simply evolving the simulations for longer may help to resolve thisdiscrepancy.Figure 9 shows the mass accretion rates (in M (cid:12) yr − ) onto theplanets in our simulations (left: Sim 1, middle: Sim 2, right: Sim3) as a function of normalized time. We recall that by definition inSPH all the material accreted by the sink is merely added to theplanet mass. In reality, part of the “accreted mass” may end up inan unresolved circumplanetary disc that buffers the material onto theplanet, so that our computed accretion rates may be considered asupper limits.All three simulations show the same behaviour: the accretion rateon the star is ∼ × − M (cid:12) yr − , an order of magnitude higher thanthe observed value. The fact that our sink radius for the star (2 au) isfar larger than the actual radius for the star can be the reason for thisoverestimate. However, after relaxing out the initial conditions, thisvalue slowly decreases with time, suggesting that gas is continuouslybeing accreted on the star. Our simulations also show that the outerPlanet c has a higher accretion rate ( ∼ − M (cid:12) yr − ) with respectto the inner Planet b, ( ≤ − M (cid:12) yr − ). This occurs naturally inour simulations because the outer planet is directly in contact withthe external mass reservoir, in contrast to the inner planet, that sitswithin the cavity.Moreover, while the accretion rate of Planet c is (almost) constantwith time, because it is actively accreting material (as the morpho-logical features — seen in dust and gas — emphasize), the accretionrate of Planet b rapidly drops as the cavity is depleted, "starving" theinner planet. This occurs because the outer giant planet prevents thereplenishment of the inner disc.The accretion rate we obtain for Planet b is consistent with theobserved value of ∼ × − 𝑀 (cid:12) yr − . By contrast, we appearto strongly overestimate the accretion rate onto Planet c (which is × − 𝑀 (cid:12) yr − ). As mentioned above, we stress that the “accretionrate” we measure is actually a measure of how much mass enters theHill’s sphere of the planet. Most likely, a circumplanetary disc would form (which would be unresolved in our simulations) which wouldact as a buffer for the accretion flow onto the planet. We also neglectany radiative feedback from the planets. Without a proper modellingof such disc, we cannot draw any firm conclusions regarding theplanetary accretion rates. However, regardless of the absolute valuesof the accretion rates, their relative ordering, with the outer planetaccreting at a higher rate than the inner planet, should be preservedand is a natural consequence of the inner disc depletion. Our simulations reproduce the main features found in ALMA 855 𝜇 m and VLT/SPHERE 2.11/2.15 𝜇 m observations (see Fig. 3). Inparticular, the complex morphology observed in scattered light atthe inner edge of the ring is explained by a combination of post-processing biases and the high inclination of the disc, with a likelyminor contribution from a spiral accretion stream (Section 3.1 andAppendix A).We find that the best scenario (Sim 3) gives final values for theplanet masses of ∼ −
10 M J , with Planet b located at ∼ ∼ −
40 au, in a 5:2 mean motion resonance,in agreement with observational results (Keppler et al. 2019) andprevious models (Bae et al. 2019).The initial mass of the inner planet (Planet b, 2.5-3 M J ) is lighterthan that of the outer planet (Planet c, 0.5 M J ). This could eitherrepresent a situation where the two planets are born at a comparabletime with very different masses or the effectively equivalent casewhere the inner planet forms earlier, migrates and reaches a highermass at the time where the second planet forms in the outer disc.We show that giant planets can reach an orbital resonance while theystill are embedded in the disc and still accreting gas, supporting Baeet al. (2019)’s findings for PDS 70. The actual value of the resonanceappears to depend on the disc and planet masses, but further studiesare necessary to better characterize this phenomenon. This behaviourhas also been proposed for Jupiter and Saturn in the Solar System(see e.g Masset & Snellgrove 2001), with the caveat that PDS 70 andour Solar System are very different due to the huge differences in thegiant planet masses and positions.The presence of two giants planets on stable orbits sculpts thehost disc, opening a wide cavity in the dust and gas density andblocking the gas reservoir and the mm dust ring in the outer partof the disc. Filtration mechanisms could be at play, as found in Bae MNRAS , 1–14 (2020) C. Toci al. orb M [ M y r ] StarPlanet bPlanet c orb StarPlanet bPlanet c orb StarPlanet bPlanet c
Figure 9.
Accretion rates (in M (cid:12) yr − ) onto the sinks for all the simulations (left: Sim 1, middle: Sim 2, right: Sim 3) as a function of the normalized time.Orange lines are the star accretion rate and red and pink lines are Planet b and Planet c accretion rates respectively. Horizontal red and pink dashed lines showthe inferred values from Planet b and c respectively from Haffert et al. (2019), while the horizontal red and pink dotted lines show the inferred values fromPlanet b and c from Hashimoto et al. (2020). et al. (2019): large grains are trapped beyond their orbits, whilesmaller grains can filtrate, changing the amount of solids availablein the inner and outer part of the disc, impacting on the formationof terrestrial planets and possibly on the composition of the coresof the giants. Indeed, this can eventually lead to the formation oftwo different and spatially separated mixtures of dust and gas, thatmay have influences on the future planetary formation and chemicalcomposition. Focusing on the gas, the flux values of the synthetic COmoment zero maps are compatible with the observations of Keppleret al. (2019), however we fail to reproduce the compact inner ringfound in the observations, probably because of our use of a 2 ausink radius for the star. Considering the azimuthally averaged fluxprofiles, the ring we find has the correct flux value but is locatedcloser to the star (0.3 arcsec in our models and 0.4 arcsec in theobservations), pointing out that probably there is too much gas in thecavity. However, finite resolution effects could be at play, due to thefact that the planets are actively emptying the cavity, decreasing alsothe SPH resolution.When two giant planets in resonance shape their host disc, therelative distance between their positions appears to be particularlyrelevant to determine the dust morphology. Indeed, in the case ofPDS 70, where the orbital ratio is close to 2:1 no dust ring is ob-served between the orbits of the planets. By contrast, in the caseof HD 169142, another source where at least two candidate giantsprotoplanets could be present (Fedele et al. 2017; Pérez et al. 2019),a prominent dust ring is observed between the likely position of thetwo giants, that are expected to be in a larger orbital ratio. Numericalresults (Toci et al. 2020) confirm the presence of a long lived dustyring in this case. Apparently, whether or not a dust ring between thetwo planets is present depends on the order of the mean-motion res-onance, with closer orbits not allowing a long-lived inter-planetarydust ring.Another result from our simulations is that the relative positionbetween the planets and the dust or gas material in the disc plays akey role in inferring their presence. While ultimately only a directmeasurement of a planet can confirm its presence, detection of ac-cretion streams or sub-mm circumplanetary discs can help pinpointthe candidate planet location.It is unclear whether the observed mm-continuum emission near c(refereed to as “spur”) corresponds to circumplanetary disc emission,or whether it traces the wake from c passing through that location. Thelatter interpretation would be consistent with all the non-detections of circumplanetary discs (e.g. Ricci et al. 2017), the sub-mm imageof PDS 70 shown in Isella et al. (2019) after subtraction of thesymmetric component (showing a stream-like residual), and the sub-mm spiral-like feature seen within an annular gap in the HD 163296disc (Isella et al. 2018). The shift between the IR and sub-mm signalsnear b could possibly be explained in a similar way.One of our aims was to test if the results of our simulations aresensitive to the initial conditions for the planets. In the small sampleshown, the initial positions of the planets, as well as their initialmasses, appear not to have a large impact on the global shape ofthe disc. However, a more detailed study of the parameter space isneeded to understand the robustness of the modelOur simulations did not allow us to firmly constrain the planetaryaccretion rates. Numerical effects, due to a high numerical viscos-ity given by low resolution of the circumplanetary disc around thesink particles, would lead to a shorter life-time for the circumplane-tary disc, drastically increasing the mass accretion rate. Particularlysensitive to these effect is the case when the sink is continuouslyinteracting with a reservoir of material, as in the case of Planet c.However, we do expect that the relative ordering of the two plan-etary accretion rates found in our simulations should be preservedat higher resolution, with the outer planet accreting at a higher rate.Current estimates of the planetary accretion rates in PDS 70 (Haf-fert et al. 2019; Hashimoto et al. 2020) appear to go in the oppositedirection, although these are affected by large uncertainties. In this work we performed 3D SPH simulations of the PDS 70 pro-toplanetary disc. Because two giant planets have been directly ob-served in this source, we initialised our simulations with two planetsembedded in a disc of dust and gas, with smaller initial mass andlarger orbital ratio than observed. In contrast to previous modellingattempts, we allowed our planets to grow in mass and to migratewithin the disc. Our simulation setup mimics a scenario in which theinner and the outer planets form and grow at different locations in thedisc and are captured in a mean-motion resonance while migratingand accreting mass. We found that the two planets quickly lock intoresonance with a ratio that is compatible with observational results aswell as with previous 2D numerical simulations (Bae et al. 2019). Wealso found excellent agreement between our mm and scattered-lightimages and the literature, reproducing the main observed features
MNRAS000
MNRAS000 , 1–14 (2020) ynamical model of PDS 70 caused either by the presence of the outer planet or by geometricbiases related to high-contrast IR observations and post-processing.Planet-disc and planet-planet interactions are at play in PDS 70,carving the disc in the dust and gas reservoir. Planet migration andthe resonant locking of the orbits are key processes. Millimetre-emitting dust grains trapped in a gas pressure bump at ∼
75 au causethe observed ring seen in continuum emission, supporting the resultsfrom Bae et al. 2019.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
ACKNOWLEDGEMENTS
REFERENCES
Absil O., et al., 2013, A&A, 559, L12Amara A., Quanz S. P., 2012, MNRAS, 427, 948Andrews S. M., et al., 2018, ApJ, 869, L41Avenhaus H., et al., 2017, AJ, 154, 33Ayliffe B. A., Bate M. R., 2009, MNRAS, 397, 657Bae J., et al., 2019, ApJ, 884, L41Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Mantelet G., Andrae R.,2018, AJ, 156, 58Ballabio G., Dipierro G., Veronesi B., Lodato G., Hutchison M., Laibe G.,Price D. J., 2018, MNRAS, 477, 2766Barge P., Ricci L., Carilli C. L., Previn-Ratnasingam R., 2017, A&A, 605,A122Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362Benisty M., et al., 2015, A&A, 578, L6 Benisty M., et al., 2017, A&A, 597, A42Calcino J., Christiaens V., Price D. J., Pinte C., Davis T. M., van der MarelN., Cuello N., 2020, arXiv e-prints, p. arXiv:2007.06155Casassus S., Pérez S., 2019, ApJ, 883, L41Christiaens V., et al., 2019a, MNRAS, 486, 5819Christiaens V., Cantalloube F., Casassus S., Price D. J., Absil O., Pinte C.,Girard J., Montesinos M., 2019b, ApJ, 877, L33Crida A., Morbidelli A., Masset F., 2006, Icarus, 181, 587Cuello N., et al., 2018, Monthly Notices of the Royal Astronomical Society,483, 4114–4139Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015,MNRAS, 453, L73Dipierro G., et al., 2018, MNRAS, 475, 5296Dong R., et al., 2012, ApJ, 760, 111Dong R., Zhu Z., Rafikov R. R., Stone J. M., 2015, ApJ, 809, L5Fedele D., et al., 2017, A&A, 600, A72Fedele D., et al., 2018, Astronomy & Astrophysics, 610, A24Gomez Gonzalez C. A., et al., 2017, AJ, 154, 7Gonzalez J. F., Laibe G., Maddison S. T., Pinte C., Ménard F., 2015, MNRAS,454, L36Gregorio-Hetem J., Hetem A., 2002, MNRAS, 336, 197Guerri G., et al., 2011, Experimental Astronomy, 30, 59Haffert S. Y., Bohn A. J., de Boer J., Snellen I. A. G., Brinchmann J., GirardJ. H., Keller C. U., Bacon R., 2019, Nature Astronomy, 3, 749Hashimoto J., et al., 2012, ApJ, 758, L19Hashimoto J., et al., 2015, ApJ, 799, 43Hashimoto J., Aoyama Y., Konishi M., Uyama T., Takasao S., Ikoma M.,Tanigawa T., 2020, AJ, 159, 222Isella A., et al., 2018, ApJ, 869, L49Isella A., Benisty M., Teague R., Bae J., Keppler M., Facchini S., Pérez L.,2019, ApJ, 879, L25Keppler M., et al., 2018, A&A, 617, A44Keppler M., et al., 2019, A&A, 625, A118Kuhn J. R., Potter D., Parise B., 2001, ApJ, 553, L189Laibe G., Price D. J., 2014, MNRAS, 440, 2136Lodato G., Price D. J., 2010, MNRAS, 405, 1212Lodato G., et al., 2019, MNRAS, 486, 453Long Z. C., et al., 2018, ApJ, 858, 112Lubow S. H., Seibert M., Artymowicz P., 1999, The Astrophysical Journal,526, 1001Marois C., Lafrenière D., Macintosh B., Doyon R., 2006, ApJ, 647, 612Masset F., Snellgrove M., 2001, MNRAS, 320, L55Mentiplay D., Price D. J., Pinte C., 2019, MNRAS, 484, L130Mesa D., et al., 2019, A&A, 632, A25Milli J., Mouillet D., Lagrange A. M., Boccaletti A., Mawet D., Chauvin G.,Bonnefoy M., 2012, A&A, 545, A111Monaghan J. J., 2005, Reports on Progress in Physics, 68, 1703Müller A., et al., 2018, A&A, 617, L2Murakawa K., 2010, A&A, 518, A63Muto T., et al., 2012, ApJ, 748, L22Pecaut M. J., Mamajek E. E., 2016, MNRAS, 461, 794Pérez S., Casassus S., Baruteau C., Dong R., Hales A., Cieza L., 2019, AJ,158, 15Pinilla P., Birnstiel T., Ricci L., Dullemond C. P., Uribe A. L., Testi L., NattaA., 2012, A&A, 538, A114Pinte C., Ménard F., Duchêne G., Bastien P., 2006, A&A, 459, 797Pinte C., Harries T. J., Min M., Watson A. M., Dullemond C. P., Woitke P.,Ménard F., Durán-Rojas M. C., 2009, A&A, 498, 967Pinte C., et al., 2018a, A&A, 609, A47Pinte C., et al., 2018b, ApJ, 860, L13Pinte C., et al., 2019, Nature Astronomy, 3, 1109Price D. J., 2012, Journal of Computational Physics, 231, 759Price D. J., Laibe G., 2015, MNRAS, 451, 813Price D. J., et al., 2018, Publ. Astron. Soc. Australia, 35, e031Ragusa E., Dipierro G., Lodato G., Laibe G., Price D. J., 2017, MNRAS,464, 1449Ricci L., et al., 2017, The Astronomical Journal, 154, 24Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33MNRAS , 1–14 (2020) C. Toci al. -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6
Separation (") -0.6-0.4-0.20.00.20.40.6 S e p a r a t i o n ( " ) Figure A1.
Differential image obtained by subtracting the predicted scatteredlight image at 2.11/2.25 𝜇 m when the location of planet c is 180 degreesfrom the observed one to the predicted image where the planets are at theirobserved position. This highlights the faint spiral accretion stream inducedby planet c.Siess L., Dufour E., Forestini M., 2000, A&A, 358, 593Suriano S. S., Li Z.-Y., Krasnopolsky R., Shang H., 2017, MNRAS, 468,3850Szulágyi J., Mayer L., Quinn T., 2017, MNRAS, 464, 3158Szulágyi J., Plas G. v. d., Meyer M. R., Pohl A., Quanz S. P., Mayer L.,Daemgen S., Tamburello V., 2018, MNRAS, 473, 3573Teague R., Bae J., Bergin E. A., Birnstiel T., Foreman-Mackey D., 2018, ApJ,860, L12Testi L., et al., 2015, ApJ, 812, L38Toci C., Lodato G., Fedele D., Testi L., Pinte C., 2020, ApJ, 888, L4Uyama T., et al., 2020, AJ, 159, 118Veronesi B., Lodato G., Dipierro G., Ragusa E., Hall C., Price D. J., 2019,MNRAS, 489, 3758Woitke P., Kamp I., Thi W. F., 2009, A&A, 501, 383Woitke P., et al., 2016, A&A, 586, A103Youdin A. N., Goodman J., 2005, ApJ, 620, 459Zhang S., et al., 2018, ApJ, 869, L47Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ, 813, 88van der Plas G., et al., 2019, A&A, 624, A33 APPENDIX A: SPIRAL ACCRETION STREAM ININFRARED IMAGES
Given both the biases induced by angular differential imaging (ADI;Figure 3, third column) and the presence of small dust grains lead-ing to bright intra-cavity signals regardless of planet c’s position(Figure 4), it is unclear which scattered light signals result from theinteraction of planet c with the disc. In order to highlight this specificsignature we show in Figure A1 the image obtained when subtractingthe scattered light image with planet c located 180 deg away fromits observed position (and no planet flux injected) to the scatteredlight image obtained when the planets are injected at their observedposition angle and contrast. These images are obtained considering convolution with the observed point spread function and the effect ofa coronagraph, but without geometric bias induced by ADI. The faintspiral accretion stream associated to c can be seen both outward andinward from its location. In particular, the blob that is found inwardsuggests that part of the signal from Feature ii (Figure 3) could tracethe inner spiral arm of PDS 70 c. It is worth noting that the spiral armseen along the edge of the cavity at the bottom part of the image isnot a secondary arm but results from the subtraction of the primaryspiral associated to planet c when it is placed at 180deg difference intrue anomaly compared to observations.Spiral arms have been observed in scattered light observationsin a number of protoplanetary discs (e.g. Muto et al. 2012; Benistyet al. 2015; Uyama et al. 2020), and have been interpreted as possibleplanet signposts (e.g. Dong et al. 2015; Zhu et al. 2015). Howeverthey have not been reported in this system, despite the presenceof confirmed embedded protoplanets. The test presented in this ap-pendix suggests that the lack of reported spirals is most likely dueto the unfavourable viewing geometry of the disc, lack of angularresolution in the deprojected radial direction and faint additional fluxcompared to the bright forward-scattered edge of the cavity. Our con-clusion meets the one of Zhu et al. (2015) who showed that a highinclination can significantly hinder the visibility of spiral arms, inparticular in polarised light due to the phase function of the degreeof polarisation. APPENDIX B: DIFFERENTIAL IMAGING EFFECTS
The primary focus of this work was to study the dynamics of thePDS 70 protoplanetary system. However, the close match betweenour radiative transfer predictions and reported observations in theinfrared allows us to examine the biases in reported images of PDS70 caused by observational and post-processing techniques.Using a similar method to our post-processed K-band total inten-sity scattered light images (Section 2.5), we also predicted the maptracing the azimuthal component of polarised intensity ( 𝑄 𝜙 ) in 𝐽 band (1.25 𝜇 m) for comparison to the coronagraphic IRDIS obser-vations presented in Keppler et al. (2018). We used the Q and UStokes polarisation maps returned by mcfost, and applied steps ii to v detailed in Section 2.5: convolution with the observed PSF ( ∼ 𝑄 𝜙 map, as is now standard for polarimetric dif-ferential imaging observations (PDI; e.g. Kuhn et al. 2001; Avenhauset al. 2017): 𝑄 𝜙 = 𝑄 cos ( 𝜙 ) + 𝑈 sin ( 𝜙 ) , (B1)with 𝜙 = arctan 𝑥 − 𝑥 𝑦 − 𝑦 . (B2)Figure B1 (left column) shows the predicted 𝑄 𝜙 maps for sim-ulations 1 to 3. The size of the IR cavity is well reproduced insimulations 2 and 3. However we note that (1) the peak signal of thecavity edge lies along the PA of semi-major axis of the disc instead ofa PA ∼ MNRAS000
The primary focus of this work was to study the dynamics of thePDS 70 protoplanetary system. However, the close match betweenour radiative transfer predictions and reported observations in theinfrared allows us to examine the biases in reported images of PDS70 caused by observational and post-processing techniques.Using a similar method to our post-processed K-band total inten-sity scattered light images (Section 2.5), we also predicted the maptracing the azimuthal component of polarised intensity ( 𝑄 𝜙 ) in 𝐽 band (1.25 𝜇 m) for comparison to the coronagraphic IRDIS obser-vations presented in Keppler et al. (2018). We used the Q and UStokes polarisation maps returned by mcfost, and applied steps ii to v detailed in Section 2.5: convolution with the observed PSF ( ∼ 𝑄 𝜙 map, as is now standard for polarimetric dif-ferential imaging observations (PDI; e.g. Kuhn et al. 2001; Avenhauset al. 2017): 𝑄 𝜙 = 𝑄 cos ( 𝜙 ) + 𝑈 sin ( 𝜙 ) , (B1)with 𝜙 = arctan 𝑥 − 𝑥 𝑦 − 𝑦 . (B2)Figure B1 (left column) shows the predicted 𝑄 𝜙 maps for sim-ulations 1 to 3. The size of the IR cavity is well reproduced insimulations 2 and 3. However we note that (1) the peak signal of thecavity edge lies along the PA of semi-major axis of the disc instead ofa PA ∼ MNRAS000 , 1–14 (2020) ynamical model of PDS 70 b c iiiiii iv Figure B1.
Near infrared observations of PDS 70 (top row) compared to the processed IR predictions obtained from simulations 1, 2 and 3 (second, third andlast row, respectively).
Left column:
Polarimetric Differential Imaging predictions in 𝐽 band (1.25 𝜇 m). Middle column: average median-ADI images over theK1 and K2 bands (2.11/2.25 𝜇 m). Right column: average PCA-ADI images over the H2 and H3 bands (1.59/1.67 𝜇 m), obtained with 5 principal componentssubtracted. MNRAS , 1–14 (2020) C. Toci al.
Another minor difference between prediction and observation is thepresence of an intra-cavity ring-like emission that is not observedin PDI. We hypothesise that a better match would be obtained if thesimulations were run for more orbits of the companions as this wouldfurther deplete the cavity of small grains. In a similar way, this wouldlikely reduce the intensity of feature ii in the median-ADI images,and make it more in line with the fainter signal (with respect to theplanets) seen in the observation.We have also processed our infrared images using the smart prin-cipal component analysis (sPCA) algorithm implemented in vip(Amara & Quanz 2012; Absil et al. 2013; Gomez Gonzalez et al.2017), similar to Keppler et al. (2018). This is shown in the right col-umn of Figure B1. We followed the exact same steps as described inSection 2.5, apart from step viii where we used sPCA-ADI instead ofmedian-ADI. Compared to the median-ADI images (centre column),the bulge in the middle of the forward-scattered edge of the cavity(feature i ) is attenuated, but several parallel arcs appear along the SWand NW parts of the edge. The number and intensity of these arcs in-crease with increasing number of principal components. We find theclosest match to the observed image for 4 to 10 principal componentssubtracted. While increasing the number of principal components in-duces more geometric biases to extended disc signal, it nonethelessappears to help isolating the signal from planet c lying along theedge. Indeed, the blob gets flanked by characteristic negative ADIside lobes in azimuth, of increasing intensity with increasing numberof principal components subtracted. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000