TW Hya: an old protoplanetary disc revived by its planet

Dark rings with bright rims are the indirect signposts of planets embedded in protoplanetary discs. In a recent first, an azimuthally elongated AU-scale blob, possibly a planet, was resolved with ALMA in TW Hya. The blob is at the edge of a cliff-like rollover in the dust disc rather than inside a dark ring. Here we build time-dependent models of TW Hya disc. We find that the classical paradigm cannot account for the morphology of the disc and the blob. We propose that ALMA-discovered blob hides a Neptune mass planet losing gas and dust. We show that radial drift of mm-sized dust particles naturally explains why the blob is located on the edge of the dust disc. Dust particles leaving the planet perform a characteristic U-turn relative to it, producing an azimuthally elongated blob-like emission feature. This scenario also explains why a 10 Myr old disc is so bright in dust continuum. Two scenarios for the dust-losing planet are presented. In the first, a dusty pre-runaway gas envelope of about 40 Earth mass Core Accretion planet is disrupted, e.g., as a result of a catastrophic encounter. In the second, a massive dusty pre-collapse gas giant planet formed by Gravitational Instability is disrupted by the energy released in its massive core. Future modelling may discriminate between these scenarios and allow us to study planet formation in an entirely new way -- by analysing the flows of dust and gas recently belonging to planets, informing us about the structure of pre-disruption planetary envelopes.


INTRODUCTION
Recent advent in the high resolution imaging of protoplanetary discs via scattering light techniques and mmcontinuum with ALMA yielded many exciting examples of sub-structures in the discs, such as large scale assymetries (Casassus et al. 2013), spirals (Benisty et al. 2017), dark and bright rings and/or gaps (ALMA Partnership et al. 2015;Long et al. 2018;Andrews et al. 2018;Dullemond et al. 2018), clumps and young planets (Mesa et al. 2019). In fact it is believed that most protoplanetary discs have sergei.nayakshin@le.ac.uk substructures; those that do not may simply be those that have not yet been imaged at high enough resolution (Garufi et al. 2018). The fact that these substructures are seen in very young discs, and that the masses of these discs appear to be insufficient to form planetary systems that we observe (Greaves & Rice 2010;Manara et al. 2018;Williams et al. 2019), indicates that planets form very rapidly, possibly faster than 1 Myr.
On the one hand, these observations bring new challenges. The presence of mm-sized dust in a few Myr old discs is surprising because of the rapid inward radial drift of the grains (Weidenschilling 1977;Birnstiel et al. 2012). While the radial drift can be slowed down by invoking very mas-sive gas discs, M gas (0.1 − 0.2)M (e.g., Powell et al. 2017Powell et al. , 2019, other arguments point against that. Veronesi et al. (2019) shows that the prevalence of annular rather than spiral features in many of the observed discs indicates that the mm-sized particles have rather large Stokes numbers, limiting disc masses to only ∼ 1 M J . Furthermore, hydrodynamical simulations and population synthesis show that planets with properties deduced from these observations should both grow in mass and migrate inward very rapidly (Clarke et al. 2018;Mentiplay et al. 2019;Lodato et al. 2019;Nayakshin et al. 2019;Ndugu et al. 2019). This would make the detection of these planets statistically very unlikely; this paradox is resolved if the disc masses are ∼ 1 M J so that planets do not migrate rapidly (Nayakshin 2020).
Here we focus on the protoplanetary disc in TW Hydra (Kastner et al. 1997). At the distance of about 60 pc, this protostar is the closest one with a protoplanetary disc, and is probably the best studied. Despite its advanced age of t * ≈ 10 Myr (Weinberger et al. 2013), TW Hydra continues to accrete gas at a respectable rate (∼ 2 × 10 −9 M yr −1 ). Its disc is the prototype for discs with cavities possibly carved by growing planets, with early observations indicating gaps on the sub-AU to a few AU scales (Calvet et al. 2002;Eisner et al. 2006). Recent high resolution ALMA observations found several axisymmetric gaps in the continuum sub-mm dust emission of TW Hya, one on the scale of ∼ 1 AU, and two more at ∼ 24 and 41 AU Huang et al. 2016;Tsukagoshi et al. 2016;Huang et al. 2018).
However, the currently unique feature of TW Hydra is the first ever image of a potential planet in the ALMA 1.3 mm dust continuum emission (Tsukagoshi et al. 2019, T19 hereafter). The few AU-scale emission excess is significant (12σ over the disc background intensity) and has an azimuthally elongated shape. T19 associated it with a circumplanetary disc of a growing Neptune mass planet. Most curiously the putative planet is located not inside the previously discovered gaps but at 51.5 AU from the star, right on the edge of a well known cliff-like rollover in the dust disc (Andrews et al. 2012Hogerheijde et al. 2016). In contrast, emission of many molecular tracers and from microscopic dust have broad peaks at 50 − 60 AU, and extending to distances as large as 100 − 250 AU (Kastner et al. 2015;Bergin et al. 2016;Teague et al. 2018).
Below we model the time evolution of the gas and dust components of TW Hydra's disc and compare the resulting dust continuum emission with ALMA bands 7 and 6 (820µm and 1.3 mm, respectively) and EVLA 9 mm observations of the source. We investigate several possible scenarios to try and explain both the disc morphology, with its sharp rollover of the dust disc, and the image of a putative planet positioned at the edge of that rollover. We find multiple lines of argument that challenge the standard quasi steady-state scenario for this disc. We propose that this disc has been recently re-invigorated by mass injection from a disrupted planet.
The paper is structured as follows. In §2 we summarise the observations of TW Hydra relevant to this paper. In §3 we give simple analytical arguments that challenge the quasi steady-state scenario. In §4 we present our numerical methods. In §5 we apply these to the quasi steady state disc scenario, finding significant and additional to §3 problems with it. We test a phenomenological Dust Source model in §6, in which a low-mass object on a fixed circular orbit ejects dust into the surrounding, initially dust-free, disc. This scenario proves promising but is not physically self-consistent. In §7 we propose that a collision between a ∼ (30 − 40) M ⊕ pre-collapse Core Accretion planet with a dust-rich envelope and another massive core could unbind the envelope and make the planet a dust source. In §8 we show that a gas giant planet formed in the Gravitational Instability scenario could be disrupted from within by its massive luminous core provided that the planet is very metal rich (Z ∼ 0.1) and its mass is M p 2 M J .
Since the planet is formed in this scenario very early on (presumably at t 0.1 Myr), survival of the planet on a very wide orbit by t ∼ 10 Myr requires that TW Hydra had neither gas nor dust disc before the disruption of the planet. In §9 we tackle the issue of the "blob" morphology in the Tsukagoshi et al. (2019) ALMA image, e.g., its spatial extent and its elongation along the azimuthal direction. We conclude with an extended discussion in §10.

DUSTY PUZZLES OF TW HYDRA
TW Hydra is an oddity in many regards. The mass of the mm-sized dust in TW Hya is surprisingly large, estimated at ∼ 80 M ⊕ (after correcting to the new GAIA distance; Andrews et al. 2016). Williams et al. (2019) shows that the mean dust mass of class I sources (generally thought to be much younger than 1 Myr) is ∼ 4 M ⊕ , whereas the mean mass of class II sources is ∼ 0.8 M ⊕ . Thus, despite being older by a factor of a few than an average type II disc, TW Hya's dust disc is ∼ 100 times more massive in dust than the mean 1 .
The mass of H 2 gas disc in TW Hya is suspected to be very high, although there are significant uncertainties. Bergin et al. (2016) inferred TW Hya H 2 disc mass from HD line observations to be 0.05M but this value was revised to between 0.006 − 0.009M by Trapman et al. (2017) with an alternative model and more data. However, as shown by Jones et al. (2012), for a disc evolving viscously, its mass at age t * is given by M d ≈ ξ d M * t * ≈ 0.1M (ξ d /5), where ξ d a few. Here the accretion rate inferred for TW Hya is M * = 1.8 × 10 −9 M yr −1 (Ingleby et al. 2013). Independently of this argument, Powell et al. (2019) use the dust radial drift constraints to estimate the gas disc mass in TW Hya as ∼ 0.11M 2 . Such a high mass is ∼ 2 orders of magnitude higher than the median gas disc mass for type II discs (Manara et al. 2018). Fig. 1 shows the summary of TW Hydra observations that are relevant to this paper. The top left panel shows 1 We shall see later that TW Hydra's dust mass may be as small as ∼ 15 M ⊕ if DIANA (Woitke et al. 2016) opacities ( §4.4) are used. Williams et al. (2019) 1.3 mm opacity is quite similar to that used by Andrews et al. (2016). The higher DIANA opacities would reduce the dust mass estimates made by Williams et al. (2019)  the image of the source in the 1.3 mm continuum presented in Tsukagoshi et al. (2019). The two annular gaps are well known from the previous ALMA observations Tsukagoshi et al. 2016;Huang et al. 2018). The white box in the left panel of Fig. 1 is centred on the location of the excess. The half widths of the excess are ∼ 0.5 AU in the radial and ∼ 2.2 AU in the azimuthal directions. The bottom panel of fig. 1 shows with the green and red curves the ALMA azimuthally averaged dust continuum intensity profiles at 1.3 mm and 820µm, respectively, and the EVLA 9 mm intensity with the blue curve 3 . The radial profiles of the ALMA images were extracted by averaging the full azimuthal angle except for the position angles of the millimeter blob (P.A.=222 • − 236 • and 238 • − 252 • ). For the profile of the EVLA data, we took an average over the full azimuthal angle. The profiles were made after the image deprojection under the assumption that the disk inclination is 7 • and the position angle is −30 • from the North. We see that the ALMA continuum emission plunges by an order of magnitude from R ∼ 52 AU to R = 70 AU. This behaviour is not apparent in the EVLA data, but because the beam size is ∼ 15 AU it is possible that the intrinsic disc I ν also has a similarly sharp rollover at the same location.
While dust discs much smaller than the gas discs as traced by CO emission are not uncommon, and while there are other dust discs with sharp outer cutoffs (e.g., Trapman et al. 2019), TW Hya is certainly the best resolved disc of this kind. Hogerheijde et al. (2016) searched for a broken power-law fit to ALMA data at the wavelength of 820µm. They found the power-law index of p ≈ −0.5 before the break, and Σ dust ∝ R −8 beyond the break. We show this Σ dust (R) profile with the dot-dash violet curve in Fig. 1. Note that the location of the T19 planet coincides very well with the break in dust I ν and Σ dust . The gas surface density profile is much harder to constrain although it is clear that the gas disc extends at least to 200 AU. The three black curves in the bottom panel of fig. 1 show three gas surface density models from the literature. While the total gas disc mass varies by about an order of magnitude between these, it is clear that the gas distribution is not as compact as that of the dust.

