Ejection and Disruption of Giant Planets

Consequences of the Ejection and Disruption of Giant Planets


The discovery of Jupiter-mass planets in close orbits about their parent stars has challenged models of planet formation. Recent observations have shown that a number of these planets have highly inclined, sometimes retrograde orbits about their parent stars, prompting much speculation as to their origin. It is known that migration alone cannot account for the observed population of these misaligned hot Jupiters, which suggests that dynamical processes after the gas disc dissipates play a substantial role in yielding the observed inclination and eccentricity distributions. One particularly promising candidate is planet-planet scattering, which is not very well understood in the non-linear regime of tides. Through three-dimensional hydrodynamical simulations of multi-orbit encounters, we show that planets that are scattered into an orbit about their parent stars with closest approach distance being less than approximately three times the tidal radius are either destroyed or completely ejected from the system. We find that as few as 5 and as many as 18 of the currently known hot Jupiters have a maximum initial apastron for scattering that lies well within the ice line, implying that these planets must have migrated either before or after the scattering event that brought them to their current positions. If stellar tides are unimportant , disk migration is required to explain the existence of the hot Jupiters present in these systems. Additionally, we find that the disruption and/or ejection of Jupiter-mass planets deposits a Sun’s worth of angular momentum onto the host star. For systems in which planet-planet scattering is common, we predict that planetary hosts have up to a 35% chance of possessing an obliquity relative to the invariable plane of greater than 90 degrees.

Subject headings:
Planet-star Interactions — Hydrodynamics — Gravitation — Chaos — Stars: rotation — Ultraviolet: planetary systems

1. Introduction

The search for planets about other stars has led to the discovery of dozens of planets with unusual properties. As both radial velocity and transit surveys are biased towards planets that are both massive and close to their parent stars, the region of parameter space corresponding to planets with a mass larger than that of Neptune and a semi-major axis au is particularly well-explored (Ida & Lin 2004; Shen & Turner 2008; Zakamska et al. 2010). This has led to the discovery of many giant planets known colloquially as hot Jupiters and Neptunes, which are thought to have formed far away from their parent stars but were then later transplanted to their observed positions by currently undetermined means. Many of these exoplanets come so close to their parent stars that they toe the line between destruction and survival, with some observed exoplanets in danger of being destroyed on a relatively short timescale (Li et al. 2010). Additionally, the inclination distribution of the hot Jupiters seems to demonstrate significant misalignment between the planet’s orbit and the stellar spin axis (Triaud et al. 2010; Schlaufman 2010), a surprising result that may require a dynamical process that acts after the protoplanetary disk dissipates.

There are three primary physical processes that can deposit a planet on an orbit that is very close to its parent star: disk migration, the Kozai mechanism, and planet-planet scattering. Disk migration can yield hot Jupiters, but as the star collapses from the same cloud as the protoplanetary disk that encircles it, it is difficult to explain the observed orbit misalignments using this mechanism alone (though see Foucart & Lai 2010; Watson et al. 2010, for discussions regarding star-disk interactions). The Kozai mechanism (Kozai 1962) can lead to the large eccentricities required to produce close-in planets, but it can only operate in systems in which a massive planet or secondary star are present, and the mechanism may be mitigated by general relativistic effects that become important before tidal dissipation is large enough to circularize the orbit (Takeda & Rasio 2005; Fabrycky & Tremaine 2007). Planet-planet scattering can produce both the observed semi-major axis and inclination distributions, and can deposit planets close enough such that tides can circularize the orbits in a time that is less than the system age. Additionally, the object that acts as a scatterer can have approximately the same mass as the scattered object itself and still yield a hot Jupiter in a significant fraction of systems (Ford & Rasio 2008), negating the need for a non-planetary companion in the system.

Previous hydrodynamical work has only focused on the planet’s first close fly-by (Faber et al. 2005, hereafter FRW), and does not investigate how prolonged tidal forcing over many orbits affects a planet’s chances for survival. In this paper we have performed hydrodynamical simulations of multiple passages of a Jupiter-like planet by a Sun-like star, bridging the gap between numerical and analytical work that have focused on extremely close and extremely grazing encounters respectively. We find that scattering planets into star-grazing orbits is more destructive than previously thought, with Jupiter-like planets being destroyed or ejected at distances no smaller than 2.7 times the tidal radius . As some exoplanets are currently observed to have semi-major axes less than twice this critical value, their initial eccentricities may be required to have been substantially smaller than unity if planet-planet scattering is the mechanism responsible for bringing them so close to their host stars. This strongly suggests that planet-planet scattering alone cannot explain the complete observed population of close-in Jupiter-like exoplanets, and that the process must operate along with one of either the Kozai mechanism, disk migration, or both. These three processes likely act in concert to produce the observed population of hot gas giants, with the relative importance of each process being a function of the system’s initial conditions.

If planet-planet scattering is common enough to explain the existence of hot Jupiters, we predict that there should be two signatures of disruption that are readily detectable with today’s instruments. Firstly, we find that the parent star can have its spin significantly altered by the accretion of material removed from the planet as a result of the disruption, producing a star that can be significantly misaligned relative to any remaining planets. Secondly, we find that most planet disruption events lead to the planet’s ejection from the host system prior to the planet being completely destroyed, and that this ejected planet can remain almost as bright as its host star for centuries.

In this paper we focus on the results of numerical hydrodynamical simulations that have been used to attempt to ascertain the true radii of destruction and ejection for Jupiter-like exoplanets, and the consequences of these planet-removing processes on their stellar hosts. In Section 2 we review the history of the analytical and numerical work done to characterize the orbital evolution of a planet that comes within a few tidal radii of its host star, and then we detail our particular numerical approach to modeling tidal disruption. We report the results of our simulations in Section 3. In Section 4 we discuss the implications of our results, with special attention paid to the viability of various mechanisms for producing hot Jupiters, and the observational signatures of planetary disruption and ejection. We summarize the shortcomings of our models and the possible fates of a Jupiter-like exoplanet in Section 5. Appendix A is provided to detail our algorithm used to simulate multiple orbits and for presenting tests of the algorithm’s conservative properties.

2. Modeling Planetary Disruption

Figure 1.— Final orbital energy scaled to the initial orbital energy (left panel), change in orbital energy scaled by (left sub-panel), and spin angular momentum scaled by the characteristic angular momentum of Jupiter (right panel) as functions of periastron distance after a single near-parabolic encounter between a Jupiter-like planet and a star. The solid lines show the results from this work, with decreasing thickness corresponding to increasing maximum refinement, square markers show the results of FRW, and the dashed lines shows the prediction of PT using FRW’s analytical solutions for and . Note that our results and that of FRW are consistent for . At the highest resolution, our results are consistent with that of PT until .
Figure 2.— The double rainbow curves show orbital trajectories for both the planet (left) and star (right) during the encounter for different values of the periastron distance . The star’s trajectory is magnified by a factor of to make its motion apparent. While the orbits do precess slightly over the course of the simulations, this plot shows the orbits with the precession removed. Note that the run (light green) experiences a particularly ejective encounter on its 2nd periastron passage.

2.1. Previous Disruption Models

Tidal encounters between a point mass and an extended, initially spherical object have been described with progressively more detailed analytical models and simulations. The first analytical models of tidal dissipation were laid out in Press & Teukolsky (1977) (hereafter PT), which assumes that the tidally excited body retains its spherical shape and that the motions induced by the encounter can be described with spherical harmonics. This model works quite well for encounters where tides are weak, i.e. more than a few tidal radii away from the perturbing body, where the assumption of sphericity is valid. If mode-mode coupling is weak, the initially excited modes have no mechanism to share energy with each other during the encounter itself, and the PT formalism can accurately predict the amount of energy injected into the extended object.

In the PT formalism, the amount of energy and angular momentum injected into an object is derived by taking the convolution of the object’s modes of oscillation and the frequency decomposition of the tidal field, which drops the time-dependence from the governing equations. However, if the incoming object is already oscillating, the phase of these oscillations relative to the time of periastron is important, which by construction cannot be resolved by the time-independent PT formalism. The absence of time-dependence also requires that the base-state is cylindrically symmetric relative to the line connecting the extended object to the perturber, which can only be accomplished for an initially non-rotating object by using a spherical geometry. As the amplitude of the mode oscillations approaches the size of the object itself, the body must become highly non-spherical and therefore the PT formalism must break down.

These two shortcomings lead to the development of the “affine” model for tidal encounters (Carter & Luminet 1983, 1985), which allows the initially spherical object to deform into a triaxial ellipsoid during and after the encounter. One axis of the ellipsoid is fixed to be perpendicular to the orbital plane, with the other two axes lying within the plane at a right angle to each other and with arbitrary orientation angle relative to the first axis. This feature allows for the object to be followed dynamically, which means that the incoming state of the object’s oscillations can affect the object’s response to future encounters.

