Giant planet formation in stellar clusters: the effects of stellar flybys
M. M. Fragner  R. P. Nelson
Astronomy Unit, Queen Mary University of London, Mile End Road, London E1 4NS, UK
Received 6 April 2009 / Accepted 23 July 2009
Abstract
Context. The majority of stars in the Galaxy are thought to have formed within stellar clusters, resulting in occasional close encounters between stars during the epoch of planet formation. Encounters between young stars which possess protoplanetary discs will cause significant modification of the disc structure, and perturb any planets forming within the disc.
Aims. The primary aim of this work is to examine the effect of parabolic stellar encounters on the evolution of a Jovianmass giant planet forming within a protoplanetary disc. We consider the effect on both the mass accretion and the migration history as a function of encounter distance.
Methods. We use a gridbased hydrodynamics code to perform 2D simulations of a system consisting of a giant planet embedded within a gaseous disc orbiting around a star, which is perturbed by a passing star on a prograde, parabolic orbit. The disc model extends out to 50 AU, and parabolic encounters are considered with impact parameters ranging from 100250 AU.
Results. In agreement with previous work, we find that the disc is significantly tidally truncated for encounters <150 AU, and the removal of angular momentum from the disc by the passing star causes a substantial inflow of gas through the disc. The gap formed by the embedded planet becomes flooded with gas, causing the gas accretion rate onto the planet to increase abruptly. Gas flow through the gap, and into the inner disc, causes the positive inner disc torques exerted on the planet to increase, resulting in a sustained period of outward migration. For weaker interactions, corresponding to an encounter distance of 250 AU, we find that the planetdisc system experiences minimal perturbation.
Conclusions. Our results indicate that stellar flybys in young clusters may significantly modify the masses and orbital parameters of giant planets forming within protostellar discs. Planets that undergo such encounters are expected to be more massive, and to orbit with larger semimajor axes, than planets in systems which have not experienced parabolic encounters.
Key words: accretion, accretion disks  planets and satellites: formation
1 Introduction
Extrasolar planets are now being discovered with a broad range of masses and orbital element distributions (Udry & Santos 2007; Fischer et al. 2008). Accompanying these discoveries, there has been a significant resurgence in theoretical studies aimed at understanding how planetary systems form and evolve. The widely accepted picture of planet formation that proceeds via the formation and growth of planetesimals, which collide to produce fully formed planets, is a model which is essentially local, and makes no significant reference to the environment within which the host protostar and its protoplanetary disc are embedded. Most stars in the Galaxy, however, are thought to form in stellar clusters consisting of hundreds to thousands of stars, and close encounters between these stars are expected to significantly modify the structure of the disc (Korycansky & Papaloizou 1995; Larwood 1997; Clarke & Pringle 1991). It is likely that the impact of this, and other environmental factors, will need to be included in planet formation models eventually if the observed distributions of planetary characteristics are to be reproduced.
The question of how often young planetary systems
and their host stars suffer a
close encounter with passing stars in a cluster has
been addressed by Adams et al. (2005). They performed Nbody simulations
to determine the probability for a flyby to occur with a
given distance of closest approach (impact parameter).
Their models included the effects of cluster gas (present for the first 5 Myr)
and its dispersal after this time, and the encounter distances
were computed over a total evolution time of 10 Myr.
As expected, the presence of gas increases the cluster
density, and results in an increased encounter rate.
Adams et al. (2005) considered both virial (hot) and subvirial (cold)
clusters. For typical clusters sizes of 0.53.65 pc,
and stellar members between
N=1001000, it was found that
approximately one encounter with an impact parameter of 100 AU
occurred every 10^{6} years for subvirial clusters.
The cold configuration represents the most interesting case,
as it corresponds to young clusters in which forming planetary
systems are likely to be in a stage of early evolution,
such that the gas discs have not yet been dispersed, and
gas accretion onto the planet is still in progress.
The stellar density in the Adams et al. (2005) model was taken from
catalogs by Lada & Lada (2003) and Carpenter (2000), and
their range of cluster sizes corresponds to about 60of the observational sample.
If we adopt a value of N=500 for a typical cluster size,
and a typical disc life time of
years, then
we can crudely estimate that 1% of stars will experience an
encounter with impact parameter 100 AU.
Only half of these will be prograde, and very
few will be coplanar with the protoplanetary disc. Our 2D
simulations clearly cannot capture the warping and twisting of the disc
that would result from an inclined encounter, but
an inclination angle of 45
still results in a disc
experiencing 70% of the maximum possible gravitational perturbation
in its horizontal plane. The tidal truncation and spiral wave excitation
that we capture in 2D simulations will therefore remain as important phenomena
in the full 3D case even when the inclination angle is
much larger than the disc opening angle.
In a similar study of stellar encounters,
Malmberg et al. (2007) performed Nbody simulations for clusters
without gas and with virial initial velocities.
Their model did not include the gas component, so the
life time of their clusters is substantially increased,
and their simulation evolution times are
about a factor 10 longer than those of Adams et al. (2005).
For an initial cluster membership number of N=700 and
size 0.38 pc, they find an encounter rate of about
10 encounters in 10^{6} years after 5 Myr of evolution,
for impact parameters 100 AU.
This larger encounter rate is due to the
higher star density used in the simulations.
Once again, this encounter rate must be halved when considering
prograde encounters only. Taking this into account, and the
reduced influence of very high inclination encounters,
estimates of the fraction of stars experiencing
encounters which may induce significant disc truncation and
spiral wave excitation are in the range 0.251.5% for a disc
life time of 5 Myr.
An analytical estimate of the stellar
encounter rate was presented by Binney & Tremaine (1987):
where n is the star density in the cluster, is the velocity dispersion, is the close encounter distance and is the total mass of bodies involved in the encounter. The term in brackets is the gravitational focusing factor. For large encounter distance one can neglect the focusing factor and . We find that this formula agrees roughly with the encounter rate obtained by both Adams et al. (2005) and Malmberg et al. (2007).
In this paper we address the question of how such encounters will modify the evolution of a giant planet forming in a protoplanetary disc. It is well known that a giant planet forming in a disc surrounding a single star will form a tidallytruncated gap, and will migrate inward via type II migration on a time scale of a few 10^{5} years (e.g. Lin & Papaloizou 1986; Nelson et al. 2000). The planet will also accrete gas which slowly diffuses through the gap at a rate which is approximately one Jovian mass per 10^{5} yr. A closeencounter between the protoplanetary disc and a passing star will significantly perturb the disc, causing it to be tidally truncated and initiating an inward flow of gas. We use a gridbased hydrodynamics code (NIRVANA) to perform 2D simulations that examine how the evolution of a giant planet forming in a disc is modified by such encounters. We are particularly interested in calculating how the mass accretion and migration rates are modified.
Our simulations start with a Jovian mass planet on a circular orbit at 5 AU, which has formed a gap and is slowly undergoing inward type II migration in a disc of radius 50 AU. The effects of stellar perturbers with distances of closest approach to the planethosting star of between 100250 AU are examined. For close encounters in particular we find that accretion of gas can be substantially enhanced, and the inward migration can be reversed such that the giant planet undergoes a sustained period of outward migration.
The paper is organised as follows. In Sect. 2 we present our model, and in Sect. 3 the numerical method is described. The results of our simulations are presented in Sect. 4, and in Sect. 5 we present our conclusions.
2 Basic equations
The system we consider consists of a central star,
a thin protoplanetary accretion disc,
whose heighttoradius ratio ,
a
giant planet whose orbit plane coincides with the disc midplane, and
a stellarmass perturber on a prograde parabolic orbit whose plane
is coincident with that of the planet. It is therefore convenient to
consider a two dimensional system in which the equations
describing the disc dynamics are vertically integrated.
We adopt cylindrical polar coordinates (r, ,
z)
whose origin is located at the position of the central star,
where the rotation axis of the disc coincides with the z axis,
and the velocities are
denoted
,
,
0).
The vertically integrated continuity equation is given by:
(1) 
The vertically integrated radial momentum equation is:
and the angular momentum equation reads:
Here denotes the surface density, P is the height integrated pressure, is the gravitational potential, and f_{r} and are the viscous forces in the r and direction, respectively:
T_{ij} is the (i, j) component of the height integrated viscous stress tensor. We adopt a locally isothermal equation of state:
with , where is the constant aspect ratio and is the local Keplerian velocity. The gravitational potential experienced by the disc is given by:
(4) 
where is the mass of the central star, and the last two terms result from our choice of a non inertial reference frame centered on the primary star. In this problem we have two bodies in addition to the central star: the planet (i=1) and the secondary stellar perturber (i=2). The gravitational potential due to body i, denoted by , is given by:
(5) 
where M_{i} is the mass of the secondary star or planet. We use a softening parameter, b_{i}, whose value is b_{i}=0.6 H_{i}, and H_{i} is the effective disc thickness at the location of the companion object, i.
The planet and the secondary star
experience forces from all the other bodies and the
disc. The equation of motion for each of these objects is:
=  
(6) 
The last two terms are indirect terms, which arise because we work in a non inertial frame centered on the primary star. The second term describes the gravitational interaction between the planet and secondary. The third term represents the acceleration from the disc. Note that in order to conserve angular momentum, the torque experienced by the disc due to the planet or secondary star is used to calculate the force on that body, which is why the third term appears in a slightly unfamiliar form. We note further, however, that the forces due to the disc arising from material within the Hill sphere of the planet were neglected, resulting in a slight asymmetry between the disc and planet accelerations.
3 Numerical Methods
The system of equations described in Sect. 2 is integrated using NIRVANA, a gridbased hydrodynamic code (Ziegler & Yorke 1997). The advection scheme uses a second order accurate monotonic transport algorithm (van Leer 1977), and conserves mass and angular momentum globally. In all simulations the number of grid cells employed was ( . The equations of motion for the planet and stars are integrated using a standard firstorder Euler scheme.
3.1 Boundary conditions
A zero gradient outflow boundary condition was applied at
the outer radial boundary, and at the inner radial boundary
we applied a boundary condition which allows outflow, but at
a rate that is limited to be smaller
or equal to 10 times the viscous flow
rate, where we define the viscous mass flow rate
to be
Limiting the outflow rate in this way was found to be necessary as standard open or reflective boundary conditions were found to generate unsatisfactory results for this particular application.
3.2 Units
The equations are integrated in dimensionless form with the gravitational constant G=1. The inner radius is located at r=1 and the outer radius of the computational domain is at r=150. We assume that the inner disc edge is located at a distance of 0.5 AU from the star when converting between code and physical units. The mass of the central object is assumed to be 1 . When discussing simulation results we use the planet initial orbital period as our unit of time.
3.3 Initial conditions
We initialise the simulations using a power law distribution
for the surface density profile, modulated by a
taper function that reduces the surface
density near the outer edge of the disc:
.
Effectively, the power law profile extends between
the inner disc edge and
,
beyond which
the taper function becomes important and reduces the
surface density down to a small value. Once the density
falls to
,
the disc joins onto
a region where the density uniformly has this small value.
This setup provides the disc with plenty of space to respond to
the impact by the parabolic perturber, thus reducing
the influence of the outer boundary condition.
The taper function takes the form
(7) 
where . We chose the effective disc edge scale height to be , as was done by Korycansky & Papaloizou (1995). Correspondingly, the initial azimuthal velocity profile is chosen to be a modified Keplerian, which ensures the disc is in radial equilibrium initially:
(8) 
The term involving the aspect ratio represents the effect of the pressure gradient generated by a disc. The additional term h represents the effect of the taper function and can be written as:
h  =  
=  (9) 
Here it becomes clear why the small but nonzero quantity was introduced. It was chosen to be about 5 times the grid spacing.
For the uniform density region beyond , the azimuthal velocity was taken to be Keplerian.
The proportionality factor for the surface density prescription, , is determined from the total disc mass, which was taken to be at the beginning of the simulation. The radial velocity was chosen to be zero initially, and the disc aspect ratio is . The viscosity coefficient was chosen to be a constant throughout the disc, corresponding to at the planet position.
Running of the simulations consisted of two phases. First, a Jupiter mass planet was evolved for 535 orbits, starting at r_{1}=10, in order to create a clean gap. Accretion onto the planet was switched off during this phase. Second, a solar mass stellar perturber was launched (usually at a distance r_{2}=1000) and was evolved on a prograde, parabolic orbit in the plane of the disc. We denote the distance of closest approach of the perturber (or impact parameter) as q, which we express in units of the initial disc outer radius. Thus a run with q=2, for example, has a distance of closest approach equal to 200 measured in code units, as the disc has an initial radius equal to 100.
In the case of model 1, which we discuss below and which did not include a parabolic perturber, we simply continue the run on from t=535 orbits, but switch on gas accretion.
3.4 Planet accretion routine
Due to low resolution in the planet Hill sphere, we adopt
a simple treatment of the gas accretion process (Kley 1999).
For disc cells which lie within a distance of
,
where
is the planet Hill sphere radius, the amount of mass
removed from the disc at each time step and added to the
planet is given by:
where is the time step size. This prescription leads to an exponential decrease of the mass contained in the inner region of the Hill sphere, and we get a halfemptying time of:
The accretion parameter, , was allowed to take values [0,0.05,0.5], corresponding to no, slow and fast accretion.
3.5 Calculating specific torques on the planet
The specific torques on the planet (body 1) are
due to the direct acceleration from the disc, the direct
acceleration due to the secondary perturber and the
indirect gravitational force.
Mathematically this can be expressed as:
(11) 
The specific disc torque can be expressed as:
=  
=  (12) 
such that the specific disc torques are separated into contributions from the outer and inner disc, and from the corotation region. In the following section we introduce the total corotation torque, , and here we note that . We write the specific torques due to the direct gravitational interaction with the secondary perturber (body 2) and the indirect force as:
(13) 
The last two terms are due to the indirect force. Note that the first two terms in this expression almost cancel out, as . Therefore this term is dominated by the indirect force due to the acceleration of the central star by the disc.
3.6 Calculating torques in corotation region
An issue that we explore in this paper is the role of corotation torques in driving the migration of the planet. The parabolic flyby is expected to truncate the disc and drive a significant mass inflow. It is possible that as mass flows through the coorbital region of the planet, it induces a positive corotation torque, inducing a period of outward type III migration as described in Masset & Papaloizou (2003).
In order to measure the torque exerted on material located in the
corotation region in the numerical simulations,
we define an annulus whose boundaries extend by one Hill radius
on either side of the planet semimajor axis.
We thus split the disc into three parts, an inner and outer disc,
and the annulus  which we refer to as the corotation region.
Integrating Eq. (3) over an annulus with outer radius
and inner radius
gives:
(14) 
with:
(15) 
where J denotes the total angular momentum content in the annulus and ( ) is the angular momentum flux through the outer (inner) boundary of the annulus. Note that the integrand is evaluated at the radial boundaries of the corotation region only. T^{dF} denotes the angular momentum change due to a possible differential angular momentum flux. represents the viscous interaction between the annulus and the inner and outer disc, which occurs along the surfaces of the boundaries, and is the torque arising from the indirect forces. In practice we find that these latter two terms are negligible. The quantity is the gravitational torque exerted by any body ion the corotation region:
(16) 
In practice we find that the torque due to the parabolic perturber (i=2) also has negligible magnitude. Also there is an associated loss of angular momentum due to gas accretion onto the planet, which we will treat as an effective torque and denote as :
(17) 
where the integral is to be performed over the surface area enclosed within a distance of from the planet, as described in Sect. 3.4. Therefore the entire angular momentum budget can be expressed as:
(18) 
In order to estimate the error associated with this estimate, we also calculate:
(19) 
where Q is a measure of the error. In our simulations we find that Q remains small, but varies between 10^{3}10^{2}. The main contributor to the error arises from the exclusion of the Hill sphere region when calculating the torques on the planet, since we did not include this when calculating .
4 Results
We now present the simulation results. We begin by discussing the effect of the perturber on the global disc structure, and how the stellar flyby modifies the global angular momentum content of the disc. We then present a calibration run of a giant planet embedded in an otherwise unperturbed disc in model 1 to demonstrate agreement with previous results on migration and accretion rates. Following this we present results which demonstrate the effects that a stellar flyby can have on the evolution of a giant planet forming in a protostellar disc. The runs we have performed, and their associated model parameters, are summarised in Table 1.Table 1: Table of runs.
In model 2 the planet was maintained on a fixed circular orbit. In the second column of Table 1 the close encounter distance in units of the outer edge of the physical domain (50 AU) is listed. The third column of Table 1 contains the variation of the accretion parameter, where a value of 0.5 corresponds to fast, 0.05 to slow and 0to no accretion.
4.1 Effect of parabolic perturbers on the global disc structure: models 28
We now consider simulations for which there is an encounter between a stellar perturber on a prograde, parabolic orbit and a circumstellar disc. We focus here on the effects that the flyby has on the disc structure. It is expected that close encounters for which the impact parameter will produce significant changes in the global structure of the disc. Flyby scenarios with gaseous discs have been investigated without an embedded planet by Ostriker (1994) using linear theory, Korycansky & Papaloizou (1995) using linear theory and two dimensional simulations performed using a finite difference code, and by Larwood (1997) using three dimensional SPH simulations. These authors all agree that the disc looses a significant fraction of its angular momentum due to the encounter. The orbital motion of a stellar perturber on a parabolic orbit scans a broad (formally infinite) range of angular frequencies, and Korycansky & Papaloizou (1995) demonstrate that the peak response of the disc occurs for frequencies corresponding to inner Lindblad resonances which are located near the edge of the disc. More specifically, most of the torque experienced by the disc is applied at the point of pericentre passage. We begin by examining the effect of a parabolic encounter on the global structure of the disc in model 2. To initiate this run we took the results of model 1 after t=535 planet orbits, and launched a solarmass perturber from an initial position of (x, y)=(750, 0), chosen so that its pericentre distance from the central star is r_{2}=200 code units (or 100 AU in physical units), such that q=2. At the time when the perturber is launched, the planet has opened a deep gap in the disc, as may be observed in the first panel of Fig. 1. The influence of the perturber is apparent in the second panel of Fig. 1, which corresponds to a time of t=576 orbits when the perturber has reached pericentre (located at x=93.4, y=178.4). A lopsided twoarmed spiral wave is launched from near the outer edge of the disc, generating a perturbation with strong contributions from azimuthal mode numbers m=2 and m=1. The m=2 inner Lindblad resonance associated with the perturber's angular frequency at pericentre is located near the disc edge, and inward propagating m=2 spiral waves are launched in the disc. The strong m=1 contribution comes from the fact that the closest approach occurs on just one side of the disc, such that the perturbing potential has a significant non resonant m=1 component.
For a parabolic (zero energy) orbit, the angular frequency of the
perturber at closest approach in a frame centered on the primary is:
(20) 
where is the initial disc outer radius, and the second equality arises because we assume G=1, M_{*}=1, and M_{2}=1. The potential due to the secondary can be written:
where the first term in brackets describes the effect of working in the primary frame. For a circular orbit, r_{2} and are timeindependent in a frame corotating with the perturber, but for a parabolic orbit, r_{2} and become time dependent. Thus if one wants to consider contributions to the potential in the Fourier domain we must write in terms of frequencydependent quantities by means of a Fourier transform in time:
(21) 
Thus there is an infinite spectrum of frequencies which contribute to the potential of the perturber. We have calculated the Fourier transform of the potential , for reasons discussed below, using an analytical expression for the perturber trajectory from (Larwood 1997):
(22) 
with
(23) 
With this it is possible to determine the potential experienced in a ring at any radial location, r, in the disc for any chosen time, allowing us to evaluate the discrete Fourier transform of , and hence the coefficients .
Figure 1: This figure shows surface density contours for different stages of the encounter during model 2. The top left panel shows the initial state at the time when the perturber is introduced, just after the gap clearing phase after 535 planetary orbits. The top right panel corresponds to when the perturber has reached pericentre. In the bottom left panel the disc is being truncated, and waves travel all the way to the planetinduced gap. The lower right panel shows the disc a long time after the encounter has occurred. 

Open with DEXTER 
The fact that we have an infinite range of frequencies means that in principle m=2 waves are triggered at inner Lindblad resonances everywhere in the disc by the perturber. It is expected, however, that the disc responds dominantly to a frequency whose inner Lindblad resonance is located near to the edge of the disc where the gravitational perturbation is large. For each simulation we performed with different values of q (q=2, 2.5, 3, 5), it was found that the dominant frequency of the spiral waves launched were associated with Lindblad resonances located at different disc locations. An interesting question is what determines the location from which the dominant wave is launched? One factor in determining this is the fact that the disc is of finite radial extent, and a surfacedensity taper is applied at the outside edge, such that strong forcing by the companion's gravity there does not necessary excite a wave which propagates with large amplitude into the main body of the disc. By plotting the product , where is the initial unperturbed surface density, we find for most of our models that the radial location in the disc where this product reaches a maximum agrees very well with the position of the Lindblad resonance associated with the dominant m=2 wave. The value of used when plotting corresponds to the frequency of a wave which would be launched at an mfold Lindblad resonance located at r. For the model 2, with q=2, we find that the dominant wave excited near the edge of the disc had a frequency , with corresponding Lindblad resonance at . This coincides with the maximum of .
The disc response is highly non linear in model 2. Fourier analysis of the disc surface density shows that in addition to the dominant wave excited near the disc edge, a higher frequency m=2 wave is more pronounced in the disc inner parts, whose Lindblad resonance is located at . The implication is that the highly non linear wave excited at the very disc edge is damped through the process of truncating the disc to a radius of 50. The wave which survives and continues to propagate inward is launched near to the radial location where the disc is eventually truncated. The truncation of the disc down to a radius of (half its original size) can be observed in the third and fourth panels of Fig. 1, and arises because the non linear wave excited at the disc edge has a lower angular frequency than the local disc material.
The third panel of Fig. 1 also shows that spiral waves excited in the disc travel all the way to the gap generated by the planet, and the associated angular momentum loss by the disc causes significant mass flow through the gap and into the disc that lies interior to the planet. This flooding of the gap is temporary, however, as the planet eventually clears the gap on a time scale of a few hundred orbits after the perturber has passed by the disc. The fourth panel corresponds to a time when the perturber is no longer influencing the disc, and the spiral waves that it excited have propagated toward the disc centre and have dissipated. It is apparent that the embedded planet has been able to reform the gap, but the outer gap edge has developed an eccentric shape (with eccentricity 0.2.) which undergoes slow retrograde precession. The origin and evolution of this disc eccentricity are discussed later in the paper.
During pericentre passage of the secondary, mass gets pulled out through the outer boundary, and is lost from the disc. Figure 2 displays the relative mass and angular momentum loss by the disc in units of the initial quantities. As one can see, the disc loses about 6% of its initial mass, and there is an associated advective angular momentum loss through the outer boundary which is included in the solid curve showing total angular momentum loss. We also plot the angular momentum change due to torques only, as this result can be compared with Korycansky & Papaloizou (1995). Their setup differed from ours, as we use an isothermal equation of state and open boundaries, and their calculations adopted closed boundary conditions and a polytropic equation of state. In spite of these differences, we obtain similar results to those authors, and find that the angular momentum loss due to torques by a perturber passing the disc with a pericentre distance equal to twice the disc radius is approximately 11 percent. Korycansky & Papaloizou (1995) obtained the value 13.3 percent.
Figure 2: Angular momentum and mass of the disc normalised by the initial values for model 2. 

Open with DEXTER 
The results for model 6, with impact parameter q=2.5, are similar to those described for model 2. The disc response is highly nonlinear, and the disc is significantly truncated down to a radius of from an initial value of r=100. A dominant m=2 wave is excited at the disc edge whose frequency corresponds to a Lindblad resonance located at , and this again coincides with the maximum of . Further into the disc, interior to the final truncation radius, more prominent m=2 waves with higher frequencies are observed, corresponding to Lindblad resonances located near r=75.
A significant difference in behaviour is observed for the discs with , as the disc response transitions from being highly nonlinear to being quasilinear, with very modest tidal truncation. These are models 7 and 8, for which the impact parameters were q=3 and q=5, respectively. In the case of model 7, we find that a dominant m=2 wave is excited in the disc, but with a frequency which corresponds to a Lindblad resonance located at , somewhat inside the outer disc edge. This again corresponds to the maximum of the product . Unlike in the nonlinear cases discussed above, this wave propagates all the way in to the inner disc without being fully dissipated.
Figure 3: Disc angular momentum versus time for the following runs: model 5  black (solid) line; model 6  red (dotted) line; model 7  green (dashed) line: model 8  blue (dashdotted); model 1  magenta (tripledotdashed) line. 

Open with DEXTER 
The results of model 8 with q=5 are intriguing. We observe that a weak m=2 spiral wave is excited in the disc, and propagates all the way to the centre of the disc where the planet is. The dissipation of this wave causes a modest modification of the disc surface density in the vicinity of the planet and gap. This wave has a frequency whose inner Lindblad resonance location is predicted to be at . The product does not accurately predict the location of this Lindblad resonance, which is outside of the main body of the disc where the density is very small, and as such we would not expect to see a wave in the disc which has been excited there. It is possible that this wave is excited because of finite resonance width effects.
The change in angular momentum for each of the disc models 5, 6, 7 and 8 are shown in Fig. 3. As expected we see that the closest encounters lead to the largest loss of angular momentum from the disc, ranging from approximately 12% for q=2, down to approximately 1% for q=3. The case with q=5 is almost indistinguishable from model 1, for which there was no stellar perturber, only an accreting planet in the disc. The results are in good agreement with those of Korycansky & Papaloizou (1995) who found angular momentum changes of 13.3% for q=2, approximately 1% for q=3, and effectively 0% for q=5.
4.2 Planet calibration run: model 1
We ran a model in which the stellar perturber was absent in order to calibrate the migration and gas accretion rate of a giant planet that is initially of . The initial stages of this simulation consisted of running the model with a non accreting planet for 535 orbits, and after this time we switched on accretion so that gas entering the planet Hill sphere and being accreted did so after viscously diffusing through the tidallytruncated gap.
At the end of the simulation (t=1880 planet orbits) we
obtained an accretion rate of
in code units.
This corresponds to
/orbit,
and is similar to the accretion rates found by
Kley (1999), Nelson et al. (2000), and Lubow et al. (1999) (whose results
were in the range
to
/orbit).
As noted by Lubow et al. (1999), accretion onto the planet is highly efficient,
such that we find that
This compares well to results found by Kley (1999) who obtained , and Lubow et al. (1999) found using a slightly different disc and accretion model. The evolution of the planet mass in this case is compared against runs for which there is a stellar flyby later in this paper, and may be observed in Fig. 17. In this model the planet mass increases by about 50 per cent over the duration of the simulation between the times 535 to 1880 orbits.
As expected for planets in the type II migration regime, the migration time scale (defined by ) is found to be years (see Fig. 8), in good agreement with Nelson et al. (2000), D'Angelo & Kley (2003) and a simple estimate of the viscous time scale. Note that due to the initial gap clearing phase the planet has already migrated to r=9.43 after 535 orbits. The eccentricity stays below 0.01 throughout the run, as shown in Fig. 14.
Figure 4: Specific torques on planet for model 2, which are denoted by the following line styles:  black (solid) line;  blue (dotted) line;  red (dashed) line;  green (tripledotdashed) line. 

Open with DEXTER 
4.3 Planet on fixed circular orbit: model 2
Before we examine the results of simulations which consider the migration and gas accretion by a giant planet embedded in a disc which is perturbed by a stellar flyby, we examine a simulation in which the planet is maintained on a fixed circular orbit. As we will see in later sections, the evolution of both disc and planet in these systems is complicated, especially in those runs for which the flyby induces a strong perturbation in the disc, so a planet on a fixed circular orbit provides a simpler starting point for which we can analyse the torques acting on the planet, and the evolution of the disc.
The simulation we present here is model 2 (see Table 1), which involves a non accreting Jupiter mass planet ( M_{1}=10^{3}), and is a continuation of model 1 in which a non accreting giant planet opens up a gap in an otherwise unperturbed disc. We restart model 1 at time t=535 orbits, by which time the giant planet has migrated from r=10 to r=9.43and has opened a clean gap. Migration of the planet is switched off, and a stellar companion is introduced with a closest approach distance to the central star equal to twice the disc radius (i.e. q=2). The global effect of the flyby on the disc has been described in Sect. 4.1.
Figure 5: Mass flow through boundaries of corotation region, where a negative value corresponds to inflow. The black (solid) line corresponds to the outer boundary, and the red (dotted) line to the inner boundary. 

Open with DEXTER 
The specific torques experienced by the planet due to the disc are shown as a function of time in Fig. 4. This figure shows the specific torques due to the inner and outer disc, those originating from the corotation region that we defined in Sect. 3.6, and the total torque. Note that the torques have been smoothed over a temporal window of about six planetary orbits.
Figure 6: The top left panel shows how the gap outer edge penetrates the corotation region, represented by the white solid lines. The top right panel illustrates the mass weighting procedure used to estimate the orbital elements of the disc material located at the gap outer edge. The dashed lines are the boundaries of the region, from where the velocity and density information was taken, and the solid line shows the ellipse thus obtained which represents the gap outer edge. The lower left panel shows the azimuthally integrated surface density in the corotation region. The lower right panel displays the gap structure around the planet. Due to eccentricity near the gap outer edge, the integrated surface density in the outer disc is reduced, but is increased in the corotation region, which is denoted by the dashed vertical lines. 

Open with DEXTER 
The simulation results are shown from the point in time when the stellar perturber was introduced (at t=535 orbits), after the planet has opened a gap. Each of the torques acting on the planet, shown in Fig. 4, is seen to be approximately constant until orbits, which is just after the point of closest approach of the stellar flyby, whose effect is to truncate the disc severely and induce an inward mass flow through the disc. Shortly after t=580 orbits we see that the corotation torque, , spikes upward from an initial value . This can be understood in the context of type III migration theory (Masset & Papaloizou 2003), which predicts that a strong inward gas flow through the orbital location of a planet can induce a positive torque due to the gas interacting gravitationally with the planet in the corotation region as it passes from one side of its orbit to the other. As we will see below, the gas flow in our simulations does not remain on circular orbits, and becomes rather complicated, such that a simple application of type III migration theory is probably not valid for this problem. But we should still expect that an externally induced mass flow through the corotation region will lead to a positive corotation torque acting on the planet. Figure 5 displays the mass flow measured at the outer and inner edges of the corotation region, with the convention that a negative mass flow corresponds to flow directed inwards. Between the times 600670 orbits it is clear that there is a significant mass flow from the outer into the inner disc through the corotation region (though it can also be seen that some of the material flowing in through the outer boundary of the corotation region does not make it all the way through this region and out through the inner boundary). The correlation between the negative mass flow in Fig. 5 and the positive corotation torque in Fig. 4 suggests that the positive corotation torque we measure is induced by gas flow through the planet orbit. It is also clear, however, that the positive corotation torque is very shortlived, and the long term evolution of the planet torques are dominated by other effects.
Over the same time scales over which we see the development of a positive corotation torque, we also see that the inner disc torques exerted on the planet increase (upper line in Fig. 4) because of mass flow into the inner disc, thereby increasing its mass and gravitational influence. We can also see that there is a long term trend for the outer disc torques to become less negative, and this is due to the modification of the outer disc structure caused by the flyby and interaction with the planet. The surface density in the outer disc has decreased due to the flyby and associated mass flow into the inner disc, as shown in the lowerright panel of Fig. 6.
In addition to this longterm weakening of outer disc torques, and strengthening of inner disc torques, we see that the torque measured as a corotation torque remains highly variable. To a large degree, this is due to the fact that the outer edge of the gap has become eccentric after the flyby, and during periods of time when the disc pericentre penetrates the corotation region, the outer disc exerts a strong negative torque on the planet. Of course, this torque is not a corotation torque induced by mass flow through the orbit of the planet, but is in fact a resonant torque between the outer eccentric disc and planet. We postpone detailed discussion about the origin and long term evolution of this eccentric disc until Sect. 4.4, where we discuss the evolution of a migrating planet embedded in a disc perturbed by a stellar flyby.
The eccentric outer disc is shown in closeup by the upperright panel of Fig. 6, and on the time scale of the runs presented here the magnitude of the disc eccentricity varies between values of 0.10.2, causing the outer gap edge to vary its distance from the planet. To illustrate this effect, we calculated the massweighted orbital elements of gas that is located in an annulus extending from the planet semimajor axis (r=9.43) out to a radius r=12.3, in order to get a representative ellipse which describes the shape of the outer gap edge. This is depicted in the upperright panel of Fig. 6. The pericentre distance of this representative ellipse is plotted as a function of time in Fig. 7, and by comparing with Fig. 4 we can see that when the pericentre distance decreases below the designated corotation torque, and the total torque, becomes negative. In particular this occurs at times approximately equal to 680, 780 and 1020 orbits.
To summarise the rather complicated interaction between the planet and the perturbed disc after the flyby, we find the following:
 (i)
 tidal truncation of the disc causes significant mass flow through the planet orbit from outer to inner disc;
 (ii)
 this mass flow can induce a shortlived positive corotation torque on the planet;
 (iii)
 the flow of mass from outer to inner disc causes a buildup of mass in the inner disc, and thus increases significantly the positive Lindblad torques exerted by the disc on the planet;
 (iv)
 the negative outer disc Lindblad torques are weakened by the change in disc structure, and reduction in surface density of the outer disc, induced by the flyby and interaction with the planet. This is such that over long times the positive inner disc torques are stronger than the negative outer disc torques;
 (v)
 on the time scale of the run presented here, a longlived and time dependent eccentric outer disc is generated and maintained. Close encounters between the outer eccentric gap edge and the planet induce strong, negative temporally varying torques on the planet.
Figure 7: The pericentre distance of the gap outer edge versus time for model 2. 

Open with DEXTER 
4.4 Migrating planet: model 3
In this section we investigate the effect of a stellar perturber on a migrating, nonaccreting planet embedded in a disc. As with model 2, which was discussed in the previous section, the initial conditions for this simulation were taken from model 1 at time t=535 orbits after the Jovian mass planet had opened up a clean gap in the disc. At this point in time the orbital radius of the planet . The perturber was introduced on a parabolic orbit with closest approach to the disc equal to 200 code units, such that the ratio of pericentre distance to disc radius was equal to two (i.e. q=2).
Our discussion of the results for model 2 presented in Sect. 4.3 suggest that a flyby with q=2 will result in net positive torques being exerted on the planet, first because the perturberinduced gas flow through the planet orbit will induce a shortlived positive corotation torque, and second (and more importantly for the longterm evolution) because structural changes in the disc will change the balance between inner and outer disc torques. This should drive outward migration of the planet. Figure 8 shows that this expectation is fulfilled in model 3 whose semimajor axis evolution is depicted by the solid black line. The outward migration lasts for more than 2000 planetary orbits, up until the end of the simulation.
Figure 8: Semimajor axes of planets as a function of time for various models. The line styles are as follows: model 3  black (solid) line; model 1  green (dashdotted) line; model 4  upper black (dotted) line; model 5  red (dashed) line. 

Open with DEXTER 
Figure 9: Specific torques experienced by the planet for model 3. The line styles are as follows:  black (solid) line;  blue (dotted) line;  red (dashed) line;  green (tripledotdashed) line;  magenta (dashdotted) line. 

Open with DEXTER 
Figure 10: This figure shows the disc surface density at different times during model 3. The top left panel shows the initial state of the disc after gap clearing. The top right panel shows that the gap outer edge has increased its semimajor axis and eccentricity. The lower left panel shows that the gap circularises after its initial eccentricity growth. The lower right panel shows the corresponding azimuthally averaged profiles. The vertical lines represent the position of the planet at the different times indicated. 

Open with DEXTER 
Accompanying the positive spike in the corotation torque we observe a sharp increase in the positive torque exerted by the inner disc on the planet, which is caused by mass flow into the inner disc. This increase in torque is seen to decline over time scales of 1000 orbits as the inner disc accretes onto the central star and the planet migrates outward. But long term changes to the structure of the outer disc cause the associated negative torque to decrease in magnitude. These changes are induced by tidal truncation, induced mass inflow, and eccentricity evolution. The result is that the total torque experienced by the planet remains positive for the duration of the simulation, driven mainly by the increase in inner disc density and mass.
As discussed in Sect. 4.3, the
outer disc in model 2 is found to become
eccentric shortly after the flyby, and the same effect is observed in
model 3 for the migrating planet. Analysis of the disc structure shows that
the disc eccentricity is localised near the outer edge of the gap between disc
radii in the range
in both models 2 and 3.
The surface density at different times for this simulation is
shown in Fig. 10. The first three panels show images of the
surface density distribution, and illustrate the disc eccentricity
evolution. The bottom right panel shows the azimuthally averaged
surface density as a function of radius, and shows the effect of the
eccentric disc on the gap width around the planet.
The origin of the disc eccentricity appears to be as follows:
the flyby induces a significant and rapid inward gas flow
through the disc. A fraction of this inflowing gas
approaches the planet on modestly non circular orbits and undergoes a
close gravitational interaction with the planet, which exerts a positive
torque on the gas near to its pericentre passage. The equations for
the change in semimajor axis and eccentricity of a particle
due to gravitational interaction with a planet can be written
(e.g. Murray & Dermott 1999):
where is the ratio of planet to central star mass, f is the true anomaly of the perturbed particle, E is the eccentric anomaly, is the radial acceleration and is the tangential acceleration. These expressions demonstrate that a fluid particle experiencing a significant torque at pericentre through interaction with the planet should increase both its semimajor axis and eccentricity (since , and ). Using the procedure adopted in Sect. 4.3 to calculate the mass weighted orbital elements of the gas near the outer edge of the gap, we have plotted the evolution of the semimajor axis and eccentricity of the outer gap edge in Figs. 11 and 12, respectively. Between the times we see that both of these quantities increase. The expressions given by Eq. (24) clearly do not govern accurately the complicated behaviour of a gas disc interacting with a planet, but they do indicate the route by which fluid elements can be scattered onto eccentric orbits with larger semimajor axes by an embedded planet. The fact that accretion discs are known to support longlived normal modes with azimuthal mode number m=1 (Papaloizou 2002) suggests that a coherent eccentric mode may arise from the processes described above.
Figure 11: Semimajor axis of planet  black (solid) line; semimajor axis of gap outer edge  blue (dotted) line; pericentre distance of outer gap edge  red (dashed) line. 

Open with DEXTER 
Figure 12: Planet eccentricity  black (solid) line; eccentricity of gap outer edge  blue (dotted) line. 

Open with DEXTER 
Once a significant disc eccentricity has been set up, then
it is expected that the planet and eccentric disc will undergo
secular interaction (in addition to continuing interaction
at Lindblad resonances), causing their eccentricities to vary
and their orbits to precess. The time evolution of the
planet semimajor axis and eccentricity are plotted along with those of
the disc near the gap outer edge in Figs. 11 and 12,
and the longitudes of pericentre for both disc and planet are
plotted in Fig. 13. We can make rough estimates for the
expected direction and rates of precession, in addition to
the eccentricity variations of disc and planet, using secular
perturbation theory, treating the disc (both interior and exterior to
the planet orbit) as smeared out rings with appropriate values for
the semimajor axes and eccentricities. We note that such
a theory, which neglects important effects such as the disc
pressure, cannot provide accurate predictions for the
observed evolution. But it is precisely this failure which
can be used to deduce which of the neglected physical processes
are important in driving the evolution of the the discplusplanet
system. Working to secondorder in the eccentricities, secular theory
for the time variation of outer disc and planet longitudes of pericentre
and eccentricities gives (e.g. Murray & Dermott 1999):
In the above expressions, the a_{x}, refer to the semimajor axes, where the subscript x is one of: od for the outer disc; id for the inner disc; p for the planet. Similarly the e_{x} terms denote the eccentricities, the terms denote the longitudes of pericentre, and the m_{x} terms denote the masses. The and are Laplace coefficients, which we evaluate using power series up to fifth and sixth order in , respectively. We denote and .
Figure 13: Longitude of pericentre for the planet  black (solid) line; longitude of pericentre for the gap outer edge  red (dotted) line. 

Open with DEXTER 
At time t=1000, we can see from Fig. 13 that the outer disc and planet are close to apsidal alignment, with the planet longitude of pericentre just leading that of the outer disc such that . From Eq. (25) we predict that the outer disc should precess in a prograde sense, and the planet in a retrograde sense (at the observed rate), using values for the disc and planet parameters that have been obtained from the simulations. In principle we could have read off values from Figs. 1113, but these apply only to material very near the outer edge of the gap (in fact the values plotted in Figs. 1113 were calculated in order to define the orbital elements of the fluid located at the gap edge), and the planet interacts gravitationally with material further out in the disc. We therefore extend the outer disc region used to calculate the orbital elements of the gas so that it covers the interval , where the larger value corresponds to the orbital distance at which the disc eccentricity becomes small. This procedure leads to the following values, which we used in Eq. (25): ; ; ; ; ; ; ; ; ; ; ; . Equation (25) predicts prograde precession for the outer disc and retrograde precession for the planet, but Fig. 13 shows that both the outer disc and planet precess in a retrograde sense at approximately the same rate, with their apsidal lines close to being aligned. This suggests that the disc and planet are in a joint mode in which the retrograde precession of the outer disc is being driven by pressure perturbations associated with the eccentric mode, and the planet precession rate is determined by secular gravitational interaction with the outer and inner disc (Papaloizou et al. 2001; Papaloizou 2002).
Examining the disc at t=1500 orbits we see that the semimajor axis and eccentricity of the outer disc located near the gap edge are decreasing. In addition the precession rate of both outer disc and planet are close to zero as the direction of precession is about to change. This change in precession direction is apparently due to the strength of the eccentric disc mode decreasing (indicated by the decrease in disc eccentricity) such that the pressure driven retrograde precession is overcome by the prograde precession induced by the planet gravity. This reduction in outer disc eccentricity also allows the planet to precess in the prograde sense as the term proportional to in the expression for d in Eq. (25) becomes small. Using the values ; ; ; ; ; ; ; ; ; ; ; (obtained as described above) in Eq. (25), secular theory predicts that the outer disc eccentricity at t=1500 should be increasing very slowly and the planet eccentricity should be decreasing, since . What we observe, however, is that and both decrease, indicating that disc eccentricity is damping due to some other process not associated with the secular exchange of angular momentum. This is most likely due to the disc viscosity, whose effects may be enhanced by the large eccentricity gradient which builds up in the disc at earlier times and which can increase the shear rate locally in the disc.
If we now consider the disc and planet evolution at t=1800 orbits, we see that the more rapid prograde disc precession has caused to change from being positive to being negative. As expected from Eq. (25), which predicts that the planet eccentricity should now begin to grow due to interaction with the outer disc, we see that starts to increase. Using the values ; ; ; ; ; ; ; ; ; ; ; in Eq. (25), however, secular theory predicts that the outer disc eccentricity should decrease at a slow rate at this point, which is the opposite of what we observe in Fig. 12. Clearly the disc eccentricity at this point is growing because of another effect, unrelated to the secular interaction between disc and planet. Previous work which has considered the growth of disc eccentricity has examined the role of non linear mode coupling, where coupling between the m=1 component of the planet potential and an m=1disturbance in the disc leads to the excitation of an m=2 wave at the 3:1 outer Lindblad resonance, whose pattern speed is equal to half the orbital angular velocity of the planet. The removal of angular momentum from the disc by this wave causes its eccentricity to grow. This process was first examined by Papaloizou et al. (2001), and was shown to operate for planets whose masses were greater than around 10 Jupiter masses. More recent calculations by Kley & Dirksen (2006) and D'Angelo et al. (2006) suggest that eccentric disc growth may occur for lower mass planets in the Jovian mass range. In order to investigate this issue further we have performed a Fourier analysis of the disc surface density distribution near to the 3:1 resonance and have observed an m=2 wave with the appropriate pattern speed, suggesting that the above explanation for the disc eccentricity growth may be valid. We note, however, that this wave is likely to be a combination of waves generated by the non linear mode coupling described above, and a wave excited directly by the planet at its 3:1 eccentric Lindblad resonance, since the planet has a small, nonzero eccentricity. Growth of the disc eccentricity may be assisted by the fact that the corotation resonances, normally associated with eccentricity damping, may be reduced in efficacy by the large gap size and modified gap structure that arises because of strong interaction with the planet by inflowing gas shortly after the flyby, when the eccentric disc was set up initially.
Figure 14: Planet eccentricity for various models. The line styles are as follows: model 3  black (solid) line; model 4  black (dotted) line; model 5  red (dashed) line; model 1  green (dashdotted) line. 

Open with DEXTER 
Over long time scales we see from Figs. 1113 that the disc and planet undergo a significant period of evolution in which their eccentricities grow and decay periodically, but with a long term trend of decreasing amplitude of variation. During this time the planet continues to migrate outward toward the gap edge, and the inner disc accretes onto the central star. This process should eventually lead to a situation where the negative outer disc torques begin to dominate over the positive inner disc torques, such that inward migration should ensue. The computational resources required to run the simulation for this length of time, however, are very substantial and are above our current capability. What is clear, however, is that a planet that forms in a disc around a star within a stellar cluster which experiences a close encounter with another member of the cluster may undergo a significant period of outward migration.
4.5 Variation of the planet accretion rate: models 35
We now consider a set of runs (models 35) where we varied the accretion rate of gas onto the planet embedded in a disc with a stellar perturber whose pericentre distance was equal to twice the disc radius (q=2). Model 3 was discussed in Sect. 4.4, and models 4 and 5 are identical apart from the inclusion of gas accretion. We distinguish between slow (  model 4) and fast (  model 5) accretion in Eq. (10), corresponding to half emptying times of t_{1/2}=0.07 and t_{1/2}=0.007 orbits, respectively. (In fact both of these accretion times are small so we should really call them ``fast'' and ``very fast'' accretion runs).
Figure 15: Specific torques on planet for model 4. The various line styles are denoted as follows:  black (solid) line;  blue (dotted) line;  red (dashed) line;  green (tripledotdashed) line;  magenta (dashdotted) line. 

Open with DEXTER 
We begin by considering model 4, corresponding to slow accretion. The change in semimajor axis is shown by the upper black (dotted) line in Fig. 8, where we see that the migration is directed outward and occurs at a rate that is slightly faster than was obtained for the non accreting planet (model 3). The torques experienced by the planet in model 4 are shown in Fig. 15. During the early phases of evolution the torques are very similar to those already discussed for model 3 and presented in Fig. 9. We see that the corotation torque spikes upward shortly after the flyby, and we also see that the positive inner disc torques increase rapidly at the same point in time due to gas flow into the inner disc. There is also a slow decrease in the negative outer disc torques which arises because of the gas flow from outer to inner disc and the change in disc structure near the planet after the flyby. The main difference between model 3 and 4 arises after approximately 900 orbits, when we see that the negative outer disc torque quickly reduces in magnitude, and this occurs because the outer disc becomes more eccentric in model 4 than model 3 due to the planet mass growing and being more able to drive the disc eccentricity upward. The increased disc eccentricity reduces the pericentre distance of the gap outer edge, and allows the planet to accrete gas from the outer disc at a fairly rapid rate. The reduction in outer disc mass and torques causes a significant imbalance between inner and outer disc torques (even though the inner disc is also being accreted by the planet, but at a slower rate), and explains why the slow accreting planet is able to migrate more rapidly than the non accreting planet in Fig. 8.
Figure 16: Specific torques on planet for model 5. The various line styles are denoted as follows:  black (solid) line;  blue (dotted) line;  red (dashed) line;  green (tripledotdashed) line;  magenta (dashdotted) line. 

Open with DEXTER 
We now consider the orbital evolution of the fastaccreting planet in model 5, whose semimajor axis as a function of time is shown by the red (dashed) line in Fig. 8. We see that over the duration of the run, the planet in model 5 migrates more slowly than those in models 3 (non accreting) and 4 (slowly accreting). The specific torques experienced by the planet for model 5 are shown in Fig. 16, and there are a number of differences which explain why the more rapidly accreting planet migrates outward more slowly. We see that the positive spike in the corotation torque at orbits is reduced in this model because the planet accretes more of the gas that tries to flow from outer disc to inner disc after the flyby. Consequently there is no sudden increase in the positive inner disc torques observed for models 3 and 4, since little mass makes it through into the inner disc. As in model 4, we find that the evolution of the negative outer disc torques includes a period of rapid reduction in the torque amplitude, which is even more pronounced in the case of model 5 than in model 4. This arises because the outer disc becomes highly eccentric in the vicinity of the outer gap edge, such that the outer disc pericentre reaches very close to the planet orbit. The rapidly accreting planet is able to quickly accrete gas from this eccentric outer disc, thus reducing the negative torque is exerts. The long term result of these effects is for the inner disc torque to be larger than the outer disc torque, thus driving outward migration, but the imbalance between inner and outer disc torques is reduced compared to models 3 and 4, leading to slower migration.
The planet masses are plotted as a function of time in Fig. 17. Model 4 is shown by the black (solid) line, model 5 is shown by the red (dotted) line, and the mass of the planet embedded in a disc without a perturber (model 1) is shown by the green (shortdashed) line. It is clear that the effect of the perturber is to significantly increase the accretion rate of an embedded giant planet, which is not surprising given that a close encounter drives a significant mass flow through the disc, causing the gap to be flooded. Over longer times, as the planet masses grow and tidal truncation becomes more effective, gas accretion slows down as the gap is reformed. At this point in time the accretion rates of the planets in the perturbed (models 4 and 5) and unperturbed (model 1) discs are very similar. The planet masses at these times, however, differ significantly because of the earlier evolution. Between times of 500 to 1500 orbits the planet in model 1 has reached 1.35 Jupiter masses, from an initial mass of 1 Jupiter mass, whereas models 4 and 5 have reached 2 and 2.7 Jupiter masses, respectively.
Overall, our models suggest that there should be significant differences between planets whose nascent discs have been significantly perturbed by a flyby compared with those whose discs have not been externally perturbed. In particular, planets born in highly perturbed discs should on average have higher masses and larger semimajor axes.
Figure 17: Evolution of planet masses for various runs. The line styles are denoted as follows: model 1  green (dashed) line; model 4  black (solid) line; model 5  red (dotted) line. 

Open with DEXTER 
4.6 Variation of the impact parameter: models 58
In this section we consider the effect of varying the impact parameter of the encounter, q, on the evolution of the planet. In each of the models we have allowed the embedded planet to both migrate and accrete gas (the accretion rate is set to be ``fast'' in these simulations).
Figure 18: Semimajor axes of planets from various runs. The line styles are denoted as follows: model 5  black (solid) line; model 6  red (dotted) line; model 7  green (dashed) line; model 8  blue (dashdotted) line; model 1  magenta (tripledotdashed) line. 

Open with DEXTER 
The semimajor axes versus time are shown in Fig. 18 for each of these models, and the corresponding torques experienced by the planets are shown in Fig. 20. The change in semimajor axis for the planet in model 5, with q=2, is shown by the solid line in Fig. 18, and the evolution of this planet has already been discussed in detail in Sect. 4.5. Here, we see that it is this planet which migrates outward at the fastest rate when compared with the other accreting planets for which the external perturbation was weaker (q > 2). We see that there is a general trend in models 58, such that outward migration is weakened when the external perturbation is weaker. Model 6 had impact parameter q=2.5, and the planet in this run is seen to undergo slower outward migration than the planet in model 5 (q=2). Model 7 had impact parameter q=3, and the outward migration in this case is not only slower than for models 5 and 6, but also appears to stall and reverse at time t=1400 orbits. Model 8 experienced the weakest perturbation (q=5), and the resulting migration is always inward. But we also notice that the rate of migration in this case is slightly slower than that for model 1 (no external perturbation), because even in this case the perturber excites a wave in the disc which modifies the surface density profile near the planet when it dissipates.
Figure 19: Planet mass from various runs. The line styles are denoted as follows: model 5  black (solid) line; model 6  red (dotted) line; model 7  green (dashed) line; model 8  blue (dashdotted) line; model 1  magenta (tripledotdashed) line. 

Open with DEXTER 
The torques experienced by the planets for models 58 are shown in Fig. 20. We recall from Sect. 4.5 that a significant feature of the torques for model 5 was the rapid decrease in outer disc torques due to the planet accreting gas from the outer eccentric disc. This led to sustained outward migration, due to the imbalance between inner positive disc torques and outer negative disc torques being skewed by this accretion of the outer disc. Examining the upper right panel of Fig. 20, we see that the negative outer disc torques for model 6 are also reduced, but no where near as sharply as for model 5. This is because the lower amplitude external perturbation does not lead to such a large disc eccentricity, and the outer disc is not so readily accreted. A positive imbalance between inner and outer disc torques still arises, however, such that outward migration is maintained. We note that the weaker external perturbation also causes the initial positive spike in the corotation torque to be reduced by .
Figure 20: Specific torques acting on planets for runs with varying q. The line styles are denoted as follows:  black (solid) line;  blue (dotted) line;  red (dashed) line;  magenta (tripledotdashed) line;  green (dashdotted) line. 

Open with DEXTER 
For model 7 (q=3), we notice that the initial spike in the corotation torque is so small that it does not induce an overall positive torque in the disc shortly after the flyby. In addition, there is no abrupt change in the outer disc torques, since the disc eccentricity remains relatively small in this case as the tidallyinduced inward mass flow is reduced. The small outer disc eccentricity, however, continues to favour accretion from the outer disc between times 6001000 orbits, leading to a period of outward migration. But this stalls as the planet moves outward and away from the inner disc which continues to accrete onto the central star.
As we have seen in Fig. 18, the planet in model 8 (q=5) behaves very similarly to the one in model 1, for which there was no external perturber. The long term decline in inner and outer disc torques shown in the lower right panel of Fig. 20 is due to accretion of the disc by the planet.
The mass accretion history for these models is shown in Fig. 19. As expected, the most strongly perturbed model results in the most rapid mass growth, as the gap is subject to greater flooding in that case. As with the migration, there is a global trend of smaller planetary mass growth for larger impact parameters, and model 8 (with q=5) shows very similar results to model 1 (without a perturber). Clearly, outward migration and enhanced accretion are strong functions of the impact parameter.
5 Conclusion
In this paper we have examined the influence of a parabolic stellar encounter on the evolution of a giant planet forming within a protostellar disc. Previous work by Adams et al. (2005) and Malmberg et al. (2007) has examined the frequency with which such encounters will occur within a young stellar cluster, and their estimates suggest that between 0.251.5% of stars should experience a prograde encounter with impact parameter 100 AU and inclination relative to the disc midplane over a 5 Myr protostellar disc life time, with the precise encounter rate determined by the stellar density.
The results of our simulations suggest that encounters whose distance of closest approach are <3 times the disc radius (i.e. q < 3) have a severe influence on the disc structure, due to the non linear response of the disc. This is in agreement with the previous work of Korycansky & Papaloizou (1995) and Larwood (1997). The disc models we have considered in this work have outer radii of 50 AU, such that encounters with pericentre distances <150 AU cause significant changes to the disc structure. The main effect of the encounter is significant shrinkage of the outer disc radius, and excitation of an inward propagating m=2 spiral wave at an inner Lindblad resonance which corresponds closely to the orbital angular frequency that the perturbing star has at pericentre.
The distance of closest approach which results in significant modification of the disc structure also corresponds, in our work, to the point at which significant changes in the evolution of the embedded giant planet occur. This is because the strong tidal truncation of the disc results in significant loss of angular momentum by the disc material. This induces a substantial inflow of gas through the disc. When this mass flows through the planet orbit it induces a shortlived episode of outward typeIII migration, causes the tidallytruncated gap around the giant planet to be temporarily flooded with gas, and increases the mass of the disc which lies interior to the planet orbit. Thus we find that for impact parameters , the mass accretion rate onto the giant planet is significantly increased, and a period of prolonged outward migration can be induced. This suggests that giant planets, which have been subject to a stellar encounter when forming in a protoplanetary disc, should have higher masses and larger semimajor axes, on average, than planets which have not been subject to such an encounter. The rate of outward migration, and the increase in accretion rate, is found to scale roughly with the inverse of the closest encounter distance, such that impact parameters of q=2 and 2.5 were found to cause the largest changes in giant planet evolution, with long periods of sustained outward migration which lasted for the duration of our simulations (approximately 2000 giant planet orbits). A q=3 encounter resulted a shorter period of outward migration, and a q=5 encounter resulted in almost no change in the evolution of the planet.
We have not considered discs with radii larger than 50 AU in this work, and so cannot comment directly on the effect that stellar encounters will have on giant planets forming in discs with larger radii. What is clear, however, is that discs provide a conduit through which the influence of a passing star can be communicated through to a forming planet via the inward propagation of a non linear spiral wave and associated inward mass flow. Observations indicate that disc radii can be significantly larger than 50 AU. For example, some of the discs which have been imaged in in Orion have radii of a few hundred AU (McCaughrean & O'Dell 1996). Given that the encounter frequency scales quadratically with the encounter distance when gravitational focusing is ignored (see Sect. 1), it is reasonable to speculate that the effect of passing stars may be even more important than we have stated in this paper during giant planet formation, especially if typical disc sizes are in excess of 50 AU.
There are a number of improvements that need to be made to the models we have presented before any serious attempt could be made to incorporate the environmental effects we have discussed into a meaningful planetary population synthesis calculation. The first is that we have only considered prograde, coplanar encounters, whereas we would expect the relative angular momentum vectors associated with encounters between passing stars and discs to be isotropically distributed within a young stellar cluster. The absence of Lindblad resonances in the disc during retrograde encounters means that these are likely to have a minimal effect on a forming planet, provided that the passing star does not actually pass through the disc. During an early phase of this project we ran simulations of retrograde encounters, and the results of these indicate that such perturbations will have minimal impact on planet formation. Inclined prograde encounters will have a weakened ability to tidally truncate the disc and induce inward mass flow, so we would expect their influence to be somewhat smaller than that described in this paper. We note, however, that inclined encounters will excite bending waves in the disc, causing the disc orbital plane to deviate away from that of the planet, at least for a short time after the impact until discplanet interactions can bring the system back into alignment. This will obviously modify both the migration and accretion history of the planet, but in ways that we are unable to quantify at present.
An additional improvement to the model presented here would include a broader survey of the available parameter space. This would include simulations of discs with varying outer radii, and planets with different masses and semimajor axes. Also, multiple planets in the process of forming coevally may be induced to undergo close interaction through the changes in migration outlined in this paper, leading to a rich variety of possible outcomes through gravitational scattering. These and other issues will be the subject of future work and publications.
Acknowledgements
The simulations presented in this paper were performed using the QMUL HPC Facility purchased under the SRIF initiative. We acknowledge comments received from an anonymous referee which led to improvements to this paper.
References
 Adams, F. C., Proszkow, E. M., Fatuzzo, M., & Myers, P. C. 2005 (In the text)
 Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, MNRAS, 341, 213 [NASA ADS] [CrossRef]
 Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton: Princeton Univ. Press), 541 (In the text)
 Carpenter, J. M. 2000, AJ, 120, 3139 [NASA ADS] [CrossRef] (In the text)
 Clarke, C. J., & Pringle, J. E. 1991, MNRAS, 249, 584 [NASA ADS]
 D'Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540 [NASA ADS] [CrossRef] (In the text)
 D'Angelo, G., Lubow, S. H., & Bate, M. R. 2006, ApJ, 652, 1698 [NASA ADS] [CrossRef] (In the text)
 Fischer, D. A., Marcy, G. W., Butler, R. P., et al. 2008, ApJ, 675, 790 [NASA ADS] [CrossRef]
 Frank, J., King, A., & Raine, D. 1992, Accretion Power in Astrophysics (Cambridge University Press), 2nd Edition
 Goldreich, P., & Tremaine, S. 1979, ApJ 233, 857
 Kley, W. 1998, A&A, 338, L37 [NASA ADS]
 Kley, W. 1999, MNRAS, 303, 696 [NASA ADS] [CrossRef] (In the text)
 Kley, W., & Dirksen, G. 2005, A&A, 447, 369 [NASA ADS] [CrossRef] (In the text)
 Korycansky, D. G., & Papaloizou, J. C. B. 1995, MNRAS, 274, 85 [NASA ADS]
 Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57 [NASA ADS] [CrossRef] (In the text)
 Larwood, J. D. 1997, MNRAS, 290, 490 [NASA ADS]
 Lin, D. N. C., & Papaloizou, J. C. B. 1993 (Tucson: University of Arizona Press), 749
 Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 [NASA ADS] [CrossRef] (In the text)
 Malmberg, D., De Angeli, F., Davies, M. B., et al. 2007, MNRAS, 000, 1 (In the text)
 Masset, F. S. 2001, ApJ, 558, 453 [NASA ADS] [CrossRef]
 Masset, F. S. 2002, A&A, 387, 605 [NASA ADS] [CrossRef] [EDP Sciences]
 Masset, F. S., & Papaloizou, J. C. B. 2003, ApJ, 588, 494 [NASA ADS] [CrossRef] (In the text)
 McCaughrean, M. J., & O'dell, C. R. 1996, AJ, 111, 1977 [NASA ADS] [CrossRef] (In the text)
 Ostriker, E. C. 1994, ApJ, 424, 292 [NASA ADS] [CrossRef] (In the text)
 Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 [NASA ADS] [CrossRef] (In the text)
 Papaloizou, J. C. B. 2002, A&A, 388, 615 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263 [NASA ADS] [CrossRef] [EDP Sciences]
 Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef]
 Throop, H. B., & Bally, J. 2007, AAS
 Udry, S., & Santos, N. C. 2007, A. R. A. A., 45, 397 [NASA ADS]
 van Leer, B. 1977, J. Comput. Phys., 23, 276 [NASA ADS] [CrossRef] (In the text)
 Ward, W. R. 1991, in 22d Lunar and Planetary Science Conference (Houston: LPI), abstract, 1463
 Ward, W. R. 1997, Icarus, 126, 261 [NASA ADS] [CrossRef]
 Ziegler, U., & Yorke, H. W. 1997, Comput. Phys. Commun., 101, 54 [NASA ADS] [CrossRef] (In the text)