ANALYTICAL PRELIMINARIES: THE STEADY-STATE SCENARIO
In the classical protoplanetary disc paradigm, the discs are formed during the collapse of the parent molecular cloud and persist until dispersed after ∼ (3 − 10) Myrs (Alexander et al. 2014). This Steady State scenario is explored numerically in §5, here we scrutinise it with more transparent analytical arguments.

Disc model
The accretion rate onto the star in the viscous steady state is (Shakura & Sunyaev 1973) where v K = Ω K R = (GM * /R) 1/2 is the local circular Keplerian speed, Σ is the local gas surface density, ν = αc s H is the disc viscosity, α is the viscosity parameter, H is the disc vertical scale height, and c s is the gas sound speed. We assume that the disc midplane temperature is given by where T 0 = 18 K is the gas temperature at the planet location, R 0 = 51.5 AU (this temperature profile is very similar to the midplane temperature profile derived in Huang et al. 2016). For reference, v K = 3.7 km/s, c s = 0.25 km/s, and H/R ≈ 0.07 at R = R 0 . In the whole of §3 we also assume that Equation 1 constrains the product of αΣ 0 . We shall use additional considerations to constrain Σ, which then limits α via eq. 1.

Gas accretion constraints on the disc mass
Numerical integrations of viscously evolving discs show that the accretion rate onto the star is large at early times and then decreases (Jones et al. 2012). The product of the accretion rate and current time t remains approximately constant initially and later also decreases (because M * drops with time). Usefully for us here, the disc mass at age t * is given by where ξ d ∼ a few is a dimensionless number (Jones et al. 2012), and we used M * = 1.8 × 10 −9 M yr −1 and t * = 10 Myr.

Dust drift constraints
As previous authors (Andrews et al. 2012) we find that particles of size a 1 mm are required to explain ALMA observations. The radial dust drift velocity is (Whipple 1972 where η = 1.25 given our model disc pressure profile, and the Stokes number is given by with ρ a = 2.1 g cm −3 (Woitke et al. 2019) being the grain material density. For a p = 3.5 power-law grain size distribution with maximum grain size a, the appropriate grain size to use in eqs. 5 & 6 in the small St 1 regime is the mean grain size, a mean = a/3 (see §4.3). Demanding the dust particle drift time scale to be 10 Myr at the location of T19 planet we arrive at the minimum gas surface density at ∼ 50 AU of With the observed gas accretion rate of M * ≈ 1.8 × 10 −9 M yr −1 , we simultaneously have that The disc mass enclosed within radius R is then This is larger than M d = 0.11M found by Powell et al. (2019), mainly because they assumed that TW Hydra is younger (5 Myr).
The mass in eq. 9 is uncomfortably large for many reasons. Firstly, the typical gas mass of a ∼ few Myr old disc is estimated at ∼ 10 −3 M (Manara et al. 2018). Secondly, such a massive disc should be strongly self-gravitating since the Toomre (1964) when evaluated at R = R 0 . We expect the disc to show spiral density structure, and in fact be fragmenting for such a low Q. Indeed, a 3D Phantom SPH calculation confirmed that the disc with the structure introduced at §3.1 and extending to 200 AU (as observed) fragments due to self-gravity. Another interesting inference from the dust drift constraints is the dust-to-gas ratio in the TW Hydra disc. According to Hogerheijde et al. (2016), Σ dust ≈ 0.2 at R ∼ 50 AU (although this depends on the dust opacity and size distribution). Hence, ε 0 = Σ dust /Σ ≈ 0.003. This is only slightly larger than that derived by Woitke et al. (2019) for TW Hydra, who obtained ε 0 ≈ 0.002. If the observed emission of the T19 feature is due to a circumplanetary disc around the planet, we would then expect the disc to have a dust-to-gas ratio smaller than ε 0 because dust in the circumplanetary disc is expected to drift into the planet faster than the gas component of the disc does (similarly to the protoplantary disc case). Therefore, we conclude that the minimum gas mass of the circumplanetary disc (CPD) is where M cpd−dust ≈ 0.03 M ⊕ is the dust mass of the T19 feature (Tsukagoshi et al. 2019). This CPD mass is surprisingly high and would require the planet itself to be much more massive, e.g., many tens of M ⊕ to maintain dynamical stability. This high M p in combination with high M d is ruled out due to planet migration ( §3.6) and spectral constraints ( §5).

Dust particles maximum size
Here we discuss three processes that limit grain growth in our numerical models below. We also use these results to argue that the quasi steady state scenario cannot naturally explain the relatively small grain size in this old disc (cf. further §10.1, item (vi)). Birnstiel et al. (2012) conclude that dust particle collisions due to radial drift are not likely to be a dominant mechanism of dust size regulation. For the steady state scenario of TW Hydra disc in particular, the drift velocity can be estimated by requiring the dust to not drift all the way into the star in 10 Myr: Grains of similar sizes collide at ∼ 1/2 of this velocity, e.g., at just about 1 cm s −1 . This velocity is two orders of magnitude smaller than the grain breaking velocities inferred from experiments and typically considered in the field, v br ∼ 10 m s −1 . Note however that this process can be dominant for planet-disruption scenarios where the relevant time scale is ∼ 10 5 yr. The grain growth time is finite. Birnstiel et al. (2012) introduce the "drift limit" to the maximum grain size as the grain size to which the grains grow before they are efficiently removed by the radial drift: at the location of the T19 planet (we used here Σ dust = 0.2 g/cm 2 from Hogerheijde et al. (2016)). This value is quite large, indicating that this process is also unlikely to stem grain growth in TW Hydra. Finally, gas turbulence also sets a maximum grain size (Weidenschilling 1984), which in the limit of small Stokes number yields Here we used the turbulent viscosity parameter α t , which may in general be different from the Shakura & Sunyaev (1973) α viscosity parameter introduced previously. The latter parameterizes the efficiency of the angular momentum transfer and subsumes in itself both gas turbulence and the effects of possible large scale magnetic torques (Bai & Stone 2013;Bai 2016). Since α ≥ α t , we estimate the minimum size set by turbulence via setting α t = α: where α −4 = 10 4 α (cf. equation 8). For the typical values of v br ∼ 10 m s −1 employed in the literature (e.g., Birnstiel et al. 2012), grains should grow much larger than ∼ 1 mm and then drift inward too rapidly (cf. §3.3). Since the observed gas accretion rate on TW Hydra sets the constraint αΣ 0 = const (cf. §3.1), we see that the scaling in eq. 14 is This shows that a less massive disc would naturally result in more reasonable maximum grain sizes, avoiding the unwelcome necessity to demand that v br is as small as ∼ 0.2 m/s.

The sharp dust rollover in TW Hydra
TW Hydra displays a cliff-like rollover in the dust density distribution at separation R ∼ − 52 AU (Andrews et al. 2012;Hogerheijde et al. 2016;Tsukagoshi et al. 2019). The observed rollover can be fit with a power-law Σ dust ∝ R −8 at R 50 AU (Hogerheijde et al. 2016). In the context of the standard paradigm for protoplanetary discs, the separation of the dust and the gas may be expected due to the radial drift of the dust and the viscous spreading of the gas disc (e.g., Powell et al. 2019;Rosotti et al. 2019;Trapman et al. 2019).
Let us consider the dust radial drift time scale dependence on distance R: where a(R) is the radius-dependent grain size, and we assumed the low Stokes number limit. Steady-state discs usually have ΣR ≈ const, as we also assumed in eq. 3. Thus a decreasing grain size with increasing R means that t dr increases with R. This in turn implies that dust density gradients are erased over time. Consider two radii in the disc, R 1 < R 2 , and Σ 1 > Σ 2 . Since the drift time scale at R 1 is shorter than that at R 2 , the dust surface density at R 1 drops with time faster than it does at R 2 , and hence the ratio Σ 1 /Σ 2 decreases with time. Therefore, if the TW Hydra dust distribution does evolve from some initial distribution, then that distribution must have had an even steeper rollover than the currently observed ∝ R −8 one. It seems rather contrived to demand such a sharp initial dust edge at time t ≈ 0.

Planet migration constraints
The type I migration time scale at the location of T19 planet is very short: where the dimensionless factor γ (Paardekooper et al. 2010) evaluates to γ = 2.5 in our disc model. The planet migration time is less than 1% of TW Hydra's age. Assuming that we are not observing the system at a special time, TW Hydra should hatch ∼ 100 such planets over the course of its protoplanetary disc lifetime for us to have a decent statistical chance to observe it. This does not appear reasonable; most likely the gas disc is much less massive.

METHODOLOGY
We model the time dependent evolution of dust and gas components in a 1D azimuthally averaged viscous disc with an embedded planet that can optionally lose a part of its mass to the disc.

Gas and planetary dynamics
Our code builds on the work of Nayakshin & Lodato (2012), who modelled the time-dependent evolution of a viscous gaseous disc in azimuthal symmetry, with the disc interacting with an embedded planet via gravitational torque and optionally exchanging mass. Without the mass exchange, the planet modifies the disc structure near its orbit only via these torques; the reverse torques from the disc onto the planet force it to migrate, usually inward. The corresponding equations for the gas surface density Σ are eqs. 36-37 (without the mass exchange term for now), and eq. 45 for the orbital radius evolution (migration) of the planet in Nayakshin & Lodato (2012). We do not solve for the thermal balance of the disc here, assuming that the disc is passively heated and the midplane temperature is given by eq. 2. Our neglect of disc viscous heating is physically reasonable since we consider much larger orbital separations and much smaller stellar accretion rates than did Nayakshin & Lodato (2012). If the planet is massive enough, a deep gap in the gas surface density profile opens up due to the planetary torques on the disc, and the planet then migrates in the type II regime (see Nayakshin 2015). We use a logarithmic bin spacing in radius R from R in = 2 AU to R out = 300 AU with, typically, 250 radial zones. The mass exchange term is not present in the massive quasi-steady disc scenario ( §5) and thus will be discussed later when needed.
The initial gas surface density profile is given by where the exponential rollover R exp = 100 AU, which is reasonable given the observed radial extent of the CO gas disc.

Dust disc evolution
Following Dipierro & Laibe (2017) we extend the code to include the dust component in the disc, although setting = Σ dust /Σ = 0 in their relevant equations as we assume the dust to be in the test particle regime. The time-dependent radial drift and turbulent diffusion equation for the dust is where v dr is the full dust drift velocity given by sum of the radial drift velocity (eq. 5), the additional component due to the gas radial flow and the gravitational torque term from the planet (see eq. 16 in Dipierro & Laibe 2017, for the full expression). D is the turbulent diffusion coefficient for dust, which is related to the turbulent gas viscosity ν = α turb c s H via where the Stoke number is given by eq. 6.

Grain size evolution
Eq. 20 is designed to follow the evolution of dust particles of a fixed size (Stokes number). It is desirable to extended the method to a distribution of grain particles. It is currently prohibitively expensive to model numerically the grain particle size evolution together with the spatial evolution of grains. A physically reasonable approximation, commonly employed in the literature, is to assume that the dust follows a power-law size distribution at all locations in the disc, dn/da g ∝ a −p g , with p = 3.5 for grain sizes between a minimum and a maximum, a min ≤ a g ≤ a. The maximum grain size is allowed to evolve in space and time due to grain growth and collisions (cf. Birnstiel et al. 2012;Vorobyov & Elbakyan 2019;Rosotti et al. 2019). The minimum grain size is of a little consequence at mm wavelengths (e.g., see the top left panel in fig. 3 in Woitke et al. 2019) and so we fix a min = 0.05µm. Note that for p = 3.5, most of the mass is at the largest sizes of dust particles, and the mean grain size is a mean ≈ a/3. By following the dynamics of particles with size a mean we then follow the dynamics of the bulk of the dust (as also done, for example, by Rosotti et al. 2019).
We start the simulations with the maximum grain size being small everywhere in the disc, a(R) = 1µm. We then allow the grains to grow until the maximum grain size reaches either one of the three well known maximum grain size constraints -the turbulent fragmentation, the radial drain, or the radial drift fragmentation limits )as described in §3.4.