Because the excitation of normal modes in a dynamically inert object can only yield an increase in the object’s energy and angular momentum budgets, an initially spherical object would always find itself in a more-tightly bound orbit after the encounter. An important facet of the problem that the affine model can investigate is how the fundamental modes excited in previous encounters are de-excited in future encounters, which can lead to a positive change in the orbital energy. The inclusion of de-excitation of modes adds a chaotic component to the problem, with the orbital evolution behaving as a “random-walk” process (Kochanek 1992; Mardling 1995).

While mass loss has been included in a nested-ellipsoid version of the affine model presented by Ivanov & Novikov (2001), asymmetrical mass loss is not treated as the models assume mirror symmetry. Asymmetric mass loss is expected in real disruptions as the tidal field is stronger on the side of the object facing the point mass. This is due to the steep power-law dependence of the strength of the tidal field and the non-linear evolution of the tide raised on the object, with both the velocity and amount of mass lost being larger on the side of the object closest to the point mass. Numerical simulations by FRW show that this asymmetric mass loss leads to a substantial deviation from the nested-ellipsoid treatment, especially when the closest approach distance is . Within this distance, the asymmetric removal of mass can lead to a positive increase in the orbital energy, which can cause the object to be completely ejected from the system if .

The disruption of polytropes in highly eccentric encounters has been investigated numerically by both Lagrangian (Nolthenius & Katz 1982; Bicknell & Gingold 1983; Evans & Kochanek 1989; Kobayashi et al. 2004; Faber et al. 2005; Rosswog et al. 2008, 2009; Ramirez-Ruiz & Rosswog 2009; Lodato et al. 2009) and Eulerian methods (Khokhlov et al. 1993a, b; Frolov et al. 1994; Diener et al. 1997; Guillochon et al. 2009), with the principle focus being on stars or compact objects that are disrupted by point-like gravitational sources. However, most hydrodynamic simulations that have focused on the long-term survival of these systems only consider when the two objects have a mass ratio close to unity (Lee et al. 2010; Lorén-Aguilar et al. 2009), or for non-disruptive encounters (Rathore 2005). Recently, Antonini et al. (2010) performed low-resolution multiple-passage simulations in the context of the galactic center, allowing the exploration of a large parameter space at the expense of accuracy.

To summarize, most analytical models of tides in an astrophysical context focus on objects which do not lose any mass at their closest approach, and conversely most numerical work has focused on encounters where the object is completely destroyed. The intermediate regime, where objects lose some mass but are not completely destroyed in their first passage, is largely uncharacterized. If the planet survives the initial encounter, some of the mass that is removed from the planet can become bound to the planet again as it recedes from periastron to a region with a weaker tidal field. The return and subsequent re-accretion of this material is not treated at all in analytical models. As we will describe in Section 3.2, the structure of the planet can be significantly altered by the re-accretion of the planetary envelope, and the inclusion of this re-accretion into tidal disruption theory is necessary to determine if a planet will ultimately survive.

2.2. Our Approach

While Lagrangian codes are well-suited for treating problems where a hydrostatic object moves rapidly with respect to a fixed reference frame, their relatively poor scalability makes it difficult to follow the evolution of an object for many dynamical timescales with sufficient spatial resolution. On the other hand, maintaining near-hydrostatic balance in rapidly advecting frames using Eulerian methods can also be quite challenging. Robertson et al. (2009) show that performing Eulerian simulations in a boosted frame with unresolved pressure gradients leads to a spurious viscosity term that tends to damp out instabilities. This has been colloquially referred to as the non-Galilean invariance (GI) of the Riemann problem (Springel 2010; Tasker et al. 2008). While Robertson et al. showed that GI issues can be avoided with increased resolution, the requisite spatial resolution can be impractical for three-dimensional simulations that are evolved over thousands of dynamical timescales, or for when the bulk velocities are many times larger than the internal sound speeds within the simulation. The problem is compounded in the outermost layers of Jupiter-like planets, which exhibit particularly steep pressure gradients.

To side-step the GI issues, our simulations are performed in the rest-frame of the planet, which limits the typical fluid velocities to the planet’s sound speed for tidal disruption calculations (Guillochon et al. 2009). We model the planet as an , polytrope, with the fluid being described by a polytropic equation of state , where the adiabatic index . This gives a reasonable approximation to the interior structure of Jupiter-like planets (Hubbard 1984), as long as the core is not a significant fraction of the planet’s mass. Our simulations are constructed within the framework of FLASH (Fryxell et al. 2000), an adaptive-mesh, grid-based hydrodynamics code. We treat the star as a point-mass because the distortion of the star itself contributes negligibly to the planet’s orbital evolution in the case where (Matsumura et al. 2008). The relative positions of the star and the planet are explicitly tracked over the course of the simulation using two virtual particles (star) and (planet), which are evolved alongside the hydro calculation using a Burlisch-Stoer integrator (Press et al. 1986) with maximum error constrained to be less than . The planet’s self-gravity is calculated using a multipole expansion of the fluid, where we take .

The flux-calculation step in any hydro code depends on the time-dependent force applied to the fluid over the duration of the step. Because the gravitational field is applied as a source term in most multi-dimensional hydro codes that include self-gravity, conservation of energy is not explicitly achievable when self-gravity is included, as the potential at timestep is unknown and must be estimated. In the case of self-bound objects in hydrostatic equilibrium, small perturbations can lead to a systematic drift in the total energy of the system. For hydrostatic objects, we find that an increase in spatial resolution reduces this drift by an amount that approximately scales with the number of voxels used to resolve a given region. This is because the relative forces acting on neighboring cells decrease as the voxels are more closely packed together, and because the sharp density gradients present in the 1D hydrostatic profile are better resolved.

The potential at must be estimated because the value of depends on , which is not known until the hydro step has been completed. This extrapolation is one of the main sources of error in the code because the obtained acceleration is only an estimate based on the previous evolution of . And because the extrapolation is only first-order accurate in time, smaller time-steps do not improve the accuracy of the results (Springel 2010). The error can be somewhat reduced by using a higher-order extrapolation that includes previous time-steps, but our tests show that including the time-step only leads to a few percent reduction in error relative to the first-order approximation. Additionally, it imposes additional memory and disk overhead to save the full potential field from the previous time-steps.

For cases where the object is nearly in hydrostatic equilibrium the gravitational field does not rapidly change with time. However, a tidal disruption can lead to variations in the gravitational field of order unity over a fraction of the planet’s dynamical timescale. As a result, the default method used by FLASH to apply the gravitational field can lead to substantial changes in the total energy. Motivated by this, we implement a novel method where the potential contribution from the fluid and from tidal forces are separated, negating much of the error associated with the extrapolation of the potential. Details of this algorithm are provided in Appendix A.

3. Simulation Results

3.1. Single Passage Encounters

Figure 3.— Slices through the orbital plane shortly after each periastron passage for the simulation where . All plots show log . The upper, red color-coded figures show a wide view of each encounter, with white corresponding to g cm and black corresponding to g cm, while the lower, blue color-coded figures show a close-up view of the core, with white corresponding to the maximum density and black corresponding to g cm.
Figure 4.— Same as Figure 3, but with frames 2-10 corresponding to apastron. The first frame shows the initial conditions.

For comparison purposes, we first ran a suite of simulations with physical initial conditions identical to that of FRW. The planet is assumed to have a radius and mass , where and are the radius and mass of Jupiter. The planets are disrupted by a star with , with the orbits of incoming planets are set to have an apastron separation . Our lowest-resolution models have maximum spatial resolutions of , where is the width of the smallest grid cells, corresponding to , slightly better than the peak spatial resolution of FRW. Our results agree very well with FRW’s results for periastron passage distances (Figure 1).

Because the amount of energy stored in the oscillations is where is the binding energy of Jupiter and is the amplitude of the mode, the simulations converge to the analytical solution only when the spatial resolution is fine enough to resolve . As is evident in both FRW and our lower-resolution runs, the measured amount of energy dissipation is larger than the analytical predictions of PT when the oscillatory amplitude is smaller than the minimum grid scale. To test for convergent behavior, we ran higher-resolution simulations with double and quadrupole the resolution of our fiducial test for the more grazing encounters. For the and runs, the improved spatial resolution allows us to recover the analytical solution. It is also apparent that our simulation is closer to convergence than the lower-resolution models, but the estimated resolution required for true convergence () means that recovering the predicted results of PT for such a grazing passage is currently beyond our computational ability.