All Tables
Table 1: Table of runs.
All Figures
Figure 1: This figure shows surface density contours for different stages of the encounter during model 2. The top left panel shows the initial state at the time when the perturber is introduced, just after the gap clearing phase after 535 planetary orbits. The top right panel corresponds to when the perturber has reached pericentre. In the bottom left panel the disc is being truncated, and waves travel all the way to the planetinduced gap. The lower right panel shows the disc a long time after the encounter has occurred. 

Open with DEXTER  
In the text 
Figure 2: Angular momentum and mass of the disc normalised by the initial values for model 2. 

Open with DEXTER  
In the text 
Figure 3: Disc angular momentum versus time for the following runs: model 5  black (solid) line; model 6  red (dotted) line; model 7  green (dashed) line: model 8  blue (dashdotted); model 1  magenta (tripledotdashed) line. 

Open with DEXTER  
In the text 
Figure 4: Specific torques on planet for model 2, which are denoted by the following line styles:  black (solid) line;  blue (dotted) line;  red (dashed) line;  green (tripledotdashed) line. 

Open with DEXTER  
In the text 
Figure 5: Mass flow through boundaries of corotation region, where a negative value corresponds to inflow. The black (solid) line corresponds to the outer boundary, and the red (dotted) line to the inner boundary. 