Computing the disc emission spectrum
Once we have the dust surface density and the maximum grain size distributions for all disc radii, Σ d and a(R), we compute the dust optical depth where κ(λ, a) is the DIANA dust opacity (Woitke et al. 2016) computed for the maximum grain size a and radiation wavelength λ. Both scattering and absorption opacities are included in κ(λ, a). We used an amorphous carbon fraction of 26%, higher than the standard value used by Woitke et al. (2016). We found this to be necessary to fit the relative luminosities of TW Hydra in 820 µm and 1.3 mm wavelengths. We subsequently found that this is very close to the 25% amourphous Carbon fraction derived by Woitke et al. (2019) for TW Hydra.
We compute the disc surface brightness at radius R as where T(R) is the local disc midplane temperature, and 0 < ζ s ≤ 1 is the dimensionless function describing the disc emissivity reduction due to dust scattering given by eq. 11 in Zhu et al. (2019). As shown by these authors, when τ d 1, eq. 23 reduces to the standard optically thin expression (used by, e.g., Andrews et al. 2016) for the radiation intensity emitted by the disc, which has no scattering contribution. However, when τ d > 1, the scattering albedo may produce a non-trivial and significant reduction of the dust emissivity from the blackbody function B ν ; this reduction is described by the function ζ s . As TW Hydra's disc inclination to us is very small (i ≈ 7 • , cos i ≈ 0.99), we shall for simplicity set i = 0 in this paper.

THE QUASI STEADY STATE DISC SCENARIO
As shown in § §3.2 and 3.3, the gas disc must be be very massive to both feed TW Hya and prevent the mm-sized dust drifting into the star in 10 Myr. We set the initial disc mass to 0.15M to avoid it becoming self-gravitating. For this disc we found that the viscosity parameter of α = 10 −4 yields stellar accretion rate between (1 − 2) × 10 −9 M yr −1 at time between ∼ 3 and 10 Myr, as appropriate for TW Hya. For simplicity, we artificially hold a planet on a fixed orbit at the location of excess emission in T19. For planets with mass M p ∼ 10 M ⊕ , the planet migration time scales are uncomfortably short for massive gas discs, t mig 0.1 Myr (see eq. 18). We explored the parameter space of such more self-consistent models and found that they are challenged by the data even more than the fixed planet orbit models. For brevity we do not show their results here.

10 Earth mass planet at 51.5 AU
Here we present a simulation with the following parameters: an initial disc mass M d = 0.11M (as suggested for TW Hydra by Powell et al. 2019), the initial dust to gas ratio of Σ dust /Σ = 0.03, and the planet mass M p = 10 M ⊕ . The disc viscosity parameter α = 1.5 × 10 −4 was constrained by demanding the gas accretion rate onto the star to be close  to the observed value. The turbulent viscosity parameter α turb was set equal to α. Note that lower values of α turb are unlikely based on the analysis in Dullemond et al. (2018). The maximum grain size in TW Hydra disc is at least 1 mm . To achieve this, the dust breaking velocity had to be set much lower in this simulation, v br = 0.3 m/s, than the value usually assumed in literature (v br ∼ 10 m s−1, e.g., Drażkowska et al. 2014;Rosotti et al. 2019). Fig. 2 shows the model disc intensity at wavelength 1.3 mm (top panel) and the dust disc surface density (bottom panel) at three different times. The observed ALMA intensity and the dust surface density profile deduced by Hogerheijde et al. (2016) for TW Hydra are shown with the solid green curves in the top and bottom panels, respectively. This deduced Σ dust is scaled down by a factor of 3 as the DIANA opacities we use are higher by a factor of a few. The bottom panel also shows the maximum grain size a at time t = 6 Myrs.
We observe a number of features expected from the presence of a massive planet in a disc (Rice et al. 2006;Pinilla et al. 2012;Dipierro & Laibe 2017;Zhang et al. 2018). The planet acts as a dam for the dust, so that a bright outer rim appears behind it. Inside the orbit of the planet, the dust is free to drift into the star, and hence a deep and wide gap appears there. Note that the disc intensity in the top panel is not simply linearly proportional to the dust surface density from the bottom panel. This occurs because the dust grain sizes vary with location in the disc, and even more importantly because the disc intensity saturates at the Blackbody function B ν (T) at very high optical depths τ d . As a specific example, consider radius R ≈ 30 AU in fig. 2. While Σ dust at this radius decreased by almost an order of magnitude going from t = 0.45 Myr to t = 2 Myr, the intensity of the disc emission at that radius did not vary at all.
Overall we see that the model disc intensity is very different from the one observed, and the dust surface densities are also different from the broken power-law fit of Hogerheijde et al. (2016). Although the disc intensity and Σ dust vary with model parameters, in all of the cases we experimented with the model always contradicts the observations: for a sufficiently massive planet, the dust emission should be suppressed inside the planetary orbit and that there should be a bright rim behind it. The observations show no suppression of the dust emission at or inside the orbit of the planet, and the disc intensity does not display a bright rim behind it. Fig. 3 shows a calculation entirely analogous to that shown in fig. 2, except the planet mass is set at M p = 3 M ⊕ . In this case the planet produces only a barely noticeable depression in dust surface density just around its orbit, as the observations demand. The results of this calculation are very similar to that with any smaller planet mass, 0 ≤ M p ≤ 3 M ⊕ . Fig. 3 shows that, with the planet effects on the disc reduced, the model may actually yield a flat emissivity profile in the inner disc followed by a steep decline, exactly as needed to explain TW Hydra's ALMA data. This occurs due to the already mentioned saturation of the intensity in the inner disc where it becomes optically thick. For example, the dash-dotted cyan curves (t = 6 Myr) in fig. 3 appears most promising, with the break in the disc intensity occurring right where needed. However, the saturation of disc intensity at the Blackbody function is also the reason why this scenario contradicts the data strongly. Fig. 4 shows the disc intensity at t = 6 Myr in three wavelengths in the top panel, along with the corresponding ALMA and EVLA observations. The bottom panel shows the respective disc scattering plus absorption optical depth for these wavelengths, the dust surface density Σ dust (which is the same as the cyan line from the bottom panel in fig.  3), and the maximum grain size a. While the model fits the broken power-law shape of the intensity profile of the disc in the two ALMA bands, it is too bright by a factor of ∼ 3. This cannot be "fixed" by any changes in the dust opacity model or variations in σ dust . To understand why, note the purple curve in the top panel of fig. 4 that shows the disc intensity profile in the optically thick limit, i.e., I ν = B ν everywhere. We now see that the break in the disc emissivity profile in the two ALMA bands indeed occurs where its optical  depth exceeds unity somewhat (where the dust absorption only optical depth is ∼ 1).