As in FRW and Khokhlov et al. (1993a), we also find that the change in and is slightly larger than what is predicted by the analytical models, even in the simulations that have surely converged and have minimal mass loss (i.e., ). This is almost certainly due to the dynamical tide effects neglected by the PT model. Comparing the hydrodynamic results to that of a dynamical treatment of tides shows better agreement for this range of pericenter distances (Lai et al. 1994), but note that these dynamical models only include the -mode of oscillation, and thus do not account for energy transferred to higher-order -modes or -modes. For grazing encounters, the dynamical tide is much less important, and thus we expect that the change in orbital energy and angular momentum should converge to the PT prediction for an inviscid planet.

3.2. Multiple Passage Encounters

Figure 5.— Fallback accretion stream formed as the result of the 8 encounter between a 1 star and 1 planet with initial . The large, red color-coded image shows a wide-view of the disruption, while the inset blue color-coded figure shows a close-up of the surviving planetary core.

The initial conditions of FRW are not appropriate for investigating multiple-passage encounters as the orbital timescale is many thousands of times longer than the dynamical timescale. Thus, for our second set of simulations in which we explore multiple passages, we assume that the planets are initially on orbits (Figure 2). Depending on the initial pericenter distance, the planets are destroyed in as few as one and as many as ten orbits (Figures 3 and 4), where a planet is considered to be destroyed once it has lost more than 90% of its original mass.

To ensure that our assumption of a less eccentric orbit is valid, we consider the effects of varying when . For nearly-parabolic orbits, the shape of the orbit in the vicinity of periastron changes very little with changing , with the main effect being that the average strength of the tidal force is slightly stronger for smaller before and after pericenter. This means that the critical distance for which a planet is destroyed or ejected is very slightly larger than our setup for single encounters described in the previous section where . To estimate the magnitude of this effect, we calculate the ratio of the distances at fixed true anomaly for two orbits with eccentricities and , assuming both orbits have the same


As the strength of the tidal force is , the difference in the strength of the tidal force between the two orbits with respect to the tidal force experienced at pericenter is


This expression is maximized at , where the force difference between an and encounter evaluates to . For all other values of , the force differences are much smaller than this maximum. Thus, our hydrodynamical treatment of tides in a orbit is directly applicable to all orbits with .

After an orbit in which the planet sheds mass, some of the material that is removed from the planet’s surface remains marginally bound to the planetary core. The majority of this material is then re-accreted by the planet over a few dynamical timescales. When the free-falling material encounters the surviving planetary remnant, its kinetic energy is converted to internal energy in a standing accretion shock, which results in the planet possessing a hot outer layer with temperature close to the virial temperature (Figure 5). Additionally, the material striking the remnant leads to some heating of the remnant’s outer mass shells. This effect is predominantly important in determining the envelope’s temperature early in the accretion history before the pressure of the growing hot atmosphere becomes comparable to the ram pressure and is able to halt the flow before it can reach the core’s surface (Frank et al. 2002).

Figure 6.— Pressure perturbations and spin excited by a grazing encounter with . The upper five panels show the pressure perturbation over several , with blue corresponding to regions of lower pressure and pink corresponding to regions of higher pressure relative to the base state. The lower five panels show the fluid’s angular momentum relative to , with green representing clockwise rotation and orange representing counter-clockwise. Note that while the object appears to be rotating rapidly, the presence of rotating and counter-rotating regions leads to almost an exact cancellation of the total angular momentum . The illusion of rapid rotation is in fact related to the pattern speed of the , normal modes. Because the angular momentum carried by this mode is related to the tangential component of the expansion and contraction of the planet as it oscillates, the fluid vacillates back and forth at the mode frequency. As can be seen in the figure, the fluid with possess slightly more total angular momentum than the fluid with . This leads to a spin frequency that is actually very small.

Because the re-accreted material is marginally bound to planet, the orbital trajectories characterizing the accretion streams have apocenter distances comparable to the planet’s Hill sphere and follow highly eccentric orbits. As the re-accreted material returns on a Keplerian trajectory (Kochanek 1994; Ramirez-Ruiz & Rosswog 2009), it carries a substantial amount of specific angular momentum. This leads to a rapid spin-up of the planetary remnant’s outer layers. In encounters with little or no mass loss, the planets spin slowly post-encounter as the angular momentum carried by the normal modes is almost equally distributed between rotating and counter-rotating regions (Figure 6). The process of re-accretion produces a planet that has a lower average density and more mass at larger radii, which makes the planet easier to destroy on subsequent passages.

As the re-accreted material has a temperature comparable to the virial temperature, radiative cooling may be able affect the atmosphere’s structure, an effect that we do not account for in our simulations. However, as the planet can only thermally evolve for one orbital period, the atmosphere does not have much time to cool down before it has another strong tidal encounter with the star. While the outer layers of the planet may cool relatively rapidly and lead produce a brief transient visible in the UV (see Section 4.4), the denser regions of the planet’s hot atmosphere component contains much more mass and is too optically thick to cool significantly before returning to pericenter. Assuming Thomson scattering, this corresponds to a mass of the hot atmosphere component . This means that we expect that the thermal evolution of the reaccreted material is unimportant in determining the planet’s density profile interior to the most tenuous outer layers. Thus, we expect that the dynamically relevant distribution of mass within the planet should remain relatively unchanged between encounters with the star, even for values of significantly closer to 1 than what we use in our simulations.

Figure 7.— Change in orbital energy attributed to each passage as a function of and orbit number . The orange regions show decreases in (more bound), while the cyan regions show increases in (less bound). The height of each column shows the number of orbits a planet survives before being destroyed.

This allows us to use the change in orbital energy and angular momentum from our simulations to calculate the orbital evolution for orbits with semi-major axes of several au. The change in orbital energy as a result of each encounter is shown in Figure 7. For grazing encounters with little mass loss, the change in orbital energy is negative; the planet becomes progressively more bound to the star after each encounter. The magnitude of this change is related to the state of the dynamical tide on the planet at pericenter, where the interaction between the oscillation of the fundamental modes can interact with the tidal field to reduce the amount of energy removed from the orbit, or even change its sign (see Section 3.3). As the mass loss per orbit exceeds , the trend in orbital energy change becomes positive, and planets become progressively less bound on each subsequent encounter. This is primarily a result of the asymmetrical mass loss, although the interaction with the normal modes is important for encounters where the mass loss is on the order of a few percent.

For orbits with , where au (Ida & Lin 2008) is the ice line, the total orbital energy is small relative to the self-binding energy of the planet, and a positive can lead to the surviving planet becoming unbound from the host star. The number of orbits completed by a planet before it becomes unbound is shown in Figure 8. All Jupiter-like planets that come within are expected to be ejected if , and destroyed otherwise. For the most grazing encounters and smallest eccentricities, we find that the planet can survive for as many as ten orbits before being destroyed by the host star. The general trend is that planets that come closer to their parent stars tend to be ejected or destroyed in fewer orbits, but the consequences of the encounter seem to depend heavily on the orientation of the planet’s long axis at pericenter, which as we will show in the next section is quite difficult to predict.

Figure 8.— Criteria leading to planet ejection given an initial periastron passage distance and eccentricity . The colored regions correspond to the number of orbits before a planet is ejected by its interaction with the parent star. The dashed lines show the value of given an apocenter distance equal to the ice line , assuming that . For all values of shown and for , planets are ejected before they are totally disrupted. Note that the region bounding the planets that are ejected on the first orbit has a monotonic dependence on , while all other regions exhibit more complicated structure. This is because the first passage involves a planet with no internal motions, and thus no relation between phase and the excitation or de-excitation of fundamental modes. All future passages for involve a planet that is both differentially rotating and oscillating, resulting in a large variance in the amount of energy added or removed from the orbit.

3.3. The Role of Chaos

As the planet continues in its orbit after its first encounter with the parent star, the planet exhibits both rotation and oscillation. Both the magnitude and sign of the orbital energy change is related to the phase of the planet’s dynamical tide at periastron. Because the tidal forces strongly excite the modes, the coherence or decoherence of these modes at periastron can greatly change the effects of that particular encounter (Kochanek 1992; Mardling 1995; Usami & Fujimoto 1997). As the ratio of the orbital period to the break-up rotation period for a Jupiter-like planet scattered in from the ice line is , even a very small change in the orbital parameters introduces a dramatic variability in the amount mass lost and the change in orbital energy.