Open with DEXTER  
In the text 
Figure 6: The top left panel shows how the gap outer edge penetrates the corotation region, represented by the white solid lines. The top right panel illustrates the mass weighting procedure used to estimate the orbital elements of the disc material located at the gap outer edge. The dashed lines are the boundaries of the region, from where the velocity and density information was taken, and the solid line shows the ellipse thus obtained which represents the gap outer edge. The lower left panel shows the azimuthally integrated surface density in the corotation region. The lower right panel displays the gap structure around the planet. Due to eccentricity near the gap outer edge, the integrated surface density in the outer disc is reduced, but is increased in the corotation region, which is denoted by the dashed vertical lines. 

Open with DEXTER  
In the text 
Figure 7: The pericentre distance of the gap outer edge versus time for model 2. 

Open with DEXTER  
In the text 
Figure 8: Semimajor axes of planets as a function of time for various models. The line styles are as follows: model 3  black (solid) line; model 1  green (dashdotted) line; model 4  upper black (dotted) line; model 5  red (dashed) line. 

Open with DEXTER  
In the text 
Figure 9: Specific torques experienced by the planet for model 3. The line styles are as follows:  black (solid) line;  blue (dotted) line;  red (dashed) line;  green (tripledotdashed) line;  magenta (dashdotted) line. 