3 Earth mass planet at 51.5 AU
The only physical way to make the model disc appropriately bright in the ALMA bands is to demand that it becomes optically thick not at R ∼ 50 AU but at R ∼ 30 AU. However, that would contradict the observed intensity profile. Further, 9 mm EVLA data pose a separate but physically similar challenge. To match the correct intensity level at ∼ 50 AU in this wavelength the disc needs to be very optically thin. This then implies that the EVLA emission must track the strong rise in Σ dust inward, but the observed profile is rather flat in the ∼ 30 − 60 AU region.
In fact, it is well known that the outer disc in TW Hydra must be optically thin in sub-mm and longer wavelengths from pre-ALMA photometry (the integrated disc luminosity) data: the luminosity of the source rises as ∝ ν 2.6 rather than ∝ ν 2 as expected for the optically thick disc (e.g., see §3 and fig. 1 in Pascucci et al. 2012). Furthermore, the image of the T19 excess emission is significantly smaller than the disc pressure scale height, and that too implies that the disc is optically thin in the ALMA bands (see §9.1).
Summarising  . Note that while this model matches the location and shape of the rollover well, it over-predicts the disc luminosity by a factor of ∼ 3 in the ALMA bands, as well as contradicting the EVLA data strongly.
TW Hydra's protoplanetary disc. In §9 & §10.1 we shall see that it faces many additional challenges.

A PHENOMENOLOGICAL DUST SOURCE MODEL
We now make a single but significant alternation to the physical setup of our simulations. We assume that the planet ejects dust in the surrounding disc. The simulation setup and initial conditions are exactly the same as those presented in §5.2 except that we assume a negligible amount of dust into the gaseous disc at t = 0 for simplicity. The dust mass loss rate from the planet is a free parameter of the model; we found that choosing M Z = 6 × 10 −6 M ⊕ yr −1 (with other parameters of the model unchanged) provides a somewhat promising spectral model. As in §5.2 we keep the planet mass and orbital radius fixed for now, even though this violates both mass conservation and Newton's second law.  The dust lost by the planet is deposited in a relatively narrow ring with a Gaussian profile around the planet location, with the surface density deposition rate given by where R p is the current position of the planet, and w inj = 0.15R p is the width of the Gaussian. The normalisation constant C ensures that the mass injection rate into the disc is equal to the planet mass loss rate M Z . Through numerical experimentation we found that our results are insensitive to the exact injection profile as long as it is not too broad, and while w inj R p . Fig. 5 shows the disc intensity at 1.3 mm (top panel), Σ dust and maximum grain size a (bottom panel) at several different times. Since initially Σ dust is very low, the dust growth time is long everywhere but accelerates as more and more dust appears in the disc around the planet location. The dust grows to the (observationally required) size of ∼ 1 mm before grain fragmentation stems grain growth. Initially, while the dust is small, dust particles diffuse both inward and outward (e.g., see the t = 1 Myr snapshot   Fig. 4, we note that the present model has a very sharp rollover in I nu behind the planet not because of the opacity transition at that point but because the dust surface density Σ dust (black curve) nose dives at 51.5 AU. The luminosity of this model is closer to what is observed, and can also be scaled down without a significant change in the profile (except for the innermost region where the model disc is optically thick) by a simple reduction in the free parameter M Z . Furthermore, there is a natural casual association between the location of the planet and the rollover behind its orbit in this scenario.
The λ = 9 mm EVLA intensity of the model, on the other hand, is still problematic. Furthermore, this is a phenomenological model that contradicts physics strongly. The planet mass is kept constant at M p = 3 M ⊕ , whereas the dust mass actually present in the disc at t = 5 Myr is M d ≈ 25 M ⊕ . Increasing the planet mass at t = 0 to a value exceeding M d would solve the mass budget problem, but as we saw in §5.1, a planet with mass of ∼ 10 M ⊕ would produce a very deep gap in the dust disc, contradicting the observations. Further, such a massive planet would migrate inward extremely rapidly, e.g., on the time scale of t mig ∼ 10 5 yrs for the disc model used in this section, invalidating the fixed orbit assumption. As the time scale for establishing the quasi-steady state dust distribution in this model is a few Myr (cf. the cyan curves in fig. 6), this is a fatal flaw -the planet ends up in the star faster than this steady state is reached.

Physical motivation
Motivated by the successes and failures of the model presented in §5.2, we now attempt to build a physically motivated model based on the Core Accretion scenario for planet formation (Pollack et al. 1996). Core Accretion scenario planets exist in two physically very different states. After the collapse of the gas envelope around a massive solid core (Mizuno 1980;Stevenson 1982;Pollack et al. 1996), and at the end of the runaway accretion phase, the planet mass is a few M J and its radius is only ∼ 2R J . On the other hand, before the gas accretion runaway, the planet mass is thought to be no more than 30 − 50 M ⊕ , with the solids making up the majority of this mass, and the outer radius of the gas envelope ∼ tens of R J (e.g., see fig. 2 in Mordasini et al. 2012b). Previous models assumed that most of the solids get locked into the core, separating cleanly from the gaseous envelope. However, more recent work (Lozovsky et al. 2017;Brouwers et al. 2018;Podolak et al. 2019) shows that most of solids are vaporised before reaching the core and are suspended as gas in the hydrogen-helium mixture. The opacity of these metal-rich gas envelopes may be significantly higher than that of the traditionally assumed Solar composition ones. Additionally, modern 3D calculations of gas and dust accretion onto cores indicate complex circulating patterns of flows which tend to recycle material from various depths in the planetary atmosphere (Ormel et al. 2015b,a;Lambrechts & Lega 2017;Cimerman et al. 2017). These flows make it harder for the grains to grow and sediment into the core.
Consider now a planet-planet collision energetic enough to actually unbind a Core Accretion planet. Taking a cue from the stellar collisions theory (Benz & Hills 1987), the relative velocity of the two equal mass planets at infinity, v ∞ , must exceed 2.3 times the escape velocity from the surface of the planet, v esc = 2GM p /R p . The required collision velocity to unbind two equal mass planets is hence The circular Keplerian velocity at 52 AU is less than 4 km/s. Collisions of CA post-collapse gas giants will lead to mergers with only a small amount of mass escaping (Benz & Hills 1987); collisions of pre-collapse planets may unbind them. It is also possible that a merger of pre-collapse planet and a massive core will lead to a common envelope like evolution, in which the cores spiral in closer together while unbinding the envelope (e.g., Ivanova et al. 2013). The aforementioned high opacity makes it all the more likely that the energy deposited by the cores in the envelopes will not escape via radiation but will drive the envelope loss.
We therefore explore a model in which a pre-collapse CA planet is a source of dust. This planet may spend a long time (a few Myrs) gaining its significant mass (e.g., see Pollack et al. 1996). In a massive disc studied in §6, such a planet would be lost into the inner disc within a very small fraction of this time due to planet migration. The migration time scale for planets scales as ∝ Σ −1 ∝ M −1 d (eq. 18). Therefore, to make this scenario plausible we must demand the gas disc to be significantly less massive than 0.11M . We pick rather arbitrarily a value of M d = 2 × 10 −3 M while keeping the initial shape of Σ (eq. 19) the same. The results do not depend very strongly on M d . Since the gas accretion rate in the disc is M * = 3παc s HΣ, we must increase the disc viscosity coefficient α to ensure the gas accretion rate remains the same. We thus set α = 2 × 10 −2 . The changes to the values of the disc mass and disc viscosity parameter are important for dust dynamics. A higher α implies that the planet gravitational influence on the disc in its vicinity is significantly reduced. More massive planets may be present in the disc without opening a deep gap that would contradict observations.