For our simulations where , the planet completes hundreds revolutions between periastron passages, which means extremely fine sampling of would be required to completely describe the problem. To explore the chaotic behavior we ran another multiple disruption simulation with , with the only difference being that the initial eccentricity is set to , corresponding to an orbital period that is one free-fall timescale longer than the corresponding simulation. As shown in Figure 9, multiple outcomes are possible for even slight changes in the initial conditions, with the planet surviving for a different number of orbits in the two simulations. The chaotic behavior is also evident for our multiple disruption simulation where , which is destroyed on its fourth periastron passage, whereas both the and runs are destroyed on their third passages. This chaotic behavior arises because the angle of the major axis of the oscillating planetary core can range from being perpendicular to parallel to the true anomaly at periastron, which dramatically affects the ability of a planet to survive on any particular orbit.

Figure 9.— The evolution of two multiple encounter simulations with almost identical initial conditions. Both simulations use the same initial conditions as the rest of the multiple encounter simulations presented in this paper, with . The only difference between the two simulations is the initial eccentricity: The solid curves show the outcome for initial eccentricity , while the dotted curves show the outcome for , which corresponds to an orbital period that is one free-fall time-scale longer than the case. The four plots show the planetary mass , orbital energy , spin angular momentum , and precession of the orbit as functions of . Note that while the outcome of the first encounter at is almost completely identical in both simulations, the behavior diverges on the second encounter, which leads to the planet on the orbit being destroyed in three orbits, whereas the planet on the is destroyed in four. This divergence is a result of the phase difference between the two simulations introduced by the slight difference in initial orbital period.

As the principle component of the tidal field are the harmonics, the only way to avoid chaotic behavior is to somehow remove energy from these modes before the planet returns to periastron. For fully-convective, Jupiter-like planets, this can be potentially achieved through three types of mechanisms: Viscosity (either microscopic or turbulent), coupling to other oscillatory modes, or the increase in entropy associated with sound waves steepening into shocks at either the surface of the planet or within its interior. It should be emphasized that most of the work investigating these energy-sharing mechanisms have concentrated on systems where the oscillations can be treated linearly, and thus the results of these studies can only give us a rough idea to their importance when applied to partially-disruptive encounters.

The two main differences between the linear models and the survivors of a partially-disruptive encounter are the amplitude of the oscillations, which are close to unity, and the presence of a hot, optically thick envelope that accumulates after the disruption and sits on top of the oscillating core. We know that in order for a particular mechanism to result in significant damping of the modes, the damping timescale has to be at least on the order of the orbital period , if not shorter. As all of the mechanisms for removal of energy from the modes result in a cascade to microscopic scales, systems with effective damping will experience inflation of the core, inflation of the envelope, or inflation of both regions. This inevitably leads to reduced survivability on subsequent passages.

We now briefly discuss the applicability and viability of each of these proposed mechanisms. For Jupiter-like planets, the microscopic and turbulent fluid viscosities seem to be too small to produce any significant damping on an orbital timescale (Guillot et al. 2004). Perhaps a more promising mechanism is the coupling of the primary modes to higher-order “daughter” modes, which then couple to “grand-daughter” modes, etc., in a cascade resembling the cascade of energy from large scales to small scales in turbulent fluids (Kumar & Goodman 1996). In the linear regime, the fundamental mode normally couples to the low frequency g-modes, with the degree of coupling being related to the amplitude of the primary perturbation relative to the size of the object, . But for Jupiter-like planets, which are fully convective and have a negligible luminosity, the polytropic index is equal to the adiabatic index , which makes their interiors incapable of supporting g-modes (Cowling 1941). Coupling can still occur through p-modes, which have higher frequencies than the fundamental mode, but the coupling is only effective for large displacements where the behavior becomes non-linear and for which the rate of energy-sharing is highly uncertain (Kumar & Goodman 1996).

Inertial waves may be an effective means of dissipating the modes given low-frequency tidal forcing (Ivanov & Papaloizou 2010), but the efficiency of this dissipation as the forcing frequency approaches the characteristic frequency is not well understood. And because inertial waves are most effective when the planet spin frequency is comparable to the orbital frequency, they may not be important during the first few passages before syncronicity is established. Additional dissipation may occur via the interaction between inertial waves and Hough waves (Ogilvie & Lin 2004), in which effective coupling is achieved when the wavelength of the inertial modes is comparable to the size of the radiative zone. As scattered planets tend to have large eccentricities, stellar insolation is unlikely to be important for these planets prior to circularization, but as the hot envelope produced by partially disruptive events is radiative, effective dissipation via Hough waves may still be possible.

In our simulations, we do not find that the modes decay appreciably between periastron passages, despite the fact that the planet oscillates thousands of times in the course of each orbit. We do find, however, that the re-accretion of loosely-bound material moderates the chaotic behavior somewhat. As the density in this region is times smaller than the core, the tidal radius for the hot envelope component is larger than the core’s initial tidal radius, which results in the hot envelope being easily removed on subsequent encounters. And unlike the core, the hot envelope has a large sound-crossing time relative to the periastron passage time, and thus no fundamental modes that could affect the interaction on future passages are excited in this region. This leads to the result that as the envelope becomes a larger fraction of the planet’s total mass, the behavior becomes notably less chaotic.

3.4. Debris Accreted by the Star

Figure 10.— Mass loss history of multiple passage encounters for different initial values of . Each curve is color-coded to correspond to a particular orbit number. The solid lines show the aggregate mass accreted by the star, while the dashed lines show the total mass lost from the planet.

The amount of mass removed on the first few orbits can vary over many orders of magnitude for a relatively small range of . For our closest simulated encounters, nearly half the planet ends up being bound to the host star after the first passage, while our most grazing encounters show almost no mass loss at all (Figure 10). However, as the initial passages excite fundamental modes of oscillation, and the returning debris leads to additional heating in the planet’s outer layers, the amount of mass removed geometrically increases with each subsequent passage. This leads to the result that approximately 20-50% of the planet’s initial mass is accreted by the host star by the time the planet has been completely disrupted, with slightly less material being accreted when the initial encounter is grazing.

The reason closer disruptions lead to more mass accreted by the central star has to do with the asymmetry of the tides; the force resulting from the inner tide (closest to the star) at pericenter is


times larger than the outer tide (furthest from the star). And while the effective size of the planet increases due to the excitation of oscillations and the loss of mass, subsequent passages yield diminishing returns as material has already been removed on previous encounters.

As the mass lost by the planet from a series of partially disruptive encounters is highly stochastic, the exact amount of mass accreted by the star for any given multiple-orbit encounter is difficult to predict. To parameterize the average mass accreted by the star, we consider two limiting cases. If the planet is scattered into an orbit with , the planet will be completely destroyed before it is ejected. On the other hand, if a planet is scattered from the ice line, is much smaller than the planet’s binding energy and the planet will be ejected on the orbit for which . A fit of for these two limiting cases yields


for , where the impact parameter . For the case where the planet is completely destroyed, the amount of mass accreted is nearly constant, with a slight decrease with decreasing . The decrease is substantially steeper for the incomplete disruptions, as planets on grazing orbits are less likely to transfer much mass to their parent stars prior to being ejected.

4. Discussion

4.1. The Jupiter Exclusion Zone