Open with DEXTER  
In the text 
Figure 10: This figure shows the disc surface density at different times during model 3. The top left panel shows the initial state of the disc after gap clearing. The top right panel shows that the gap outer edge has increased its semimajor axis and eccentricity. The lower left panel shows that the gap circularises after its initial eccentricity growth. The lower right panel shows the corresponding azimuthally averaged profiles. The vertical lines represent the position of the planet at the different times indicated. 

Open with DEXTER  
In the text 
Figure 11: Semimajor axis of planet  black (solid) line; semimajor axis of gap outer edge  blue (dotted) line; pericentre distance of outer gap edge  red (dashed) line. 

Open with DEXTER  
In the text 
Figure 12: Planet eccentricity  black (solid) line; eccentricity of gap outer edge  blue (dotted) line. 

Open with DEXTER  
In the text 
Figure 13: Longitude of pericentre for the planet  black (solid) line; longitude of pericentre for the gap outer edge  red (dotted) line. 

Open with DEXTER  
In the text 
Figure 14: Planet eccentricity for various models. The line styles are as follows: model 3  black (solid) line; model 4  black (dotted) line; model 5  red (dashed) line; model 1  green (dashdotted) line. 

Open with DEXTER  
In the text 
Figure 15: Specific torques on planet for model 4. The various line styles are denoted as follows:  black (solid) line;  blue (dotted) line;  red (dashed) line;  green (tripledotdashed) line;  magenta (dashdotted) line. 