Numerical results
Figs. 7 & 8 show the results for a simulation started with initial planet mass M p = 38 M ⊕ . For simplicity we assumed a uniform planet composition, with metallicity Z = 0.5. Both gas and dust are injected into the disc at a constant (and equal because Z = 0.5) rate of M Z = 2 × 10 −4 M ⊕ yr −1 until the planet mass drops to 10 M ⊕ 4 . The planet mass is reduced accordingly as it loses mass. To exemplify the weak dependence of our results on the exact dust injection profile, the width of the Gaussian is here set to w inj = 0.05R p (cf. eq. 24; this is three times narrower than in §6). Unlike the phenomenological model of §6, the planet is free to migrate, but on the account of the low disc mass it migrates very little during this simulation, from the starting radius of R = 53 AU to 51.5 AU.
On the whole we see that the dust profile, and the resulting disc intensity in the three wavelengths, is quite similar to those obtained in the phenomenological massive disc model (figs. 5 & 6). As in the latter model, the disc emissivity has a very sharp -in fact too sharp compared with the observation -rollover behind the orbit of the planet. The similarity in the results despite the difference in the gas disc mass of a factor of 50 between the two models shows that there is a certain degeneracy in the model parameters, e.g., a higher value of α could be compensated for by a higher v br .

Physical motivation and constraints
In §7.1 we argued that a Core Accretion planet may lose a major fraction of its gas-dust envelope if two conditions are satisfied: (i) the envelope is in the extended, pre-collapse state; (ii) a significant energy is injected in it, e.g., via collision with another massive core. We now detail conditions under which a GI planet disruption could be relevant to TW Hydra's observations.

The need for a very rapid primordial disc dissipation
In the Gravitational Instability (GI; Kuiper 1951) theory for planet formation, massive and very young gaseous discs fragment at separations ∼ 50 − 100 AU onto Jovian mass gas clumps (e.g., Rafikov 2005;Rice et al. 2005). Hydrodynamical simulations show that these planets migrate inward very rapidly in massive discs (e.g., Vorobyov

Chauvin et al. 2015)
. On the other hand, planet-planet scatterings may allow some GI planets to survive on wide orbits (Vigan et al. 2017), especially if the primordial disc is dispersed rapidly. For the early massive protoplanetary discs, the primary timescale on which its mass is lost (Clarke et al. 2001) is the viscous time, where α −1 = (α/0.1) and the estimate is made at R = 51.5 AU. Numerical simulations show that the α parameter due to self-gravity of the disc may reach values of order ∼ 0.1 in early massive discs (Gammie 2001;Rice et al. 2005;Haworth et al. 2020), and even α ∼ 0.2 in the magnetised discs (Deng et al 2020). Additionally, disc depletion due to external photo-evaporation may be faster than previously thought (e.g., Haworth & Clarke 2019).
For the case at hand we must require that the primordial gas disc is long gone in TW hydra. This is because the migration time scale of a GI planet with mass initially exceeding 1 M J would be much shorter than 10 Myr. Indeed, if the planet did not open a gap and migrated in the type I regime then its migration time is less than 1 Myr even for a disc as low mass as a few M J . If, on the other hand, the planet did open a wide gap and migrated in type II then the migration time scale is (e.g., Lodato & Clarke 2004) where we used M p ∼ 2 M J (this will be justified later) and M * ∼ 2 × 10 −9 M yr −1 . Therefore, the planet would have been long lost into the star if the disc was there for the last 10 Myr. We emphasise the distinction with the CA problem setting discussed in §7 brought about by the different planet formation mechanisms. In the CA scenario the planet does not need to be born at t ≈ 0. As is well known, in the classical CA model massive cores are most likely to be made at late times, e.g., at t ∼ 5 − 10 Myr (Ida & Lin 2004b;Mordasini et al. 2012a) since the process of core growth is slow. Further, due to its lower mass the type I migration time scale is longer. Therefore, there is no reason to demand a complete disappearance of the primordial gas disc before the planet disruption commences in the CA framework.