Figure 11.— Possible apastrons for the known hot Jupiters with . Each arc shows the apocenter of an orbit with pericenter , assuming that the angular momentum of the initial orbit is equal to the orbital angular momentum observed today. The open circles show each planet’s currently observed and , scaled to the tidal radius and the ice line respectively, which are determined by the planet/host star properties. Blue-filled circles are planets that are thought to be aligned with their host stars, red-filled circles are thought to be misaligned (either via direct measurement or because they have observed to have statistically inconsistent rotation rates, see Schlaufman 2010), and white-filled circles show systems with unknown orientations. The thick vertical line (labeled ) shows the minimum possible value of corresponding to the Jupiter exclusion zone, which is determined through our numerical simulations for Jupiter-like planets to be , while the thinner vertical line (labeled ) shows the maximum possible value for , which is defined by the planet that is currently observed to be closest to its classically defined tidal radius, WASP-19 b. Filled black circles show the intersection between the lines and each of the arcs of constant angular momentum shows the maximum apastron distance a planet could have been scattered from without being destroyed on its subsequent encounters with the host star. Note that the typical values are significantly smaller than , which indicates that those planets must have migrated prior to being scattered if planet-planet scattering brought them to their current positions. All data taken from the Extrasolar Planets Encyclopaedia (http://exoplanet.eu) and René Heller’s Holt-Rossiter-McLaughlin Encyclopaedia (http://www.aip.de/People/RHeller/) on December 8, 2010.

As discussed in Section 3, there exists an exclusion zone with radius within which all Jupiter-like planets are either ejected or destroyed. For our simulations we used a polytropic model of a gas giant, with a mass and radius equal to present-day Jupiter. As some gas giant planets may not have had much time to cool before being disrupted, it is likely that disrupted Jupiter-like planets are larger than our fiducial cold Jupiter model (Bodenheimer et al. 2001). This is confirmed by observations, many of the known hot Jupiters are observed to be inflated, with larger radii and lower densities as a result of the injection of entropy via stellar insolation (Fortney et al. 2007; Miller et al. 2009; Li et al. 2010), and potentially also other dissipative mechanisms (Ogilvie & Lin 2004; Laine et al. 2009; Arras & Socrates 2010). But this means that the tidal radii for these planets must be larger than what we use in our calculations, translating to planets that are more easily disrupted than cold gas giants. And while the core mass of Jupiter-like planets is uncertain, all plausible models include cores that compose such a small fraction of the planet’s total mass that they are unimportant in determining the structure of the gas envelope, and thus the dynamics of the disruption. The value of calculated in this paper is thus a lower limit for Jupiter-like planets. For gas giants with masses more similar to Neptune, the affect of an increased core mass may allow these planets to survive closer to their parent stars due to the increase in average density. We stress that our model should only be applied to planets of .

Because the spatial resolution required to resolve grazing passages becomes computationally prohibitive beyond , the true value of may lie beyond what we have been able to calculate for even cold Jupiter-like planets. Thus, the vertical cut-off line denoted by in Figure 11 almost certainly lies closer to the right of the diagram, which would add more planets to the list that could not have been directly scattered to their present positions from the ice line. However, the true cutoff value certainly lies no further than the current pericenter distance of WASP-19 b (denoted as ), currently the exoplanet with the smallest known value of , which appears to be quite inflated (Hebb et al. 2010) and probably has a tidal radius that is larger than a cold gas giant of the same mass. Hence, we can only constrain to lie within the range , or .

4.2. Implications for the Formation of Hot Jupiters

Planet 5 6
HAT-P-12 b 5.2 0.41  —  0.4 — 4
OGLE-TR-56 b 4.8 0.029  —  0.9 — 2
WASP-4 b 5.1 0.093  —  1 — 6
WASP-12 b 4.7 0.021  —  0.5 — 1
WASP-19 b 3.5 0.010  —  0.07 — 1

Because specific orbital angular momentum is very nearly conserved during even strong tidal encounters (Figure 2), the currently observed semi-major axes of hot Jupiters is at most 2 , where is the planet’s initial pericenter distance after being scattered into a disruptive orbit. As the exclusion radius is larger than the for many known exoplanets, there exists a maximum initial apastron distance from which a planet could have been scattered without having been destroyed or ejected


where and are the currently observed semi-major axis and eccentricity. This maximum apastron value is illustrated for the currently known hot Jupiters in Figure 11. Because a planet that is scattered into a disruptive orbit from the ice line possesses , the maximum allowed apastron is a highly sensitive function of for . Thus, if the current distance of an exoplanet from its host star is less than , it is highly likely that its initial orbit prior to scattering must have been significantly closer to the star than the ice line. This means that no Jupiter-like planets would be observed for initial pericenter distances smaller than this value if they were scattered from the ice line to their present-day orbits.

Five of the currently known hot Jupiters (HAT-P-12 b, OGLE-TR-56 b, WASP-4 b, WASP-12 b, and WASP-19 b) have observed semi-major axes , corresponding to (Table 4.2). For four of the five planets, the maximum initial apastron distance is less than 0.1 . As all Jupiter-mass planets are thought to form beyond , these planets could not have been directly scattered to their current locations. This implies that a migration process that resulted in a dramatic decrease in must have occurred either before or after these planets were scattered. If the true value of for cold Jupiter-like planets lies closer to the maximum possible value set by WASP-19 b (), as many as 18 of the currently known hot Jupiters with could not have survived the scattering event that deposited them at their current locations.

As both planet-planet scattering and the Kozai mechanism do not radically alter a planet’s semi-major axis, disk migration would be required to explain how close-in planets would have initial apastrons that are significantly smaller than (Lin et al. 1996; Ida & Lin 2004, 2008). As many of the currently known exoplanets have , this would imply that the migration timescale in these systems must have been substantially shorter than the disk lifetime if the planets migrated prior to the scattering event. For the misaligned systems, a dynamical process (either planet-planet scattering or the Kozai mechanism) may need to occur after the migration or while the gas disk dissipates (Matsumura et al. 2010) to produce the observed misalignment.

Alternatively, the planets may have migrated after the scattering event. For this to be the case, the scattering event had to bring the planet close enough to its parent star such that the tide raised on the star can alter the planet’s orbit in a time shorter than the system age, but not too close that the planet is destroyed or ejected by the star. Suppose that a planet is scattered from beyond such that . After circularization the pericenter distance doubles to , and the planet’s semi-major axis evolves via the interaction between the planet and the tide it raises on the surface of the star (Eggleton et al. 1998; Hut 1980). At this distance, the orbital period of the planet is almost always shorter than the spin-period of the star, except perhaps for stars with ages Myr (Irwin & Bouvier 2009). This results in a spin-up of the star, which tries to “catch up” to the Keplerian frequency imposed by the planet’s orbit, and an inward migration of the planet.

As the planet is likely to be tidally locked even prior to circularization, the timescale for evolution of the planet’s semi-major axis is entirely determined by the star’s properties and the orbital frequency (Dobbs-Dixon et al. 2004),


where is the star’s radius, is the star’s tidal quality factor, and is its rotation frequency. The fastest inward migration occurs when a star is not rotating, e.g. . As planet-planet scattering seems to be rare after yr (Matsumura et al. 2008), we can set equal to the system age to determine , the maximum tidal quality factor for which a planet can migrate from to


where is the initial orbital period. When setting , all of the known hot Jupiters with yields values for that are consistent with those expected for stars (Table 4.2).

If the true radius of disruption lies closer to , the measured values of are more restrictive, meaning that higher dissipation rates would be required to enable planet-planet scattering to viably produce hot Jupiters with . This would also imply that most observed hot Jupiters would only exist a short while longer before being destroyed by their parent star. This remaining lifetime can be calculated by solving Equation (7) using present-day values for , , and (denoted in Table 4.2). As is extremely sensitive to the exact value of , only a slight increase in is required to eliminate direct scattering as a possible genesis mechanism for a significant fraction of the hot Jupiters. If it can be shown that is substantially larger than what we have calculated in this work (i.e. ), can put definitive constraints on the mechanism responsible for producing hot Jupiters.

4.3. Stellar Spin-Up from Planetary Disruption

Figure 12.— Changes to the stellar spin as a result of accreting tidally disrupted planets. The curves show the cumulative probability that a star will possess a given angular momentum and spin inclination after 1 (solid) or 5 (dashed) planetary disruptions for different initial stellar angular momentum . Planets are assumed to have a logarithmic distribution in mass () and semi-major axes , where is the minimum distance for which a planet won’t be tidally disrupted and is the ice line. The eccentricity and inclination relative to the invariable plane are assumed to follow Rayleigh distributions, with and .

As discussed in Section 3.4, the star acquires a substantial injection of angular momentum from the disrupted planet. Our simulations give us an empirical determination of , the amount of mass gained by a star with from the disruption of a Jupiter-mass planet given the impact parameter (Figure 10). The angular momentum acquired by the star is simply equal to the specific angular momentum at the star’s radius multiplied by the amount of mass accreted, . Note that this is not the same as the total angular momentum content of the disk formed from the debris as claimed by Jackson et al. (2009), as the accretion stream does not intersect the star’s radius unless the planet’s original orbit has (Kochanek 1994).

If all disrupted Jupiter-like planets are approximately the same size and , all disruptions can be treated using our fiducial model and the inclusion of the simple scaling associated with the increase in the orbital angular momentum. This is because the only difference that arises from changing the mass of the planet and the star is the degree of asymmetry of the tides, which is third-order in the force expansion and is times smaller than the tidal force itself. Thus, in the event that the planet is completely destroyed by the star, we expect that the amount of angular momentum acquired by the star should simply scale with and


However, our results show that a substantial fraction of planets are ejected from the system before they can be fully destroyed by the star. This reduces the amount of mass accreted by the star, and thus . Therefore for any given encounter should depend on which particular orbit a planet is ejected . This parameter depends on the orbital energy , and the change in orbital energy associated with each encounter, which as we explain in Section 3.2 should also should scale simply with the masses of the planet and the star. Thus our expression for becomes


where is equal to the lowest value of for which Equation (11) is satisfied. While this gives us the amount of angular momentum acquired by the star for a set of encounter parameters, we must know the distribution of in order to estimate the probability for which a star will possess a total angular momentum and obliquity after a number of planetary disruptions .

The integrated rate of scattering (i.e., all encounters with greater than some value) has been evaluated numerically by a number of authors. Jurić & Tremaine (2008) show that up to 20% of planets present at the end of phase I of planetary formation can collide with the host star, e.g. . Ford & Rasio (2008) show that up to 16% of planets in a 3-body system can be thrown into a “star-grazing” orbit, which they define as . Nagasawa et al. (2008) show that planet-planet scattering events tend to induce transitions between Kozai states (with the duration of each state being yr) until a planet is ejected from the system or until the eccentricity of the innermost planet is damped by tides. Nagasawa et al. note that most of the close-in planets seem to be driven to their closest approaches via the Kozai mechanism, which gently drives the planets into the region where they can be circularized, as opposed to being directly scattered into such orbits. None of the models presented above include the precession associated with general relativity that would normally quench the Kozai mechanism.

Because the numerical scattering experiments do not include a hydrodynamical treatment of tides, the fates of planets that are either scattered or driven by the Kozai to are not accurately represented. What the scattering experiments do reveal is the total rate of planet-planet interactions as a function of the number of gravitating bodies in the system, and the distribution of orbits that arises from numerous planet-planet interactions. And despite the simplistic treatment of tides or the neglect of tides altogether, different prescriptions of tidal dissipation do not seem to strongly affect the distribution of the remaining planets in the system (Nagasawa et al. 2008), which means that the orbit distributions these models predict should still be appropriate to use as inputs for our post-disruption stellar spin estimates.

For systems that are not dynamically stable after the gas disk dissipates, or for systems which are driven to instability by an external perturber (Zhou et al. 2007; Malmberg et al. 2010), the models indicate that the eccentricity distribution of planets quickly evolves to a Rayleigh distribution (Jurić & Tremaine 2008)


As the relaxation to this distribution is mainly driven by strong two-body planet-planet interactions, the new eccentricity of a planet after a scattering event should also be drawn from the above distribution. Because we are interested in objects that may be unbound from the system after the encounter, objects with initially larger than 1 should be included when calculating the number of events at each . Using our definition of ,


and by making a change of variable from to , Equation (12) becomes


When , as is the case in planet disruption, the above expression simplifies to


which when integrated yields a disruption probability that scales inversely with (Rees 1988). Because we are including hyperbolic encounters, can assume both positive and negative values, with . To first order, the value of the integrand is equal for both positive and negative values of , and the total number of events where is given by


This implies that equal numbers of planets will be scattered into prograde and retrograde orbits. Under these assumptions and given our empirically determined lower limit for ejection , it is immediately clear that the rate of collisions with the central star is lower than the combined rate of ejections and disruptions by at least a factor , assuming solar and Jovian values.

Now that we have a model for the expected initial distribution of giant planets in a dynamically relaxed system, we can use Equation (14) to evaluate the expected values of and after planetary disruptions for a star of (Figure 12). Here we consider the scattered objects to be cold, Jupiter-like planets with , resulting in a nearly-constant planetary radius . This assumption is only valid if the system is older than yr and if the planet’s initial orbit is far enough from its parent star to have negligible insolation, but note that including these effects would only act to increase the amount of mass removed from the planet per orbit. We also assume that the planets are distributed uniformly in log from to , and in log from 0.1 to 10 , with Rayleigh distributions in and with and . For simplicity, we assume that any angular momentum acquired by the star through a disruption is shared equally with all parts of the star.

For a single planet disruption in a system where the star possesses initial angular momentum equal to the Sun, exceeds 30 in 15% of the stellar hosts, and 90 (i.e. the star rotates retrograde to the invariable plane) in 8% of systems. The spin rate of the star also tends to increase, with 10% of stars possessing after the disruption. If disruptions are very common (), the probability of increases to 47/22%, and 41% of stars have . The probability for enhanced values of and are all slightly smaller if we restrict disruptions to come from , as the amount of mass the star acquires from disruptions is slightly less on average, Equation (5). However, the effect is minor, with changes in the cumulative probabilities being on the order of a few percent.

If a star is unable to share the angular momentum deposited in its outer layers in a time less than its age, the star may be observed to have larger values of or than what would be expected given complete mixing. The timescale for sharing angular momentum across the tachocline in the Sun is known to be only 3 Myr (Gough & McIntyre 1998), with the timescale decreasing for increasing rotation rates in the convective region. The timescale for sharing angular momentum across the tachocline appears to increase significantly as the size of the convective zone shrinks for rapidly rotating stars (Barnes 2003; Barnes & Kim 2010), but this may be moderated by a fingering instability made possible by the larger molecular weight of the disruption debris (Garaud 2010; Rosenblum et al. 2010). A disruption in a system where the host star never has a thick convective zone could effectively erase the star’s original spin rate and obliquity, with being pulled from the inclination distribution of planetary orbits.

4.4. Observational Signatures

If we assume that the pressure profile of the hot atmosphere that accumulates on the planet has a pressure profile at all radii (Frank et al. 2002), the initial Kelvin-Helmholtz cooling timescale given a total atmosphere mass is approximately


where is the virial temperature, is the mass of the atmosphere and is the radius of the atmosphere. Because is initially very large, the atmosphere is at first completely ionized. At this temperature, is a few days, leading to a rapid thermal evolution of the planet’s outer layers shortly after most of the sundered mass returns to the planet. During this phase, the planet can briefly outshine its parent star with erg s, with most of the radiation being emitted in the UV. As the atmosphere cools it shrinks back down onto the planet’s surface until . We can then estimate the temperature at which the planet radiates for the majority of its orbit by setting equal to , which yields a temperature of a few K. This indicates that the planet’s outer layers will still be somewhat inflated before the planet returns to pericenter, and thus will be easily removed on subsequent passages.

For planets that are ejected from their host stars, the thermal evolution of their outer layers continues until the temperature reaches a few thousand Kelvin, at which point hydrogen begins to recombine, which acts as a thermostat to maintain the temperature at a relatively constant value. The recombination timescale is


where is the Hydrogen fraction. As a result of this relatively long time-scale, these ejected planets could remain quite bright () for an extended period of time, even without additional tidal forcing. If we assume that one in ten planet-hosting stars ejects a Jupiter-like planet via a partial disruption, and adopting an average star formation rate of per year, at least one ejected planet in the recombination phase should be visible in the galaxy at any one time.

5. Conclusions

5.1. Limitations and Future Directions

The principle assumptions that we have made in this paper is that Jupiter-like planets are represented accurately by a polytropic model of its structure. One advantage of this model is that disruption simulations are trivially scalable to planets of a different size by a simple correction to , assuming that the planet’s mass interior to a given radius scales self-similarly and that the fluid remains unchanged. An polytrope reproduces the mass profile of coreless 1 planet relatively well, with the difference in never exceeding 10% throughout the planet’s interior (N. Miller, private communication).

The inclusion of a core of a few tens of affects out to a few times the core radius, for which . Beyond this radius, the structure of the planet is nearly identical to the coreless/polytropic models. This means that our simulations should be an accurate representation for disruptions where of the planet’s mass is removed for Jupiter-like planets. For Neptune-like planets, where the core mass can be larger than the gas mass, the difference in is substantial all the way to the planet’s outermost layers, and thus our simulation results should not be directly applied. As the average densities of Neptune-like planets is larger than Jupiter-like planets, for Neptunes should assume a smaller value.

Additionally, we assume that throughout the simulation volume, even for regions of very low density where the fluid is completely ionized and should behave as an ideal gas () or even a radiation pressure dominated fluid () in the lowest-density regions. This transition to different values of should affect the structure of the hot envelope that forms from the re-accreted debris that surrounds a partially disrupted planet, which is dynamically unimportant but may affect the planet’s observable signature. This is not to say that a more realistic equation of state would not affect the mass loss itself. As the process of ripping material from the planet involves rapid fluid decompression, a decrease in may result in slightly altered disruption dynamics.

Ideally, one would like to extend the models we have presented here to include a more physical equation of state that can treat all components of the pre- and post-disrupted planet realistically. As the resolution required to determine for multiple-orbit encounters beyond what we have presented here is prohibitive, it seems that the exploring the affects of using a more-complete equation of state with a realistic initial planet model is the next natural step for future studies. In the case of planets with a substantial core, these modifications are necessary to determine with any confidence.

5.2. The Fates of Scattered Jupiters

The fate of a Jupiter-like planet after a strong scattering event is a function of the strength of the tidal forces it experiences at periastron. In this paper, we have determined the disruption radius for Jupiter-like planets which sets the boundary between long-term survival and rapid tidal disruption. Below, we summarize the various post-scattering outcomes in order of decreasing distance, using and the tidal radius as points of reference.

Stalled : The planet is deposited into an orbit where the rate of tidal dissipation is too small to result in a change in semi-major axis over the lifetime of the system. This planet may be in a Kozai state driven by a third body in the system, or could experience another strong scattering event, which may lead to an increase of eccentricity and subsequent circularization.

Circularization/Migration : In this region, the planet is close enough to its parent star that tidal dissipation is effective, and the planet can circularize in yr or less for moderate values of . For near-radial orbits, circularization may still be longer than the stellar age, but again the Kozai mechanism or scattering could lead to a more rapid orbital evolution. All currently observed hot Jupiters are either stalled, in the process of circularizing, or already have circular orbits. If the planet is close enough to its parent star and , the planet will raise a tide on the star and migrate inwards due to the transfer of angular momentum.

Ejection : A planet that passes within the exclusion zone will be ejected from the system if its initial orbit is radial enough such that its orbital energy is significantly smaller than the self-binding energy of the planet. Slightly less than half of the planet’s initial mass remains bound to the central star, carrying with it a large reservoir of angular momentum that can significantly alter the host star’s spin rate and axis of rotation. The ejected planet will remain as bright or brighter than its host star for a few years, eventually plateauing via hydrogen recombination as an object with a fraction of solar luminosity for a century. All Jupiter-like planets that scatter in from beyond such that will be ejected from the system.

Complete disruption : For planets that are deep within their parent star’s potential well, the planet cannot soak the change in energy required to significantly alter the orbit, which eventually leads to its complete disruption. Approximately half of the planet’s mass accretes onto the stellar host, carrying the same specific angular momentum as in the ejection case, leading to even more pronounced effects on the stellar spin. Only planets that have migrated close to their stars prior to being scattered are destroyed before they are ejected.

Collision with the central star : The planet strikes the surface of the star directly. Anywhere from half to all of the planet’s mass is absorbed by the star, with the angular momentum being carried by this material potentially being smaller than that carried by the debris from a disruption, depending on how direct the impact is. These events are approximately twice as uncommon as ejections/disruptions.

We have benefited from many useful discussions with E. Ford, K. Schlaufman, E. Quataert, F. Rasio, G. Laughlin, N. Miller, D. Fabrycky, B. Hansen, M. Rees, and A. Socrates. The software used in this work was in part developed by the DOE-supported ASCI/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. Computations were performed on the Pleiades and Laozi UCSC computer clusters. We acknowledge support from the David and Lucille Packard Foundation (JG and ER-R), NSF grants PHY-0503584 (ER-R) and AST-0908807 (DL), NASA grants NNX07A-L13G (DL), NNX07AI88G (DL), NNX08AL41G (DL), and NNX08AM84G (DL); and the NASA Earth and Space Science Fellowship (JG).

Appendix A Modified PPM Gravity Algorithm

Figure 13.— Cartoon showing the sequence of states used to compute the evolution over a single a time-step using our modified gravity solver. The initial state is showing in the left panel, the middle panel shows the intermediate state , and the last panel shows the final state . is the position of the point mass representing the star, is the planetary core’s true center of mass, is the center of mass of the complete system. The variable vector field represents the value of all the conserved quantities in a given state. Differences between the states are exaggerated for illustrative purposes.

Because the binding energy of a planet on a disruptive orbit is comparable to the planet self-binding energy, the conservative properties of a code used to investigate planetary disruption are important. As the simulation of a partially-disruptive encounter involves the simultaneous resolution of both a compact core and two debris tails which are hundreds of time larger than the core, we found that the standard methods used to calculate the gravitational potential in a tidal disruption are too computationally expensive given the required accuracies.

Our approach was to improve upon the gravity solver found within the FLASH hydrodynamics code (Fryxell et al. 2000) such that it is better suited to investigating the problem of tidal disruption, a pictorial representation of the algorithm described below is shown in Figure 13. Two centers of mass are calculated at each time-step, which are defined as


where is the “true” center of mass and is the “core” center of mass, which only includes matter above a density cut-off , where is set to . The planet’s virtual particle is fixed to spatially coincide with the vector at all times. The planet’s self-gravity is calculated using a multipole expansion of the planet’s mass about instead of , which allows us to better approximate the planet’s potential using less terms in the multipole expansion.

Our algorithm is particularly well-suited for investigating tidal disruptions. If the expansion were performed about after large tidal tails have formed within the simulation, would be mostly determined by the position of material that is far away from the core, and the region with the best resolution of the potential may lie in empty space. We can estimate the number of terms required to represent the core’s potential if the expansion is carried out about instead of . Assume the core has a radius , and that the core lies a distance from the true center of mass. The angular scale of a lobe corresponding to a spherical harmonic of degree is simply , which means that even a first-order approximation of the core’s potential requires an expansion with . Assuming of a planet’s mass is lost during an encounter and this material lies an average distance from the core at apocenter, is on order , meaning that the multipole expansion must be carried out to . This is highly impractical, and thus it is much more efficient to carry out the multipole expansion about the planet’s core, whose position is associated with the densest material in the simulation and is where the potential gradients are largest.

The potential calculated from the multipole expansion about is used both to apply forces to the fluid in the simulation domain and to the virtual star and planet particles. A consequence of not expanding about the true center of mass is that there exists a non-zero force that is applied to the core. These forces are associated with the odd- multipole terms that usually cancel when the expansion is carried out about the true center of mass. Our multipole expansion does not discard these odd- terms, which allows us to confidently represent the fluid’s potential using an expansion of relatively low order.

In the FLASH code’s split PPM formalism, the equations used for coupling hydrodynamics and the gravitational field for a cell along each of the three cartesian directions are (Bryan et al. 1995)


The acceleration is calculated by extrapolating and to obtain an estimate for


This is in turn used to calculate


where is defined as


to enforce monotonicity in .

Because the potential of the star in our simulations is approximated by an analytical expression (for our simulations, a monopole potential), we can implement the following modification such that the component of the acceleration attributed to the star can be calculated to much higher accuracy. We first make the assumption that the matter distribution remains fixed over the course of the time-step, and then integrate the virtual particle positions forward in time from to . This allows us to calculate


an estimate of the star’s contribution to the potential based on the star’s approximate final position. This estimate should be much closer to the true value than simple extrapolation as the potential at a given location has a non-trivial time-dependence. Splitting into two components (star) and (planet), Equation (A5) becomes


Additionally, a correction must be made to Equation (A6)


where is the acceleration experienced by the core due to the presence of odd- terms in the multipole expansion. Using instead of , a hydrodynamical step is performed according to Equations (A3) and (A4), which yields and thus the true contribution at of the planet to the potential, . The positions of the virtual particles are then re-integrated from to , but this time using a linear interpolation of the time-evolving potential over the time-step.

Figure 14.— Accumulation of numerical error in the total energy and angular momentum for the simulations ran to compare with the results of FRW. The solid lines show simulations for which , the dashed lines show simulations where , and the dash-dotted lines show simulations where . Note that the rate of error accumulation is greatest shortly after periastron when the planet’s self-potential is rapidly varying as a function of time. As the planet returns to a state of quasi-equilibrium, the rate of error accumulation asymptotes to the original rate.

As mentioned previously, a complication introduced by this method is that the net force applied to the point about which the multipole expansion is carried out is non-zero. This is because we do not cancel out the acceleration applied to the extended object’s true center of mass , but rather the center of mass defined by the densest material, . While our correction to the gravitational acceleration (Equation (A10)) mitigates the problem somewhat, the mass density within the simulation domain varies as a function of time, and thus the total mass that satisfies the density criteria to be included when calculating also changes as a function of time. This leads to small, but measurable, changes in the reference point over the course of a simulation which are not entirely corrected by simply subtracting off the force applied at . To ensure that angular momentum and energy are conserved, a correction is made to virtual particle positions over the course each time-step when re-integrating from to ,


In practice, and are very nearly equal, with a substantial difference only arising when the extended object is close to being destroyed. Typical corrections are significantly smaller than the size of individual grid cells in the most highly-resolved regions.

The accumulation of error for our single passage encounters is shown in Figure 14. In the deepest encounters where the self-potential changes rapidly, we find that the rate of relative error accumulation in total energy and angular momentum is no larger than per dynamical time for our simulations with the lowest maximum-resolution, with the planet’s initial diameter being resolved by 50 grid cells at . Most of this error arises from the ejected tidal tails, which are resolved at lower resolution out of necessity because of the large volume they occupy. For encounters in which all the matter within the simulation is always resolved at highest resolution (i.e. those that do not have significant mass loss), the relative error accumulation is only per dynamical time.

We also tested the error accumulation for two sets of simulations with higher maximum resolutions, with the planet’s initial diameter being resolved by 100 and 200 grid cells. These simulations show an accumulation of fractional error of and per dynamical time, respectively.


  1. affiliation: Theoretical Astrophysics Santa Cruz (TASC), Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA;
  2. affiliation: NASA Earth and Space Science Fellow
  3. affiliation: Theoretical Astrophysics Santa Cruz (TASC), Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA;
  4. affiliation: Theoretical Astrophysics Santa Cruz (TASC), Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA;
  5. Assuming , . The minimum and maximum values are calculated using the lower and upper limits for as compiled by Schlaufman (2010).
  6. Assuming , . The minimum and maximum values are calculated using the lower and upper limits for as compiled by Schlaufman (2010).


  1. Antonini, F., Lombardi, J., & Merritt, D. 2010, eprint arXiv:1008.5369
  2. Arras, P. & Socrates, A. 2010, ApJ, 714, 1
  3. Barnes, S. A. 2003, ApJ, 586, 464
  4. Barnes, S. A. & Kim, Y.-C. 2010, ApJ, 721, 675
  5. Bicknell, G. V. & Gingold, R. A. 1983, ApJ, 273, 749
  6. Bodenheimer, P., Lin, D. N. C., & Mardling, R. A. 2001, ApJ, 548, 466
  7. Bryan, G. L., Norman, M. L., Stone, J. M., Cen, R., & Ostriker, J. P. 1995, Computer Physics Communications, 89, 149
  8. Carter, B. & Luminet, J. 1983, A&A, 121, 97
  9. Carter, B. & Luminet, J. P. 1985, MNRAS, 212, 23
  10. Cowling, T. G. 1941, MNRAS, 101, 367
  11. Diener, P., Frolov, V. P., Khokhlov, A. M., Novikov, I. D., & Pethick, C. J. 1997, ApJ, 479, 164
  12. Dobbs-Dixon, I., Lin, D. N. C., & Mardling, R. A. 2004, ApJ, 610, 464
  13. Eggleton, P. P., Kiseleva, L. G., & Hut, P. 1998, ApJ, 499, 853
  14. Evans, C. R. & Kochanek, C. S. 1989, ApJ, 346, L13
  15. Faber, J. A., Rasio, F. A., & Willems, B. 2005, Icarus, 175, 248
  16. Fabrycky, D. & Tremaine, S. 2007, ApJ, 669, 1298
  17. Ford, E. B. & Rasio, F. A. 2008, ApJ, 686, 621
  18. Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661
  19. Foucart, F. & Lai, D. 2010, arXiv
  20. Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics
  21. Frolov, V. P., Khokhlov, A. M., Novikov, I. D., & Pethick, C. J. 1994, ApJ, 432, 680
  22. Fryxell, B., Olson, K., Ricker, P., Timmes, F. X., Zingale, M., Lamb, D. Q., MacNeice, P., Rosner, R., Truran, J. W., & Tufo, H. 2000, ApJ, 131, 273
  23. Garaud, P. 2010, arXiv
  24. Gough, D. O. & McIntyre, M. E. 1998, Nature, 394, 755
  25. Guillochon, J., Ramirez-Ruiz, E., Rosswog, S., & Kasen, D. 2009, ApJ, 705, 844
  26. Guillot, T., Stevenson, D. J., Hubbard, W. B., & Saumon, D. 2004, In: Jupiter. The planet, 35
  27. Hebb, L., Collier-Cameron, A., Triaud, A., Lister, T., Smalley, B., Maxted, P., Hellier, C., Anderson, D., Pollacco, D., Gillon, M., Queloz, D., West, R., Bentley, S., Enoch, B., Haswell, C., Horne, K., Mayor, M., Pepe, F., Segransan, D., Skillen, I., Udry, S., & Wheatley, P. 2010, ApJ, 708, 224
  28. Hubbard, W. B. 1984, New York
  29. Hut, P. 1980, A&A, 92, 167
  30. Ida, S. & Lin, D. N. C. 2004, ApJ, 604, 388
  31. —. 2008, ApJ, 685, 584
  32. Irwin, J. & Bouvier, J. 2009, IAU Symposium, 258, 363
  33. Ivanov, P. B. & Novikov, I. D. 2001, ApJ, 549, 467
  34. Ivanov, P. B. & Papaloizou, J. C. B. 2010, MNRAS, 407, 1609
  35. Jackson, B., Barnes, R., & Greenberg, R. 2009, ApJ, 698, 1357
  36. Jurić, M. & Tremaine, S. 2008, ApJ, 686, 603
  37. Khokhlov, A., Novikov, I. D., & Pethick, C. J. 1993a, Astrophysical Journal v.418, 418, 181
  38. —. 1993b, Astrophysical Journal v.418, 418, 163
  39. Kobayashi, S., Laguna, P., Phinney, E. S., & Mészáros, P. 2004, ApJ, 615, 855
  40. Kochanek, C. S. 1992, ApJ, 385, 604
  41. —. 1994, ApJ, 422, 508
  42. Kozai, Y. 1962, AJ, 67, 591
  43. Kumar, P. & Goodman, J. 1996, ApJ, 466, 946
  44. Lai, D., Rasio, F. A., & Shapiro, S. L. 1994, ApJ, 437, 742
  45. Laine, R., Lin, D., & Dong, S. 2009, ApJ, 685, 521
  46. Lee, W. H., Ramirez-Ruiz, E., & van de Ven, G. 2010, ApJ, 720, 953
  47. Li, S.-L., Miller, N., Lin, D. N. C., & Fortney, J. J. 2010, Nature, 463, 1054
  48. Lin, D. N. C., Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606
  49. Lodato, G., King, A. R., & Pringle, J. E. 2009, MNRAS, 392, 332
  50. Lorén-Aguilar, P., Isern, J., & García-Berro, E. 2009, A&A, 500, 1193
  51. Malmberg, D., Davies, M. B., & Heggie, D. C. 2010, arXiv, astro-ph.EP
  52. Mardling, R. A. 1995, ApJ, 450, 722
  53. Matsumura, S., Takeda, G., & Rasio, F. A. 2008, ApJ, 686, L29
  54. Matsumura, S., Thommes, E. W., Chatterjee, S., & Rasio, F. A. 2010, ApJ, 714, 194
  55. Miller, N., Fortney, J. J., & Jackson, B. 2009, ApJ, 702, 1413
  56. Nagasawa, M., Ida, S., & Bessho, T. 2008, ApJ, 678, 498
  57. Nolthenius, R. A. & Katz, J. I. 1982, Astrophysical Journal, 263, 377, a&AA ID. AAA032.066.185
  58. Ogilvie, G. I. & Lin, D. N. C. 2004, ApJ, 610, 477
  59. Press, W. H., Flannery, B. P., & Teukolsky, S. A. 1986, Cambridge: University Press
  60. Press, W. H. & Teukolsky, S. A. 1977, ApJ, 213, 183
  61. Ramirez-Ruiz, E. & Rosswog, S. 2009, ApJ, 697, L77
  62. Rathore, Y. 2005, Proquest Dissertations And Theses 2005. Section 0037, 23
  63. Rees, M. J. 1988, Nature, 333, 523
  64. Robertson, B. E., Kravtsov, A. V., Gnedin, N. Y., Abel, T., & Rudd, D. H. 2009, MNRAS, 1774
  65. Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S. 2010, arXiv
  66. Rosswog, S., Ramirez-Ruiz, E., & Hix, W. R. 2008, ApJ, 679, 1385
  67. —. 2009, ApJ, 695, 404
  68. Schlaufman, K. C. 2010, ApJ, 719, 602
  69. Shen, Y. & Turner, E. L. 2008, ApJ, 685, 553
  70. Springel, V. 2010, MNRAS, 401, 791
  71. Takeda, G. & Rasio, F. A. 2005, ApJ, 627, 1001
  72. Tasker, E. J., Brunino, R., Mitchell, N. L., Michielsen, D., Hopton, S., Pearce, F. R., Bryan, G. L., & Theuns, T. 2008, MNRAS, 390, 1267
  73. Triaud, A. H. M. J., Cameron, A. C., Queloz, D., Anderson, D. R., Gillon, M., Hebb, L., Hellier, C., Loeillet, B., Maxted, P. F., Mayor, M., Pepe, F., Pollacco, D., Ségransan, D., Smalley, B., Udry, S., West, R. G., & Wheatley, P. J. 2010, arXiv, astro-ph.EP
  74. Usami, M. & Fujimoto, M. 1997, ApJ, 487, 489
  75. Watson, C., Littlefair, S., Diamond, C., Cameron, A. C., Fitzsimmons, A., Simpson, E., Moulds, V., & Pollacco, D. 2010, arXiv
  76. Zakamska, N. L., Pan, M., & Ford, E. B. 2010, MNRAS, 1566
  77. Zhou, J.-L., Lin, D. N. C., & Sun, Y.-S. 2007, ApJ, 666, 423
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description