PDS 70: A transition disk sculpted by a single planet
The wide, deep cavities of transition disks are often believed to have been hollowed out by nascent planetary systems. PDS 70, a Myr old transition disk system in which a multi-Jupiter-mass planet candidate at 22 au coexists with a au gas and au dust gap, provides a valuable case study for this hypothesis. Using the PEnGUIn hydrodynamics code, we simulate the orbital evolution and accretion of PDS 70b in its natal disk over the lifetime of the system. When the accreting planet reaches about 2.5 Jupiter masses, it spontaneously grows in eccentricity and consumes material from a wide swathe of the PDS 70 disk; radiative transfer post-processing with DALI shows that this accurately reproduces the observed gap profile. Our results demonstrate that super-Jupiter planets can single-handedly carve out transition disk cavities, and indicate that the high eccentricities measured for such giants may be a natural consequence of disk-planet interaction.
A planetary-mass companion, PDS 70b, has recently been directly imaged by SPHERE around the pre-main sequence star PDS 70 (Müller et al., 2018). It is located about 22 au away from its host star (Keppler et al., 2018; Wagner et al., 2018), and is inside a transition disk cavity that stretches to 60 au (Hashimoto et al., 2012; Dong et al., 2012; Hashimoto et al., 2015; Müller et al., 2018). Is PDS 70b the cause of this cavity? If so, what can we tell about its properties and evolutionary history from its interaction with the disk?
Transition disk cavities have long been suspected to be emptied by gravitational disk-planet interaction. The PDS 70 system presents a valuable testing ground for this hypothesis. Planets repel material from their orbits by exerting Lindblad torques, and so open gaps in their disks. Past studies have implied that these gaps are narrow, with widths only a fraction of the orbital radius of the gap-opening planet (e.g., Crida et al., 2006; Fung et al., 2014; Duffell & Chiang, 2015; Kanagawa et al., 2016; Ginzburg & Sari, 2018). Transition disk cavities, on the other hand, are much wider, often spanning tens of au. At first glance, PDS 70b’s location at 22 au appears too close in to explain a 60 au cavity.
In reality, however, gaps, even those opened by a single planet, can be much wider than previously expected. There are three reasons for this. First, continuum observations in sub-mm (e.g., Hashimoto et al., 2015; Long et al., 2018) probe only the large dust grain distribution. Dust tends to concentrate at local pressure maxima, leaving regions interior empty, whereas the gas density can show a more gradual decline (van der Marel et al., 2016; Dong et al., 2017). This is indeed the case in PDS 70, where Long et al. (2018) have found that the gas cavity extends to only 33 au. The inner radius for small dust grains (which tend to follow the gas distribution) was found to be well inside the millimeter dust grain cavity as well (Hashimoto et al., 2015).
Second, planets can accrete and remove gas directly from their disks. Evolutionary models suggest that PDS 70b is between 2 to 17 Jupiter masses () (Müller et al., 2018), while disk modeling places the disk’s gas mass at about 7 (Long et al., 2018). Furthermore, H observations with MagAO (Wagner et al., 2018) indicate that the planet is still growing. It is therefore likely that over its lifetime, PDS 70b has consumed a significant fraction of the disk, which would conceivably deepen and widen the gap.
Third, planetary orbits can be eccentric. Radial velocity surveys find that giant planets have eccentricities distributed nearly uniformly between 0 and 0.8 (Butler et al., 2006). When the planet opens a deep gap, eccentric Lindblad resonances are known to be able to excite eccentricity (Artymowicz, 1993; Goldreich & Sari, 2003; Sari & Goldreich, 2004; Duffell & Chiang, 2015). The 1:3 resonance is particularly effective when the gap is much wider than the local scale height (Papaloizou et al., 2001). Numerical simulations have verified that the eccentricity of the disk/planet can grow to for multi-Jupiter-mass planets (Papaloizou et al., 2001; Kley & Dirksen, 2006; Bitsch et al., 2013), which should result in significantly wider gaps.
Taking into account all of these factors, we produce a self-consistent hydrodynamics simulation that tracks the accretion and migration of PDS 70b over 4.6 Myr ( orbits). Our fiducial model for the disk-planet system successfully reproduces the location and mass of the planet, along with the cavity size of the disk. §2 describes our numerical methods. §3 presents our model, §4 shows how it compares with ALMA observations, and §5 concludes and discusses future directions.
2 Hydrodynamics Simulations
We use the graphics processing unit (GPU)-accelerated code PEnGUIn (Fung, 2015) to simulate disk-planet interaction in 2D. With the exception of planetary accretion, this code is identical to that used by Fung & Lee (2018) who similarly simulated planetary migration in disks. We describe here our accretion prescription, and refer the reader to Fung & Lee (2018) for a description of the other parts of the code.
The planet accretes mass from each of the nearby cells at a rate of:
|where is the surface density of gas, the gravitational constant, the mass of the planet, and the typical timescale for material at to free-fall onto the planet. We choose , where is the planet’s Hill radius, the star’s mass, and the instantaneous distance between the planet and star. We limit this rate by a Gaussian function around the planet:|
where is the maximum accretion rate. We set it to Myr, although in practice the rate is always lower. The actual planetary accretion rate may then be computed as
integrated over the full domain.
While this prescription does not explicitly conserve angular momentum, the small means that the material consumed co-orbits with the planet and has nearly the same specific angular momentum. We note that the planet’s potential is softened with a smoothing length of times the local disk scale height, which is nearly always larger than .
2.1 Initial and Boundary Conditions
with , or approximately 0.05 times the surface density of the minimum-mass solar nebula (MMSN) at 1 au (Hayashi, 1981). Based on calculations with the radiative transfer code DALI (Bruderer, 2013, and Section 4), we impose a fixed sound-speed profile of
with . This results in a flaring disk with aspect ratio , and a temperature profile . At 22 au, this yields a temperature of 35 K for a disk with a mean molecular weight of 2.34. For disk viscosity, we choose the Shakura-Sunyaev viscosity parameter to be (Shakura & Sunyaev, 1973). This choice is motivated by the fact that PDS 70 has no clear large-scale asymmetry, implying that vortex formation at gap edges through the Rossby wave instability (e.g., Lovelace et al., 1999; Li et al., 2000; Lovelace & Romanova, 2014) is suppressed. Our numerical tests show this requires roughly .
Our initial condition enforces hydrostatic equilibrium, where the radial velocity is zero everywhere, and the azimuthal velocity is
where is the Keplerian orbital frequency. Per Keppler et al. (2018), we set .
Our simulation domain spans 5 to 100 au in radius and the full 2 in azimuth. Radial boundaries are fixed to their initial values. To ensure that our fields remain remain smooth there, we implement wave-killing zones similarly to the prescription of Fung et al. (2015), damping all fields within 2 local scale heights of the boundaries to their initial values over 10 local orbital periods.
In our fiducial run, the planet is initialized as a 10 core with an initial semimajor axis of 23 au; over the first 100 Kyr of the simulation it migrates to 22 au, where it remains over the run’s lifetime. We compare the results to three “control” runs, in which non-accreting planets of 2, 4, and 8 are fixed on circular, 22 au orbits.
The grid dimensions used for the fiducial and control simulations are 960 () 1944 () cells, spaced logarithmically in the radial and uniformly in azimuth; this corresponds to 20 cells per scale height at 22 au. In tests at both 50% higher and lower than our fiducial resolution, as well as at different viscosities, we find that the position of the planet agrees to within 8%.
We test over 100 different models, varying parameters such as the planet’s starting position ( au), disk surface density (), and viscosity (). We also enable and disable mechanisms such as planetary migration and accretion. Our investigations lead to a best-fit fiducial model, in which the planet is free to both migrate and accrete, with and . To accommodate the slow growth of the planet, we run it to 4.6 Myr ( orbits), approximately the current age of PDS 70 (Müller et al., 2018). In figures 1 and 2, we compare the fiducial run to three representative control simulations in which planet masses and orbits are fixed; these are run to 250 Kyr ( orbits), by which time they should have reached steady-state (e.g., Fung & Chiang, 2016).
By the end of the fiducial simulation, the planet has depleted the gas below out to 30 au, and produced a pressure peak at 62 au. Observationally, these would manifest as a 30 au optically thin cavity in CO and a au ring in dust continuum, which we verify with radiative transfer modeling in Section 4. Accretion and eccentricity together enable the fiducial planet to assimilate a substantial fraction of the disk and produce a wide, deep cavity that resembles PDS 70, which would not be possible unless both mechanisms are accounted for.
3.1 Disk-Planet Evolution
In agreement with the results of Kley & Dirksen (2006) and Fung et al. (2014), we find in all of our simulations that the disk becomes eccentric when the planet-to-star mass ratio ( for PDS 70). This is seen in the 4 and 8 control runs, where the disk remains visibly eccentric to the end (Figure 1). It also occurs in the fiducial run when the planet grows to , but the subsequent evolution differs because the planet experiences a back-reaction from the disk.
At the start of our fiducial simulation, the planet rapidly accretes material at a rate of , which decreases as the gap becomes wider and deeper. By 3 Myr, when the planet reaches a mass of , the disk’s surface density peak is located at 48 au—% farther than in the comparable control simulation—and the gap is depleted by a factor of . Throughout this period, the planet’s eccentricity , and its semimajor axis is nearly constant at 21.7 au.
Past 3 Myr, the planet grows to . The gap is now wide enough for eccentric Lindblad resonances to strongly influence disk-planet interaction, making the disk eccentric (e.g., Papaloizou et al., 2001). In response, the planet’s eccentricity grows exponentially, before saturating at . This has two major consequences—first, it distributes the planet’s torque over a larger radial range, making the gap substantially wider, albeit modestly shallower. Second, it enables the planet to access and consume more disk material, further enlarging the gap; see Figure 2 for a comparison. Outside the now au-wide cavity, the disk is now far enough from the planet’s gravity to remain essentially circular.
At very late times, we observe renewed growth in the planet’s eccentricity. Unlike a real disk, whose mass (and thus gravitational influence on the planet) would be photoevaporatively depleted on multi-Myr timescales, our simulated disk can never fully drain, meaning that the planet may grow in mass and eccentricity indefinitely. We therefore consider 4.6 Myr (approximate age of PDS 70) a reasonable point to terminate our fiducial run.
3.2 Accretion and Disk Mass
Disk evolution in our fiducial model is qualitatively consistent over a wide range of disk masses, which we demonstrate using two runs with identical initial conditions, except with 2 and 3 the surface density and a lower resolution, 634 () 1296 (). These low resolution runs reproduce the timescale for eccentricity excitation to within 10% of that at the fiducial resolution. As shown in Figure 3, planetary mass and eccentricity evolve similarly to our fiducial run, but the increased disk mass means that the planet can more rapidly reach the critical for eccentricity excitation.
Over Myr timescales, a lower disk mass may help stabilize the system. When the initial disk mass is of order a few , as in our fiducial run, the formation of a giant planet would remove enough mass to severely inhibit migration via disk-planet interaction. If the initial disk mass is too high, or the planet accretion rate is too low, migration could move the planet across the entire disk in very short timescales. Indeed, we have observed such fast migration in some of our tests with non-accreting, super-Jupiter planets.
4 Radiative transfer modeling
Using DALI, we compute the expected CO 3–2 and 350 GHz dust continuum emission of our fiducial model, and compare with ALMA observations. DALI is a physical-chemical modeling code (Bruderer, 2013) which self-consistently solves for gas temperature, molecular abundances, and excitation given gas and dust surface density profiles. For the gas surface density, we use the output of the fiducial run. For dust, we assume a power law with an exponential tail:
with , and , with the surface density reduced by 1000 for . Based on the pressure peak found in our fiducial simulation, we set (see Figure 4). We extrapolate the dust profile for and using equation (5), and the gas profile by a piecewise version, choosing coefficients in the inner and outer disk that ensure continuity. This procedure gives an overall gas-to-dust ratio of 24. For vertical structure, we assume an aspect ratio with and , consistent with the infrared excess in the Spectral Energy Distribution (SED) as visible in Figure 4. Images are raytraced at 113 pc (Gaia Collaboration et al., 2018). More details on this procedure are given in van der Marel et al. (2016).
We have reduced the ALMA Band 7 observations of Long et al. (2018) of the 350 GHz continuum and CO 3–2 emission using the default pipeline scripts, and cleaned the data using natural weighting, resulting in a beam size of 0.22”0.18”. Figure 4 shows that the SED and continuum ring closely match the data. The CO map shows a similar morphology, indicating that the CO emission is sufficiently optically thin at the planet’s location. The total integrated flux of the model is also comparable. However, our model predicts a bright central peak in CO, significantly stronger than that observed—likely the result of our inner-disk boundary conditions. We address this point more thoroughly in the following section (§5).
5 Conclusions and Discussions
Using 2D hydrodynamics simulations and radiative transfer post-processing, we demonstrate that the super-Jupiter PDS 70b, located 22 au from its host star, is by itself capable of carving out the au gas cavity and au dust cavity in which it resides (Figures 1 and 2). Radiative transfer modeling yields results consistent with the observed continuum emission and SED of PDS 70 (Figure 4).
Our model makes the following predictions about PDS 70:
PDS 70b is on an eccentric orbit with . This naturally results from disk-planet interaction and is required to explain the width of the observed cavity.
The planet formed in situ; the low surface density required to reproduce the observed gas cavity could not have caused significant migration in a forming super-Jupiter planet.
The disk and planet mass are comparable, implying that the planet must have consumed a large region of the disk as it evolved.
One alternative would be for PDS 70b to occupy a circular orbit, with additional, unseen planets at large radii carving out the observed cavity. Given PDS 70b’s mass, we find this scenario unlikely—suppressing eccentricity excitation would entail a disk density an order of magnitude or more lower than that observed. If additional planets exist in the system, they are likely much less massive than PDS 70b and have little impact on the disk morphology.
Our fiducial model has several areas for improvement. For instance, the optically thick CO emission it predicts from the inner disk is far in excess of that observed, implying that its simulated size—approximately 10 au—is likely too large. This size is influenced by the location of the inner hydrodynamical boundary, which we have set to 5 au because this work focuses on large-scale ring and gap, rather than the inner disk. Additionally, the details of gas-grain interaction and vertical stratification in the inner disk could significantly elevate the temperature and emission of CO in our fiducial model as compared to that in the real PDS 70. Simulating a closer inner boundary and testing alternative vertical structures would improve our fit of the inner disk, but at a substantially elevated computational cost.
With our current planetary-accretion prescription, the planet simply consumes the material in its vicinity, subject to an artificial maximum rate. Simulating thermal physics and circumplanetary disk dynamics would better constrain the mass and growth rate of the planet, and thus the long-term evolution of the disk. Such changes, however, would not affect our central conclusion—that a growing super-Jupiter planet would spontaneously become eccentric, and consume a large radial range of the disk to create a wide, deep cavity.
Another concern is that our simulations are in 2D. Running our 4.6 Myr long fiducial model in 3D would currently be impractical, requiring years in wall-clock time. 2D may be sufficient — previous work has shown that 3D gap-opening (Fung & Chiang, 2017) and planetary torques (Fung et al., 2017) are reasonably well-captured in 2D; however, it remains to be seen how closely 2D eccentricity excitation and planetary accretion would match their 3D counterparts. 2D also cannot address whether the orbit of PDS 70b is inclined, and the effects this would have on the cavity structure. In the future, detailed studies of the 3D dynamics would improve understanding of how gas giants and transition disks form.
Our findings also raise interest in the orbital evolution of giant planets themselves—in particular, why super-Jupiter planets are often observed to have significant eccentricity (Butler et al., 2006). Our simulations show that disk-planet interaction may be a cause, and we intend to investigate this more thoroughly in a forthcoming paper.
- Artymowicz (1993) Artymowicz, P. 1993, ApJ, 419, 166
- Bitsch et al. (2013) Bitsch, B., Crida, A., Libert, A.-S., & Lega, E. 2013, A&A, 555, A124
- Bruderer (2013) Bruderer, S. 2013, A&A, 559, A46
- Butler et al. (2006) Butler, R. P., Wright, J. T., Marcy, G. W., et al. 2006, ApJ, 646, 505
- Crida et al. (2006) Crida, A., Morbidelli, A., & Masset, F. 2006, Icar, 181, 587
- Dong et al. (2012) Dong, R., Hashimoto, J., Rafikov, R., et al. 2012, ApJ, 760, 111
- Dong et al. (2017) Dong, R., van der Marel, N., Hashimoto, J., et al. 2017, ApJ, 836, 201
- Duffell & Chiang (2015) Duffell, P. C., & Chiang, E. 2015, ApJ, 812, 94
- Fung (2015) Fung, J. 2015, PhD thesis, University of Toronto, Canada
- Fung et al. (2015) Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101
- Fung & Chiang (2016) Fung, J., & Chiang, E. 2016, ApJ, 832, 105
- Fung & Chiang (2017) —. 2017, ApJ, 839, 100
- Fung & Lee (2018) Fung, J., & Lee, E. J. 2018, ApJ, 859, 126
- Fung et al. (2017) Fung, J., Masset, F., Lega, E., & Velasco, D. 2017, AJ, 153, 124
- Fung et al. (2014) Fung, J., Shi, J.-M., & Chiang, E. 2014, ApJ, 782, 88
- Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
- Ginzburg & Sari (2018) Ginzburg, S., & Sari, R. 2018, MNRAS, 479, 1986
- Goldreich & Sari (2003) Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024
- Hashimoto et al. (2012) Hashimoto, J., Dong, R., Kudo, T., et al. 2012, ApJL, 758, L19
- Hashimoto et al. (2015) Hashimoto, J., Tsukagoshi, T., Brown, J. M., et al. 2015, ApJ, 799, 43
- Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
- Kanagawa et al. (2016) Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, PASJ, 68, 43
- Keppler et al. (2018) Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44
- Kley & Dirksen (2006) Kley, W., & Dirksen, G. 2006, A&A, 447, 369
- Li et al. (2000) Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023
- Long et al. (2018) Long, Z. C., Akiyama, E., Sitko, M., et al. 2018, ApJ, 858, 112
- Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805
- Lovelace & Romanova (2014) Lovelace, R. V. E., & Romanova, M. M. 2014, Fluid Dynamics Research, 46, 041401
- Müller et al. (2018) Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2
- Papaloizou et al. (2001) Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263
- Sari & Goldreich (2004) Sari, R., & Goldreich, P. 2004, ApJL, 606, L77
- Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, in IAU Symposium, Vol. 55, X- and Gamma-Ray Astronomy, ed. H. Bradt & R. Giacconi, 155
- van der Marel et al. (2016) van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58
- Wagner et al. (2018) Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJL, 863, L8