Open with DEXTER  
In the text 
Figure 16: Specific torques on planet for model 5. The various line styles are denoted as follows:  black (solid) line;  blue (dotted) line;  red (dashed) line;  green (tripledotdashed) line;  magenta (dashdotted) line. 

Open with DEXTER  
In the text 
Figure 17: Evolution of planet masses for various runs. The line styles are denoted as follows: model 1  green (dashed) line; model 4  black (solid) line; model 5  red (dotted) line. 

Open with DEXTER  
In the text 
Figure 18: Semimajor axes of planets from various runs. The line styles are denoted as follows: model 5  black (solid) line; model 6  red (dotted) line; model 7  green (dashed) line; model 8  blue (dashdotted) line; model 1  magenta (tripledotdashed) line. 

Open with DEXTER  
In the text 
Figure 19: Planet mass from various runs. The line styles are denoted as follows: model 5  black (solid) line; model 6  red (dotted) line; model 7  green (dashed) line; model 8  blue (dashdotted) line; model 1  magenta (tripledotdashed) line. 

Open with DEXTER  
In the text 
Figure 20: Specific torques acting on planets for runs with varying q. The line styles are denoted as follows:  black (solid) line;  blue (dotted) line;  red (dashed) line;  magenta (tripledotdashed) line;  green (dashdotted) line. 

Open with DEXTER  
In the text 
Copyright ESO 2009