Why did the planet not collapse in 10 Myr?
GI planets are born extended, with their radius r p ∼ 1 a few AU, and cool: their central temperature is in hundreds of K (e.g., ). If dust growth inside the GI planet is neglected, then it cools, contracts, and eventually collapses dynamically when the endothermic reaction of H 2 molecule dissociation absorbs a vast amount of thermal energy of the planet (Bodenheimer 1974). The collapse terminates in formation of a planet that is ∼ Million times denser, with radius R ∼ 2R J and an effective temperature of 1, 000 to 2, 000 K. This luminous post-collapse state is often called the "hot start" of gas giant planets (Marley et al. 2007). Similar to the post-collapse gas giant CA planets, the post-collapse GI planets are unlikely to lose mass at ∼ 50 AU from the star.
Hence we must demand that the planet remains in the pre-collapse state before the onset of the mass loss. This is surprising given the age of the system. The evolution from birth to collapse (hot start) is usually thought to be very fast. This result is rooted in the pioneering work of Bodenheimer (1974) who found planet collapse time scales 0.5 Myr for M p = 1 M J , and even shorter for higher mass planets. However, the collapse time scale is sensitive to the dust opacity model used. More recent dust opacity calculations (e.g., Semenov et al. 2003;Zhu et al. 2009;Woitke et al. 2016) indicate that dust opacity may be higher by up to a factor of ∼ 30 at T ∼ 10 − 30 K (the effective temperature of GI protoplanets) compared to the opacity employed in the 1980s (e.g., Pollack et al. 1985, see Appendix A). We find that these higher dust opacities lengthen the duration of the pre-collapse phase by a factor of 5 − 10. Further, 3D simulations of GI planets immersed in protoplanetary discs show that these planets accrete pebbles very rapidly and become significantly enriched in dust ( The time scales on which the core grows are a minimum of thousands of years but may be much longer, depending on convection and grain material/growth properties, such as v b ). If the core grows more massive than ∼ 10 M ⊕ , then the energy release during its formation can be too large for the pre-collapse planet -its envelope expands and is eventually lost (Nayakshin & Cha 2012; Nayakshin 2016; Humphries & Nayakshin 2019). This scenario for the core-initiated disruption of GI protoplanets is physically analogous to how cores of AGB stars eject their envelopes except for the energy source -the gravitational potential energy rather than the nuclear energy of the core -and the physical scales of the systems. At present, there exists no stellar/planet evolution code that takes into account all the relevant physics that we wish to explore here. For example, ; Helled & Bodenheimer (2011) present models of grain growth and sedimentation in pre-disruption isolated planets cooling radiatively. Vazan & Helled (2012) investigate how external irradiation affects contraction of these planets. These studies did not include the effects of the massive core energy release onto the planet, which is central for us here. On the other hand, Nayakshin (2015Nayakshin ( , 2016 include grain growth, sedimentation, core formation and the effects of the core energy feedback onto the envelope, but use a simplified followadiabats approach to model radiative cooling of the envelope, and a simpler equation of state than Vazan & Helled (2012) do. Further, the opacities used by the two codes are different.
Here we shall use the code of Nayakshin (2016) to understand the physical constraints on the pre-disruption planet mass, metallicity, and the mass of the core responsible for the planet disruption. These constrains will be seen to place significant limitations on the disrupted GI planet scenario (e.g., the planet mass is unlikely to exceed ∼ 2 M J ). In appendix B we compare for the first time the results of uniform planet contraction calculations (no dust sedimentation allowed) computed with this code with that of the proper stellar evolution model of Vazan & Helled (2012) at the same (Pollack et al. 1985) dust opacity. We find that the difference in the planet collapse time scale computed by the two codes is within a factor of two, which we deem sufficiently close given the much larger dust opacity uncertainties ( §A).
The thick dashed curves in the top panel of Fig. 9 show the evolution of the radius r p of a 2 M J planet circling the 5 This conclusion holds as long as grain growth and settling do not deplete a < ∼ 100µm population of grains. At higher grain sizes, Rosseland mean dust opacity may drop (cf. fig. A1). In that case higher Z planets may actually cool more rapidly (Helled & Bodenheimer 2011).  star with TW Hydra properties at 60 AU for different planet metallicities, from Z = 1 (in units of Z = 0.015) to Z = 12. All the models start with the central planet temperature of 200 K, the uniform composition and an initial grain size of a = 0.01 mm. The Zhu et al. (2009) opacity table is used for this calculation. The dust opacity is assumed to be proportional to the metallicity Z of the envelope. The bottom panel of fig. 9 presents the core mass versus time. The grain breaking velocity is here set at v br = 5 m/s. The thin curves in the top panel show the same calculations but where grain growth and core formation are artificially suppressed. Fig. 9 shows that at Solar metallicity, Z = Z , the planet contracts and collapses by t ≈ 1 Myr, whether core formation is allowed or not. The planet evolutionary time scale increases as the metallicity of the planet increases, and so does the core mass. Formation of the core in the planet speeds up evolution of the planet in all the cases. This occurs due to a lower dust opacity in the envelope as some of the dust settles and gets locked in the core. At the lower metallicities in the figure, the core masses are relatively low, so the effect of the core formation is negligible save for the dust opacity reduction. However, for metallicities Z ≥ 8Z the core mass exceeds 20 M ⊕ . The gas envelope starts to expand soon after this mass is reached and is eventually unbound. For the Z = 8Z model, the envelope is disrupted at about 8 Myr in this calculation. A higher metallicity planet meets its end sooner, at about 4 Myr, as its core is more massive and more luminous.
We have found that planets more massive than ∼ 2 M J are not likely to be disrupted by their cores. This comes about due to two factors. First of all, even at a fixed bulk composition, more massive planets contract much more rapidly, shortening the time window for the core growth (this effect exists whether the core feedback is included or not, see, e.g., . Secondly, observations (Miller & Fortney 2011;Thorngren et al. 2016) show that more massive planets are less metal enriched than their less massive cousins. Simulations of pebble accretion onto gas giant planets (Humphries & Nayakshin 2018) also lead to the same conclusion. The dust opacity of the massive planets is hence expected to fall with planet mass, exacerbating the challenge of assembling a massive core and disrupting the planet with it.

Deposition of matter in the secondary disc: methodology
There are several free parameters in this model (just like for the model in §7) which we constrain by trial and error. In the beginning of the calculation, we specify the initial planet mass, M pi , and the starting position of the planet, R pi . As per §8.1.3, we use a GI planet with initial mass M pi = 1.5 M J . By experimenting we found that the disc viscosity parameter α = 2.5 × 10 −2 results in gas accretion rate similar to the one observed in this system. Similarly, setting R pi = 56 AU resulted in the planet remnant being stranded at 51.5 AU; the planet initial bulk metallicity of Z = 6Z gave the right ALMA luminosity for the disc. We assume that at the time of disruption the planet contains a massive solid core or at least a region so metal rich that it survives the disruption of the hydrogen-rich atmosphere. We refer to the final planet mass as simply the core, and it is set to M c = 0.03 M J ≈ 10 M ⊕ for definitiveness here. The part of the planet that is lost and injected into the protoplanetary disc is termed the "envelope", and its mass is The mass loss rate, M p (t), is not known a priory. As for binary stars undergoing mass exchange (e.g., Rappaport et al. 1983;Ritter 1988), it is a function of the planet internal structure and its orbital evolution that in itself depends on the planet-disc interaction and the mass loss rate (see Nayakshin & Lodato 2012). Such a fully self-consistent calculation is beyond the scope of the current paper, and we instead specify the planet mass loss rate: where t dis = 10 5 yrs. This mass loss rate is partitioned between that for gas (H and He) and metals (mass fraction Z), thus In general we do not expect the envelope to have a uniform composition, Z(M), since dust is of course able to sediment down through the gas. Hence we expect Z(M) to be a function that decreases from the maximum in the core, which we simply set Z(M p < M c ) = 1, to some minimum. For definitiveness, we choose this functional form: The parameter M tr M e describes how large the metal rich region in the centre of the planet is, and Z ∞ is the metallicity at the outer reaches of the envelope where the term exp(−M e /M tr ) ≈ 0 as we use M tr M e . In practice, we specify the mean metallicity of the planet envelope,Z, and M tr , from which Z ∞ is computed. For the calculation presented in §8.3, M tr = 0.03 M J .
The dust and gas lost by the planet are deposited in a Gaussian ring around the planet location, as described by equation 24. The normalisation constant C ensures that the mass injection rate into the disc is equal to the planet mass loss rate (eq. 28). As per eq. 29, the injected mass is split into gas and dust. The gas surface density (thin lines) evolution shows the dominant features of the well known viscous "spreading ring" calculation modified by the continuous mass loss from the planet. The gas spreads quickly all the way to the star and to R ∼ 200 AU, as required by the observations of gas accretion and the large extent of the gas disc in TW Hydra (cf. §2). Despite a continuous mass injection into the disc, the planet manages to make a depression in the gas surface density profile around its orbit due to gravitational torques acting from the planet onto the gas. This effect is noticeable while the planet mass is in the gas giant planet regime. By t ≈ 0.1 Myr the planet has lost too much mass (the remnant mass M p = 10 M ⊕ ) to affect the gas surface density profile gravitationally at this relatively high value of α for our disc. However, the GI planet legacy lives on in the form of the significant break in gas Σ profile; the break is coincident with the planet location. Physically, the break appears because gas flows inward towards the star inside the orbit of the planet, and outward outside the orbit.

Numerical results
The thicker lines show the dust surface density profile at the respective times. We see that initially the dust surface density is narrower than that of the gas, but eventually the dust spreads. This spread is mainly inward of the planet. As in §7, at late times the dust dynamics is dominated by the radial drift: Once ejected by the planet, large dust particles are blown inward of the planet by the aerodynamical friction. The resultant dust surface density profile at t ≈ 0.1 Myr is qualitatively similar to the broken powerlaw fit (the thick green dashed line) used by Hogerheijde et al. (2016) to fit TW Hydra's ALMA dust continuum intensity profile in Band 7. The upturn in Σ just inward of the planet is due to the assumed dust composition profile within the planet (eq. 30) in which the dust concentration near the core is far greater than at the outer edge. Fig. 11 shows the resultant disc emissivity profile in the three wavelengths in the top panel and the disc properties in the bottom panel. The model fits the data reasonably well except for the rollover, which is too sharp, just like the model in §7.2. This model is optically thin, as is observationally desired. The total mass of the dust in the disc in this model is about 15 M ⊕ , much smaller than ∼ 80 M ⊕ estimated by Andrews et al. (2016) and also smaller than ∼ 35 M ⊕ estimated by Woitke et al. (2019). This is mainly due to the different dust models used in these studies. The standard DIANA dust opacity (Woitke et al. 2016) are higher by a factor of a few than that used by Andrews et al. (2016); Hogerheijde et al. (2016). Although we use the public DI-ANA opacity code to compute our dust opacities here, we use the standard a −3.5 dust size distribution whereas Woitke et al. (2019) found that a −4 power law was a better fit in their modelling.

THE ALMA IMAGE OF T1FEATURE
Here we discuss the implications of the 2D morphology of the excess emission observed by Tsukagoshi et al. (2019). We use these as additional probes of the scenarios explored in this paper.

The disc is optically thin
As found by Tsukagoshi et al. (2019), the excess has a radial half width of ∼ 0.5 AU and an azimuthal half-width of ∼ 2.2 AU. The disc pressure scale height for TW Hydra is H ∼ 3.5 AU at separation of 51.5 AU, and thus the observed feature is significantly smaller than H. We note that this immediately implies that the disc (but not necessarily the feature) is optically thin at 51.5 AU. This is because photons emitted from the disc midplane perform a random walk in an optically thick disc until they escape vertically out of the disc. Therefore, an image of a point source placed in the midplane of such a disc would be broadened by ∼ H at least. This forms an independent confirmation, in addition to the arguments spelled out in §5.2 that TW Hydra's disc is optically thin and thus the rollover in Σ dust observed behind the T19 feature could not be due to the disc becoming optically thin at this radius.

A vortex or a circum-planetary disc?
Vortices (e.g., Li et al. 2001) have been suggested to trap dust material in the protoplanetary discs and thus produce bright excess in the dust continuum emission (Baruteau & Zhu 2016). Furthermore, vortices are azimuthally elongated structures, with axis ratio ∼ 4 − 6 (Richard et al. 2013), exactly as observed. However, the radial half-width of vortices is at least H, and likely ∼ twice that (e.g., fig. 3 in Lin 2012). For the same reason a vortex would also look much more extended in the azimuthal direction (see fig. 4 in Baruteau & Zhu 2016) than observed. Just like with a gap edge created by a planet, we expect a hole in the dust density distribution inward of the vortex ( fig. 10 in Baruteau & Zhu 2016), which is not observed.
Circum-planetary discs are believed to be at most ∼ 1/3 (Bate et al. 2003;Ayliffe & Bate 2009), and more likely 1/10 of the planet Hills radius, as shown by the more recent higher resolution calculations (Wang et al. 2014;Ormel et al. 2015b). Thus, the planet would have to exceed the mass of ∼ 1.5 M J to account for just the radial size of the feature. This does not account for the much larger azimuthal extent of the T19 excess emission. Such a high planet mass would produce a noticeable gap at the disc gas surface density and the dust intensity profile near the planet even for disc viscosity as high as α = 0.025. This is not observed. Finally, Tsukagoshi et al. (2019) also points out that the total flux from the circum-planetary disc is insufficient to account for the total flux in the feature.

The dust trail of a planet losing mass
Here we consider the dynamics of grains lost by a low mass planet on a circular orbit embedded in a laminar gas disc around it. The most likely physical origin for a dust particle outflow from a planet is a thermally driven gas outflow that picks and carries the dust with it. The dynamics of grains in the planet vicinity, at radii between the planet radius and the Hills radius, r p r R H , clearly deserves a separate detailed investigation which is outside the scope of this paper. Here we are concerned with how the flow may manifest itself to ALMA on scales of much larger than r p . We perform a 2D calculation of dust particle orbits assuming their trajectories lie in the midplane (note that R H H for a low mass planet).
The dust particle size at the outflow is likely to be much smaller than the mm-sized particles that ALMA sees in TW Hydra's disc. This is because the inner region of the planet is expected to be sufficiently hot to vaporise even the most refractory dust (Brouwers et al. 2018). This is relevant because for both the CA and GI planet losing mass scenarios the time when the model fits the data best is close to the end of the mass loss phase from the planet, when the most central regions of the planet are lost into the disc.
As the outflow leaves the planet, the gas density drops, and so does its optical depth. Due to adiabatic expansion and radiation (the outflow eventually becomes transparent to radiation) the gas temperature drops with distance from the planet and the metals re-condense into dust particles. The grains then grow rapidly. In the disc geometry, the grain growth time scale is t grow ∼ (1/Ω K )(Σ/Σ dust ) . For the gas just lost by the planet the gas-to-dust ratio is not very large as the planet central regions are very metal rich in our model, therefore t grow may be expected to be of order a few orbital times in the disc, 2π/Ω K .
The dynamics of dust particles is most sensitive to the Stokes number, St, and hence we reformulate the problem in its terms. The dust particles are ejected by the planet with initial Stokes number St min = 10 −3 , and grow to a maximum size corresponding to the maximum Stokes number of St max = 0.1. We describe the particle growth process as a time-dependent St number where q = t /(t grow + t ), where time t counted from the time the grain was ejected by the planet, t grow = N g (2π/Ω) is the growth time scale, with N g = 4. For reference, The planet is assumed to be of a sufficiently small mass that we can neglect its dynamical influence on the surrounding gas. Similarly, we neglect the planet physical size compared with the scales of interest, assuming that dust is emitted from a point (planet position). The planet is on a circular orbit around the star at R = 51.5 AU. We integrate the standard 2D equations of motion for individual dust grains lost by the planet. Grains are ejected by the planet at a steady rate. As expected, the grains are blown inward by the radial drift and eventually disappear into the star but here we are interested in the dust particle morphology in the immediate vicinity of the planet to compare to the T19 ALMA image of the excess emission region.
The left panel of Fig. 12 shows the location of the dust particles integrated as described above. The coordinates are centred on the planet which is on a prograde circular orbit that follows the dashed circular path. We observe that the dust particles perform a U-turn as seen from the planet location. Since initially the dust particles are small ( St 1), they are very tightly coupled to gas at t t grow . In the frame of the planet they start to lag behind the planet because the dust is picked up by the gas in the disc and so travels with the velocity which is smaller than the planet orbital speed v K . However, as the dust particles grow, they start to drift inward of the planetary orbit. When the dust particle drifts to radius R < R such that its angular speed there exceeds that of the planet, that is, v φ (R )/R > v K (R)/R = Ω K (R), it starts to orbit around the star faster than the planet. Hence the particle overtakes the planet eventually.
This orbital motion of the dust looks like a tight Uturn around the planet. Additionally, as the particle sizes grow as they move farther and farther away from the planet, the speed differential between the planet and the dust particles increase. This leads to dust particles spending more time in the U-turn region than in the region in front of the planet.The dust density is hence larger in close proximity to the planet and decreases with distance along the dust tail.
The right panel of fig. 12 shows a simulated ALMA image of the dust dynamics smoothed by the ALMA beam for TW Hydra. We added a point source to the planet location and a uniform background. We see that the U-turn of the dust produces an emission feature elongated in the azimuthal direction, somewhat analogous to the observations shown in fig. 1. The exact brightness of the extended tail compared to the U-turn region, the pitch angle of the tail with respect to the azimuthal direction, and the length of the U-turn region do depend on the parameters on the dust growth model used. However, a good qualitative match to the shape of the Tsukagoshi et al. (2019) excess emission is obtained for a wide range of model parameters, not requiring fine tuning. The weak excess emission in the dust tail running in front of the planet and pointing inward of its orbit is a unique testable prediction of our model.

DISCUSSION
TW Hydra hosts a protoplanetary disc resolved by ALMA with a record breaking resolution of about 2 AU. The disc has a cliff-like rollover beyond about 50 AU and a significant excess emission resolved by ALMA into a blob with sizes ∼ 1AU radially by ∼ 4 AU azimuthally. The excess is suspected to be a young planet and is located at 51.5 AU, right at the edge of the dust disc. The protostar continues to accrete gas at a respectable rate of ∼ 2 × 10 −9 M yr −1 despite being one of the oldest protstars known (∼ 10 Myr old). Furthermore, its dust disc is nearly two orders of magnitude more massive than the median for Class II sources (Williams et al. 2019), most of which are younger than TW Hydra. Here we have shown that we can use these known and unique properties of the system to constrain scenarios of the protoplanetary disc evolution.

The quasi Steady State scenario
In this scenario ( §5) the protoplanetary disc in TW Hydra is primordial, e.g., ∼ 10 Myr old. Previous work (e.g., Powell et al. 2019) and analytical arguments on the presence of ∼ mm-sized dust in the disc ( §3.3), and the observed gas accretion rate ( §3.2), require the gas disc mass M d to exceed ∼ 0.1M in this case. We found this scenario to be challenged by the data and other results in the field: (i) The observed cliff-like rollover in the dust continuum emission beyond ∼ 50 AU is very puzzling for such an old disc that is known to extend to ∼ 200 AU in CO and other molecular tracers. Numerical experiments ( §5.2) and analytical arguments ( §3.5) show that one expects a powerlaw like decline in Σ dust with radius in this case. We found that the scenario may produce a sharp break in the disc emissivity if the disc becomes optically thick inside 50 AU ( fig. 4). However, for TW Hydra this over-predicts the total flux from the disc by a factor of ∼ 3, contradicts earlier conclusions from photometry data ( §5.2), and the relatively small size of the Tsukagoshi et al. (2019) feature ( §9.1). All of these observations require an optically thin disc in the ALMA bands. Further, this model disc becomes optically thin at longer wavelengths and hence predicts a steep decline in the disc intensity with radius at ∼ 50 AU, whereas the 9 mm EVLA data show a very gradual decline in that region (blue curves in fig. 4). scale for these planets is ∼ 1 to a few% of TW Hydra age ( §3.6). To observe even one of these three putative planets, we need to be quite lucky. To observe all three of these planets, we need to postulate that they were all born essentially simultaneously. This is unlikely because the rates of planet embryo assembly are strongly separation-dependent (Ida & Lin 2004a;Mordasini et al. 2012b;Lambrechts et al. 2014;Ndugu et al. 2019).
(iii) The association of the Tsukagoshi et al. (2019) feature with the dust disc rollover. Planets are expected to block the inward flow of dust, producing gaps at the location of the planets, and bright rings just beyond the gaps (Rice et al. 2006;Pinilla et al. 2012;Dullemond et al. 2018;Zhang et al. 2018). Our numerical experiments in §5.1 confirmed these well known results. These dust emission characteristics are not observed in TW Hydra, where the emission plunges instead of rising beyond the planet location.
(iv) The size and shape of the extended excess emission detected by Tsukagoshi et al. (2019). The half-sizes of the feature at 1.3 mm continuum is ∼ 2.2 AU in the azimuthal direction and ∼ 0.5 AU in the radial emission. This extent is too small for a vortex but too large for a circum-planetary disc of a Neptune mass planet ( §9.2). Increasing the planet mass above 1 M J may result in the disc large enough but its elongation along the azimuthal direction by a factor of 4 contradicts numerical simulations of circumplanetary discs which show no such elongation (Ayliffe & Bate 2009;Wang et al. 2014;Ormel et al. 2015a;Szulágyi et al. 2014). Additionally, such a high mass planet is ruled out based on the protoplanetary dust disc morphology as explained in (iii).
(v) A very small disc viscosity. The observed gas accretion rate onto the star is surprisingly low if the disc mass is really as high as 0.1M , and requires the disc viscosity parameter α 10 −4 (see §3.3). The results of modelling ALMA observations of other bright discs with annular gaps and rings show that for particle sizes a > ∼ 2 mm the disc viscosity parameter must be larger than ∼ 10 −3 ( fig. 7 in Dullemond et al. 2018).
(vi) Surprisingly small grains. Analytical arguments ( §3.4) and numerical models show that in a disc as massive as 0.1M the grains should grow rapidly to sizes much larger than ∼ a few mm, and be lost into the star too soon, unless the dust fragmentation velocity is ∼ 0.3 m/s. The latter value is significantly lower than the values obtained in laboratory experiments (e.g., Blum & Wurm 2008) and ∼ 10 m/s typically used in the protoplanetary disc literature (e.g., Birnstiel et al. 2012;Drażkowska et al. 2014;Rosotti et al. 2019).
(vii) Dusty rings rather than spirals. Veronesi et al. (2019) show for TW Hydra and other ALMA discs that gas discs can be "weighted" by understanding the response of the mm-sized grains to the planets embedded in these discs. They find that in massive discs the mm-sized grains would tend to be in spiral features driven by the planets, whereas in low mass gas discs they would conform to the shape of rings. Based on the absence of spirals and presence of rings the authors conclude that the disc masses are ∼ 1 M J .
(viii) Oddity compared to other discs. The dust mass of TW Hydra is extraordinarily large. Using the pre-DIANA opacity model, this mass is estimated at ∼ 80 M ⊕ (Andrews et al. 2012. The mean dust mass of class 0 sources was recently estimated at 26 M ⊕ (Tobin et al. 2020); these sources are a factor of ∼ 100 younger than TW Hydra. The more comparable yet still younger by a factor of several class II discs have dust disc masses almost 2 orders of magnitude lower than that of TW Hydra (Williams et al. 2019). Using DIANA (Woitke et al. 2016) opacities we find a factor of ∼ 6 lower dust mass for TW Hydra; however all the other results cited above should then be scaled down as well, leaving TW Hydra's dust mass excess just as large.

The phenomenological planet losing dust scenario
In §6 we used the same massive disc scenario as discussed in §10.1, but assumed that the protoplanetary disc is dust-free before a source of dust of unspecified nature starts ejecting dust. This produced a better match to the data, resolving qualitatively the problems listed in (i). Large dust particles drift inward, naturally explaining why the T19 feature is positioned at the dust disc rollover, alleviating (iii). However, the model violates mass conservation and second Newton's law, and does not resolves the other challenges from §10.1.

A Core Accretion planet losing dust
We argued in §7.1 that a massive core can lose its precollapse massive dusty gas envelope if a catastrophic release of energy occurs in its core, due to e.g., a merger of the core with another massive core. We dropped the assumption of a massive gas disc, which becomes unnecessary if the dust in the protoplanetary disc of TW Hydra is of a recent rather than primordial origin. In the particular example of the numerical calculation in §7.2 we considered a pre-collapse planet of the total initial mass of 38 M ⊕ to lose its half dust/half gas envelope (planet metallicity Z = 0.5) until its mass dropped to 10 M ⊕ , at which point the mass loss was turned off. The initial disc mass was set at 0.002M , which required disc viscosity of α = 0.02 to yield the correct gas accretion rate onto the star. This vastly improved the results, resolving all issues (i)-(viii), e.g., producing a reasonable match to the observed spectra with a now reasonable value for the grain fragmentation (breaking) velocity, v br = 10 m/s, a value of α in accord to the constraints from DSHARP modelling (Dullemond et al. 2018). Due to the much lower gas disc mass, the planet migration time is comparable to the age of the system, not requiring a miracle of several planets being born at the same time. The Stokes number of mm-sized grains satisfied the Veronesi et al. (2019) constraint. A disruption of a massive pre-collapse planet via a catastrophic collision with another planet is not likely to be a common outcome for the Core Accretion scenario, and this may explain why TW Hydra is such an oddity (viii). Further, recently Demidova & Grinin (2019) showed that catastrophic collision of planetary embryos in a protoplanetary discs releases enough dust to be observable with ALMA.

A Gravitational Instability planet disruption
In §8 we considered a disruption of a gas giant planet formed by the GI scenario. Since GI planets are presumably born in very young, class 0/I discs, this planet would have survived at such a wide orbit only if the disc was dissipated very rapidly in this system ( §8.1.1). This therefore requires that TW Hydra had no protoplanetary disc before the planet disruption. It is also possible that there were more GI planets early on, and that one of them was scattered on a wider orbit than the rest, boosting its chances of survival far out.
Just as with the Core Accretion planets, the pre-collapse GI planets are extended and are susceptible to their envelopes being destabilised if enough energy is injected into the planet centres. Massive solid cores (M core 10 M ⊕ ) were previously shown to be capable of disrupting the planet envelopes (e.g., Nayakshin & Cha 2012;Nayakshin 2016;Humphries & Nayakshin 2019). For a very old system such as TW Hydra, we found that only very metal rich (Z ∼ 0.1) GI planets with masses no larger than 2 M J can be disrupted via this mechanism at t ∼ 10 Myr ( §8.1.2 and §8.1.3).
In §8.3 we found that disruption of a planet with an initial mass M p = 1.5 M J , initially orbiting TW Hydra at 56 AU resulted in a gas disc quickly spreading both inward, to fuel gas accretion onto the star at rates close to those observed in the system, and outward to ∼ 200 AU. As with the CA planet losing dust scenario, the dust lost by the planet grows to mm sizes and then streams only inward of the planet due to the aerodynamical friction with the gas.
Although the physics of the models differs, the two scenarios give equally promising explanations for the observed spectra of the source and yield attractive explanations to all the points (i) -(viii) raised as difficulties of the standard quasi-steady state scenario.

The dust morphology of the T19 excess emission source
Finally, we investigated the dynamics of dust grains lost by a low mass planet in 2D in §9. We argued that dust particles must be carried away from the planet by a gas outflow, and must therefore be microscopic initially. We then argued that when released into the disc the dust will grow to larger sizes as constrained by the disc properties. The dust particles were found to perform a U-turn around the planet, first being dragged along by the gas flowing past the planet, but then overtaking the planet a little inward of its orbit when they have grown sufficiently to drift through the gas. When convolved with the ALMA beam at 1.3 mm this results in an emission excess elongated along the orbit and predicts a weak tail extending in front and a little inward of the planetary orbit. GI disrupted planet vs previous models: disc properties gas Trapman+17 gas Woitke +19 gas Gorti +11 ALMA T19 Planet dust, H+16 x 1/3 gas, t= 9.8E+4 yr dust, t= 9.8E+4 yr Figure 13. Comparison of the gas (red dotted) and dust (blue dotted) disc surface density profiles of the GI planet disruption model ( fig. 11, §8) with that of previous authors for TW Hya. While our model matches the dust density profile from Hogerheijde et al (2016) reasonably closely, our gas surface densities are lower by a factor of a few than Trapman et al (2017). Fig. 13 shows the model dust and gas surface density profiles from fig. 11 that we found to match the observed disc spectra best. These are compared with the broken power-law dust surface density model of Hogerheijde et al. (2016) and the three models for Σ previously shown in fig. 1. We see that while the dust surface density match is reasonably good (which of course is the goal of our paper), the gas surface density profile and the total gas disc mass is significantly lower than previous workers assumed or derived. Without sophisticated modelling it is unclear how serious the disagreement is. Since the disrupted planet is metal rich in our scenario, the mass of various molecular species may be sufficiently high in the model to account for their observed emission since previous workers assumed much lower abundances for the disc. However, HD line emission is not expected to be sensitive to the metallicity of the gas (e.g., Bergin et al. 2013;Woitke et al. 2019); it remains to be seen whether our much less massive disc may account for the observed line fluxes.

CONCLUSIONS
Here we focused on the first ALMA 1.3 mm dust continuum excess emission (Tsukagoshi et al. 2019) positioned right at the edge of a cliff-like rollover of the dust disc in TW Hydra. We showed that the morphology of the blob-like excess and its relation to the dust disc are best explained by a planet losing dust and gas within the excess. We argued that pre-runaway Core Accretion planets and pre-collapse Gravitational Instability planets may be disrupted and may provide the required mass injection into the system. This catastrophic event may also explain why there is a factor of ∼ 100 more dust in this very old system than the mean for (typically younger) class II protoplanetary discs. Future modelling needs to improve on the internal planet structure, mass loss dynamics, dust composition and opacity, and chemodynamical modelling. /g] Rosseland dust opacities at Solar metallicity, Z = Z (0.015) Pollack +85 Semenov +03, spheres Semenov +03, aggregates Zhu +09 DIANA amax_10mu DIANA amax_100mu Figure A1. Rosseland opacity as a function of temperature for different dust opacity models as indicated in the legend. Fig. A1 shows several models for the absorption Rosseland mean opacity, κ R , of gas at Solar metallicity from several different authors as a function of gas temperature. For temperatures T ≤ 10 3 K and gas densities typical of pre-collapse planets, κ R is strongly dominated by dust. For the present paper, the shaded temperature region, 10 ≤ T ≤ 30 K, is most important as this encompasses the expected effective temperatures for our dust-rich wide-separation planets (the energy transfer in deeper hotter planet interiors is dominated by convection anyway ). The two DIANA opacity calculations (Woitke et al. 2016) neglect grain vaporisation as this is not important in the shaded region, but include the effects of grain growth, by allowing the maximum grain size to be either 10µm or 100µm.

APPENDIX A: DUST OPACITY UNCERTAINTIES
We see that there is a factor of about 30 uncertainty between the smallest and the largest κ R . This shows that early calculations of giant planet contraction (e.g., Bodenheimer 1974) may have significantly over-estimated the luminosity of these objects. Furthermore, higher metallicity objects have proportionally higher dust opacities, further delaying planet contraction.

APPENDIX B: PRE-DISRUPTION PLANET CONTRACTION COMPUTED WITH TWO DIFFERENT CODES
As explained in §8.1.3, to model planet contraction simultaneously with dust growth and sedimentation into the core, we use the code of Nayakshin (2016). Here we compare the results of this code, which uses an isentropic (followadiabats) approximation to the energy transfer through the planet envelope, to the more accurate stellar evolution model of Vazan & Helled (2012) for the simpler case in which grain growth and sedimentation are neglected. Fig. B1 shows the evolution of planetary radius computed with the two different codes for the same opacity (Pollack et al. 1985) for several planet masses. The evolutionary tracks computed with the two codes are within ∼ 30% of each other in terms of   (2016) and Vazan & Helled (2012). For each planet mass, the two models are initialised with the same central temperature.
the absolute value of the planet radius, and within a factor of two in terms of the planet collapse time scales. We deem this sufficiently close given the much larger uncertainty that exists in the dust opacity.