From stardisc encounters to numerical solutions for a subset of the restricted threebody problem
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
email: abreslau@mpifr.de
Received: 10 March 2015
Accepted: 21 November 2016
Various astrophysical processes exist, where the flyby of a massive object affects matter that is initially supported against gravity by rotation. Examples are perturbations of galaxies, protoplanetary discs, or planetary systems. We approximate such events as a subset of the restricted threebody problem by considering only perturbations of noninteracting lowmass objects that are initially on circular Keplerian orbits. In this paper, we present a new parametrisation of the initial conditions of this problem. Under certain conditions, the initial positions of the lowmass objects can be specified as being largely independent of the initial position of the perturber. In addition, exploiting the known scalings of the problem reduces the parameter space of initial conditions for one specific perturbation to two dimensions. To this twodimensional initial condition space, we have related the final properties of the perturbed trajectories of the lowmass objects from our numerical simulations. In this way, maps showing the effect of the perturbation on the lowmass objects were created, which provide a new view on the perturbation process. Comparing the maps for different massratios reveals that the perturbations by low and highmass perturbers are dominated by different physical processes. The equalmass case is a complicated mixture of the other two cases. Since the final properties of trajectories with similar initial conditions are also usually similar, the results of the limited number of integrated trajectories can be generalised to the full presented parameter space by interpolation. Since our results are also unique within the accuracy strived for, they constitute general numerical solutions for this subset of the restricted threebody problem. As such, they can be used to predict the evolution of real physical problems by simple transformations, such as scaling, without further simulations. Possible applications are the perturbation of protoplanetary discs or planetary systems by the flyby of another star. Here, the maps enable us, for example, to quantify the portion of unbound material for any periastron distance without the need for further simulations.
Key words: methods: numerical / protoplanetary disks / gravitation / planets and satellites: dynamical evolution and stability / scattering
© ESO, 2017
1. Introduction
Various astrophysical processes exist where matter, initially supported against gravity by rotation, is affected by the gravitational force of a passingby massive object. Among them are, for example, galaxygalaxy encounters and perturbations of protoplanetary discs, debris discs, or planetary systems by passing stars (see e.g. Toomre & Toomre 1972; Farouki & Shapiro 1981; Clarke & Pringle 1993; Heller 1995; Hall et al. 1996; Hall 1997; Kobayashi & Ida 2001; Melita et al. 2002; Levison et al. 2004; Pfalzner et al. 2005b; Olczak et al. 2006; Scharwächter et al. 2007; Malmberg et al. 2011; Craig & Krumholz 2013; Punzo et al. 2014). Some of these processes happen on different scales and might involve additional forces, such as viscous friction or magnetic fields. Nevertheless, what they all have in common is that, approximately at the closest approach of the perturber, the perturber’s gravitational force dominates the event. Therefore, it is very likely that, at least, matter that is close to the perturber at pericentre passage is in all cases affected in a similar way.
For simplicity, we focus here on cases where relatively lowmass matter is orbiting a massive object (the host) and is perturbed by another massive object (the perturber). As examples, these can be a single planet, or a disc of gas or debris orbiting a star that is perturbed by another star flying by. In the case of a disc, we follow the approach of several previous studies, (e.g. Toomre & Toomre 1972; Hall et al. 1996; Kobayashi & Ida 2001; Pfalzner et al. 2005b) and neglect interactions of the disc material with itself, such as selfgravity or viscosity, to approximate all these cases as a subset of the restricted threebody problem.
The investigation of the threebody problem has a long history starting with Newton, Euler, Lagrange, and Jacobi (see e.g. Marchal 1990; Valtonen & Karttunen 2006; Musielak & Quarles 2014, and references therein). A well studied subset of the general threebody problem is the perturbation of stable twobody systems by a third body. This problem was investigated analytically (e.g. Hut 1983; Heggie & Hut 1993) and numerically (e.g. Valtonen & Heggie 1979; Hut & Bahcall 1983; Heggie & Hut 1993). The perturbation of lowmass matter owing to a flyby is related to this problem. It was analytically approximated, for example, by Ostriker (1994), Larwood (1997), and D’Onghia et al. (2010). However, to date, no satisfying general analytical description of such perturbations has been found. Therefore, previous studies investigating the influence of an encounter on planetary systems or the global properties of discs (e.g. energy, mass, size) usually performed statistical numerical simulations that cover certain subsets of the large space of encounter parameters (e.g. Hall et al. 1996; Hall 1997; Pfalzner et al. 2005b; Fregeau et al. 2006; Steinhausen et al. 2012; Hao et al. 2013; Breslau et al. 2014; Li & Adams 2015).
Also, based on numerical simulations, some approaches were made to generalise the effect of these perturbations. As such, usually the properties of the perturbed system were related to initial properties.
Clarke & Pringle (1993) and Hall et al. (1996) investigated the perturbation of protoplanetary discs in equalmass encounters and presented the probability for a disc particle to have a certain fate (remaining bound, captured by perturber, unbound) depending on its initial distance to the host, r_{init}, normalised to the perturber’s periastron distance, r_{p,peri}^{1}. The perturber’s periastron distance here and in the following means the distance between host and perturber at the moment when the perturber is in the pericentre of its orbit.
In a similar scenario, Hall et al. (1996) identified the particles that undergo the largest energy change during an encounter and located them within the disc (see their Figs. 4a, 5, and 6). They found that, shortly before periastron passage, these particles formed four distinct groups that interact with the host and perturber in different ways. However, since the disc is already deformed at the time of periastron passage, a description of the spatial distribution of these particles proved difficult.
Toomre & Toomre (1972) performed simulations of galaxy perturbations with noninteracting lowmass disc particles. Since they modelled the gravitational potential of the involved galaxies by point masses, their results are comparable to simulations of lowmass discs around stars. Similar to Clarke & Pringle (1993) and Hall et al. (1996), they showed the fates (remaining bound, captured by perturber, unbound, forming tails, and a bridge between host and perturber) of the disc particles, but, depending on “the positions they would have reached at time t = 0 in the absence of forcing”. By “absence of forcing”, they mean without the force from the perturber on the particles and in their setup, t = 0 is the moment of the perturbers pericentre passage (see their Fig. 15 and the describing text).
Using a similar approach, like Toomre & Toomre (1972), we developed a new parametrisation of the initial conditions of the problem. For a given perturbation^{2}, the resulting twodimensional initial condition space enables the creation of maps of the properties of the system, e.g. the final values. This new view on the problem provides a better understanding of the underlying perturbation process itself. Since the solution for one set of initial conditions is unique, once obtained, numerical results can be easily applied to different physical problems by simple transformations, such as scaling, and thus making additional simulations unnecessary. This reduces the computational effort of parameter studies in this context significantly.
This paper focuses on the presentation of the new method, which is described in Sect. 2, along with the numerical scheme that was used. In Sect. 3, some exemplary results of the new method are provided in the form of the abovementioned maps. The limitations of the method along with the application to real physical systems are discussed in Sect. 4.
2. Method
Our new method is basically a refinement of the approach used by Toomre & Toomre (1972) for their Fig. 15. Since they mentioned it only in passing, and without complete description, we first explain their approach and then describe how we developed it further.
To describe it simply, we consider a disc of massless particles initially orbiting counterclockwise on circular Keplerian orbits in the xyplane around a host with mass M_{h} in the centre of the coordinate system. These particles are perturbed by a perturber of mass M_{p}. The actual orbit of the perturber depends on the massratio of host and perturber, m = M_{p}/M_{h}, the eccentricity of the orbit, e_{p}, the pericentre distance, r_{p,peri}, and the plane of the orbit relative to the orbit of the particles. We further restrict this to cases where the orbit of the particles and the orbit of the perturber are coplanar and prograde, thus reducing the problem to two dimensions. The line from the host to the pericentre of the perturber’s orbit defines the xaxis of our reference system. The particles’ gravitational effect on host and perturber is assumed to be negligible.
As in some previous studies (e.g. Toomre & Toomre 1972; Clarke & Pringle 1993; Hall et al. 1996), we scale all lengths of the problem by setting r_{p,peri} = 1. The initial radius of particle i’s orbit, r_{i,init}, is then given in units of r_{p,peri}.
We assume the perturber to pass the host on a parabolic orbit. Therefore, their distance is at t = −∞ and t = ∞ in principle infinite. In the moment of pericentre passage the distance is smallest, namely r_{p,peri}, and the perturber’s influence on the host is strongest. The gravitational force of the perturber on the particle is strongest close to pericentre passage. Since it is impossible to cover the period from t = −∞ to t = ∞ in a simulation, this fact can be used to restrict the investigation of the perturbation of the particle’s orbit to a time period around pericentre passage^{3}. The perturbed particle trajectories obtained for one specific perturber orbit converge with increasing initial separation between perturber and host.
Owing to this convergence, numerical investigations of the perturbation process should start when the perturber’s influence on the particles is still negligible. At that time (t = t_{init}), the perturber’s position is given in polar coordinates by (r_{p,init}, θ_{p,init}). Owing to the above’s condition, the perturber is then far outside the initial particle orbits, r_{p,init} ≫ r_{i,init}. For a parabolic orbit the time of flight from this initial position to pericentre is given by^{4}(1)with the standard gravitational parameter for the twobody problem μ = G(M_{h} + M_{p}) and the eccentric anomaly (2)Since the particles orbit the host initially on circular Keplerian orbits with radii r_{i,init}, their initial angular velocities are given by (3)where G is the gravitational constant. The particles’ orbital periods are T_{i} = 2π/ω_{i,0}. Without the influence of the perturber, the positions of the particles are given for all times by the initial radius r_{i,init} and the angle (4)During the period which the perturber needs to reach pericentre (t_{peri}−t_{init}), a particle that is initially (at t = t_{init}) at the position (θ_{i,init},r_{i,init}), would move without the perturber’s influence to (θ_{i,init} + ω_{i,0} (t_{peri}−t_{init}), r_{i,init})^{5}. This is the position Toomre & Toomre (1972) called “the position [it] would have reached at time t = 0 in the absence of forcing”. Since the particle will never be at this position, we will call it virtual pericentre position (VPP) in the following. This position is obtained from the initial position by the transformation Figure 1 shows an example of this transformation for a disc of particles perturbed by a perturber with mass M_{p} = M_{h}. The perturber passes the host on a parabolic orbit coming from the bottom of the image, passing pericentre at (1, 0), and leaving at the top of the image (see black line in Fig. 1a). The particles are coloured by their final fates: blue particles remain bound to the host, green particles become captured by the perturber, and red particles become unbound. Figure 1a shows the particles when the perturber has nearly reached pericentre (perturber position indicated by a black dot at the bottom of the image). Owing to the influence of the perturber, the disc is already deformed at that moment.
Fig. 1 Virtual pericentre positions (VPPs) for a perturbed disc. The perturber has the same mass as the host and the pericentre of the parabolic, prograde orbit is at (1, 0). In a) and b) the disc is seen face on when the perturber has nearly reached pericentre and at t = t_{init}, respectively. In c) the VPPs of the disc particles are shown. For full description see text. 

Open with DEXTER 
Figure 1b shows the same for t = t_{init}. Because of the differential rotation within the disc, the pattern from Fig. 1a appears “wound up”. In contrast to Fig. 1a, the disc is not yet deformed by the influence of the perturber.
Figure 1c shows the particles at their VPPs. In this case, the disc is neither deformed nor is the pattern from Fig. 1a wound up like in Fig. 1b.
2.1. New parametrisation of the initial conditions
Instead of setting up simulations by sampling the initial positions of the particles (at t = t_{init}) randomly and plotting the particles later at the VPPs, as shown for example in Fig. 1c, we invert the transformation given by Eqs. (5) and (6) to parametrise the initial conditions. This is analogous to the definition of the impact parameter for the description of scattering processes, e.g. of a testparticle in the field of a single attracting body. The impact parameter describes at which distance the testparticle would pass the target without being influenced. Similarly, we define the initial particle positions by the positions where the particles would be when the perturber is in pericentre without the perturbation. For a given initial distance between host and perturber, the initial particle positions are obtained by sampling the VPP space and transforming the coordinates with
2.2. The numerical simulations
For our simulations, we sampled the VPP space spanned by x_{vpp} and y_{vpp} with massless tracer particles. Both dimensions were divided into bins and sampled with particles, each in the centre of the respective bin. For the figures presented in Sect. 3, the region x_{vpp}/r_{p,peri} ∈ [−5;5], y_{vpp}/r_{p,peri} ∈ [−5;5] was sampled with 500·500 particles. The mass of the host is M_{h} = 1, the perturber with mass M_{p} = m moves on a parabolic orbit (e_{p} = 1) with a pericentre of r_{p,peri} = 1. The trajectory of each particle was integrated individually with the adaptive LSODA integrator from the ODEPACK library (Hindmarsh 1983).
The initial position of the perturber was determined individually for each particle. To fulfil the claim that the initial influence of the perturber is negligible, the initial distance of the perturber was chosen so that the force from the perturber onto the particle, F_{p}, was small relative to the force from the host, F_{h}: (9)This maximum initial influence of the perturber determines the initial position of the perturber. For a particle with , the minimum initial distance between particle and perturber, d_{p}, follows from Eq. (9): (10)For practical reasons, the initial distance between host and perturber is set to r_{p,init} = d_{p} + r_{vpp}. From r_{p,init} and the eccentricity of the perturber orbit, e_{p}, follows θ_{p,init} according to (11)and thus the time until pericentre passage, t_{peri}−t_{init}, with Eqs. (1) and (2). After determining the time until pericentre passage of the perturber, the initial position of the particle was obtained according to Eqs. (7) and (8).
Tests have shown that for ϵ = 10^{4} the particles are sufficiently unperturbed at the onset of the simulations. The obtained accuracy is discussed in Sect. 4. With these initial conditions, the perturber is initially at much larger distances from the host than in most simulations of this kind that have been performed so far since in earlier studies the computational resources necessary for a larger number of these type of simulations were not available.
During the simulation, it was regularly tested whether the particle had settled into a stable orbit. Here, possible stable, final orbits are either bound orbits (e ≤ 1) around host or perturber, or unbound orbits (e> 1) with the centre of mass. A bound orbit with the centre of mass is only a transition state. An orbit was accepted as stable when one of the above conditions was fulfilled and the orbital elements (eccentricity, semimajor axis, argument of periapsis, etc.) that were calculated for the respective orbit did not change more than 1% between two tests. The duration between two tests was usually (t_{peri}−t_{init})/50. When the orbit was stable, the integration of the particle’s trajectory was stopped.
3. Results
Fig. 2 Maps of the final properties of the particles depending on their VPPs for the region r_{vpp}/r_{p,peri}< 5 for a massration of m = 0.1. a) and b) show whether the particles are finally bound to host or perturber, or become unbound. a) Also shows the orbit of the perturber through the Cartesian space and b) the perturber’s interaction orbit (see text for explanation). c) Shows the final eccentricities of the bound particles, and d) shows the final semimajor axes of the bound particles relative to the radii of their initial orbits. 

Open with DEXTER 
Fig. 3 As for Fig. 2, but for a massratio of m = 20.0. 

Open with DEXTER 
From our simulations, we determined the fates of the particle trajectories: whether the particles remain bound to the host, are captured by the perturber, or become completely unbound. Additionally, the final orbital elements of the particles were determined. They define the final orbits completely. Depending on the particles VPPs, these final properties can be presented as maps. Because we want to present the new parametrisation of the initial conditions, along with some consequences here, we only show some of the final properties for three sample massratios.
3.1. Lowmass perturbation
Figure 2 shows the maps for the region with r_{vpp}/r_{p,peri} < 5 for the massratio m = 0.1. This corresponds to the perturbation of a disc with initial radius of r_{disc} = 5 r_{p,peri}, or a perturbation with r_{p,peri} = 0.2 r_{disc}, respectively. Initially, the particles were orbiting the host (indicated by the black dot at (0, 0)) counter clockwise. The perturber’s pericentre is at (0, 1).
Figure 2a shows the particles’ fates: particles with VPPs in the blue region remain bound to the host, particles from the green region^{6} are captured by the perturber, and particles from the red region become unbound. The black line in Fig. 2a shows the orbit of the perturber – the perturber came from the bottom edge of the image (as indicated by the black arrow), moved to (1, 0), and left across the top edge. The white regions indicate particles for which the accuracy achieved with the used integration method is not sufficient, and which are therefore not shown (see below for further explanation and Sect. 4).
It is often assumed that there is a simple radial dependence on whether particles remain bound or not. Figure 2 shows that this is certainly not the case, but a much more complicated dependence exists. It can be seen that only particles with VPPs in a narrow band along the approaching branch of the perturber’s orbit (black line in Fig. 2a) become unbound. On the side of the departing branch of the perturber’s orbit, particles with VPPs from a broader region become unbound. The only region from where particles are captured by the perturber is close to the perturber’s pericentre (green area around (1, 0)). Owing to the low perturber massratio, the region from where particles become unbound is relatively small.
As can be seen in Fig. 2a, the shape of the perturber orbit cannot be used to explain, for example, the shape of the region from where particles become unbound. The reason is that the perturber orbit shows all positions of the perturber for a time span while the map shows the fates of particles which are at a certain position at a certain moment. Even under the assumption that the particles are not influenced by the perturber at all, particles that meet the perturber are not shown near the perturber orbit in the map. When the perturber is at the coordinates (r_{1}, θ_{1}) at the time t_{1}, a particle being exactly at this position at that time is shown in the map at the position (r_{1}, ). These positions are highlighted for all t_{1} with the black line in Fig. 2b and are referred to in the following as the perturber’s “interaction orbit”.
Figure 2b shows that the interaction orbit matches the shape of the region from where particles become unbound much better than the real perturber orbit. Comparing Figs. 2a and b, we can also see that part of the white dots in Fig. 2a are very close to the interaction orbit. This means that particles with VPPs in these regions are directly on the way of the perturber and approach the perturber very closely. This results in integration errors, which is why these particles are excluded from the results.
Besides the fates of the particles, the shapes of their final orbits are important to describe the outcome of the perturbation process. The shapes of the orbits are defined by their eccentricities and semimajor axes. These two properties can also be transferred into energy and angular momentum to obtain another, often used, view on the perturbation process.
For the particles that are finally bound to the host or perturber, the final eccentricities (e < 1) and semimajor axes relative to the initial radii of their orbits are shown in Figs. 2c and d. Since the unbound particles are of minor importance for most applications, their final orbital elements can be found in Fig. A.1.
The final eccentricities of the particles that remain bound to the host can be considered as nearly symmetrical to the xaxis (see Fig. 2c). Close to the region from where the particles become unbound (white region in Fig. 2c), the eccentricities are very high, up to e = 1. With increasing distance to the white region, the eccentricities decrease. Around the host, there is an approximately round area, where the particles are still on circular orbits, i.e. the eccentricities are close to 0.
Figure 2d shows the semimajor axes of the particles that are finally bound to the host or perturber relative to the radius of the initial particle orbit around the host. The particles from the green area (which indicates values around 1) above the negative xaxis have final semimajor axes which are close to the initial ones. Only particles from a small region close to pericentre (yellow to red region) approximately around (0.8, −0.6) finally have larger semimajor axes than initially, these particles are transported outwards. Particles from all other regions, especially along the cleared area and those particles captured by the perturber, finally have smaller semimajor axes, so they are moved inwards (cyan to blue regions). For the captured particles, this agrees with the finding of Pfalzner et al. (2005a) that the captured particles are usually on very tight orbits.
For the massratio shown, we can generalise that, as expected, the closer the VPP of a particle to the perturber’s interaction orbit (see above), the stronger it is influenced: particles with VPPs very close to the interaction orbit become unbound, particles with VPPs far away from the interaction orbit remain bound to the host with relatively low eccentricities.
Since the particles initially rotated counterclockwise, the cleared area (red area in Fig. 2a) mostly lags behind the interaction orbit. Particles which move ahead^{7} of the perturber are decelerated by its gravitational attraction. These particles lose part of their energy and/or angular momentum, and remain finally often on an eccentric orbit, bound more strongly to the host than before. In contrast, particles which move behind the perturber are accelerated. If the acceleration is strong enough they gain enough energy to reach escape velocity and become unbound. The rear (in direction of rotation) borders of the red region in Fig. 2a therefore mark those particles that gain enough energy to exactly reach escape velocity. The particles beyond these lines gain less energy and remain bound to the host on orbits of decreasing eccentricity with increasing distance from the lines. Since particles ahead of the interaction orbit are decelerated and remain bound to the host, while particles behind are accelerated and may become unbound, the perturber’s interaction orbit describes approximately the front border (in direction of rotation) of the cleared region.
3.2. Highmass perturbation
We now show the extreme case of a perturber with a relative mass of m = 20. As can be seen in Fig. 3, the results differ considerably from the previous case since most of the material from the investigated parameter space is now removed, either by becoming unbound or by being captured by the perturber. Only particles with VPPs in relatively small regions remain bound to the host (blue regions in Fig. 3a). These regions are mainly around the host and along the negative and positive xaxis. Interestingly, only very few particles with VPPs between r_{vpp}/r_{p,peri} ≈ 0.3 and r_{vpp}/r_{p,peri} ≈ 2 remain bound to the host. Most of the particles from the upper half of the parameter space become unbound, from the lower half most particles are captured by the perturber.
Figure 3b shows that the interaction orbit (see Sect. 3.1) does not make much sense in this case. This is a consequence of different dominating physical processes for different massratios.
Fig. 4 Relative orbits in the case of a highmass perturber. The centre of mass is indicated by the cross. The solid line depicts the orbit of the host and the dashed line the orbit or the perturber. The grey circle denotes the lowmass material and the arrows show the rotation directions. 

Open with DEXTER 
In the case of a highmass perturber (m ≫ 1), the perturber rests close to the centre of mass, while the host and the particles move on a wide orbit. The distance between host and perturber, and the relative velocity of the host is thereby defined by the parabolic orbit. In the case of a prograde encounter, the particles outside the host’s orbit (see also Fig. 4) have a higher velocity relative to the perturber owing to their superimposed rotation around the host. Since the perturber’s attraction is weaker outside the host’s orbit, these particles become unbound. This is the case for the particles in the upper half of your parameter space. The removal of these particles can be described as centrifugal clearing. The particles inside the host’s orbit have lower velocities relative to the perturber than necessary to overcome the perturber’s gravitational attraction at their positions. These particles mainly become captured by the perturber.
Just as the interaction orbit helped to explain the final properties of the particles in the lowmass perturbation case, this process helps to explain them for the highmass case. Especially, it explains why the particles from the upper half of our parameter space are differently affected during the perturbation than the particles from the lower half.
As in the lowmass case (see Sect. 3.1), particles from an approximately circular region around the host finally have relatively low eccentricities (e ≈ 0, see Fig. 3c). In contrast, most of the particles that are captured by the perturber have very high eccentricities.
Nearly all particles that remain bound to the host finally have semimajor axes that are significantly smaller than the radii of their initial orbits around the host (see Fig. 3d). Only particles which almost become captured by the perturber or become unbound remain bound to the host with larger final semimajor axes (e.g. the red region at ≈(3, 0.8) to (5, 0.8) in Fig. 3d). Also the particles captured by the perturber have, in general, final semimajor axes which are smaller than the radii of their initial orbits around the host. The particles from the green region around (4, 1.5) in Fig. 3a are the exception.
A signature of the centrifugal clearing can also be seen in the final periapsides with the centre of mass of the unbound particles (see Fig. A.1d). The particles that are subject to the centrifugal clearing are ejected roughly in the direction they were moving shortly before periastron passage, which is from right to left (see also Fig. 4). Their hyperbolic orbits basically originate from the positions where they were at that time. Therefore, the periapsides of these particles can approximately be explained with the initial radii of their orbits around the host plus the distance between host and perturber at the moment of ejection, which is ≈1. For example, a particle from (0, 1) has an initial orbital radius of 1 and a final periapsis of ≈2.
3.3. Equalmass perturbation
Fig. 5 As in Fig. 2, but for a massratio of m = 1.0. 

Open with DEXTER 
Comparing the results for the equalmass case (m = 1) with the lowmass and the highmass case, one can see (in Fig. 5) that the equalmass case shows features of both other cases. As in the lowmass case, particles from a band along the approaching branch of the perturber orbit become unbound. Owing to the higher perturber mass the band is broader. The rear border of this area again depicts the region where particles exactly gain escape velocity and the front border can also be described approximately with the perturber’s interaction orbit (see Fig. 5b).
The red, dropshaped region around (−2, 1) in Fig. 5a is centrifugally cleared as in the highmass case. Because of the lower perturber mass, the region is smaller. All other features, especially the red and green regions in the right half of Fig. 5a, are caused by repeated interactions of the particles with host and perturber.
In contrast to the lowmass case, the final eccentricities and semimajor axes of the particles (see Figs. 5c and d) do not show any symmetry. However, the finding that close to the cleared regions, the final eccentricities are very high (e ≈ 1), still holds. The nearly unperturbed region around the host is smaller than in the lowmass case and bigger than in the highmass case.
Most particles that remain bound to the host or are captured by the perturber have final semimajor axes, which are smaller than the radii of their initial orbits. Only particles with VPPs in the red regions around (−4, 2) and (−1.5, 2.5) in Fig. 5d are finally bound to the host with relatively large semimajor axes. These regions are continuations of the region from where the particles are centrifugally removed. Here, the energy gain is just not enough to leave the gravitational field of the host.
Although an encounter between stars of equal mass was often taken as a standard case to investigate, for example, the effect on disc properties (e.g. Clarke & Pringle 1993; Hall et al. 1996), the fate and final properties of the particles are not easy to explain. This is because this case is rather a complicated mixture of low and highmass encounters.
3.4. Final orbital elements of the bound particles
Fig. 6 Final eccentricities of particles in the central region depending on the VPPs for the massratios m = 0.1 a), m = 1.0 b), and m = 20.0 c). 

Open with DEXTER 
Fig. 7 Sample application of the particle fates for nonviscous, lowmass discs with different initial sizes perturbed by a perturber with a massratio of m = 1.0 and different pericentre distances. The scale is the same as in Fig. 5 for better comparison. For full description, see text. 

Open with DEXTER 
In the Figs. 2, 3, and 5, it can be seen that there are large regions in the investigated parameter space where particles have the same fate. Furthermore, similar initial conditions within each of these regions usually result in similar final orbits, indicated by the smooth colour gradients in the images with the final orbital elements of the particles (e.g. Figs. 2c and d). This knowledge can be used to interpolate the final orbital elements for initial conditions between the numerically integrated points.
For the inner parts of the investigated parameter space, where all particles remain bound to the host, the final orbital elements of the particles follow relatively smooth functions of the radius and the angle, θ_{vpp}. As an example, in Fig. 6, the final eccentricities are shown for some sample radii (r_{vpp}/r_{p,peri}) as a function of θ_{vpp}. In contrast to the data for Figs 2, 3, and 5, these data were obtained from simulations where the parameter space was sampled in polar coordinates. The data are not smoothed.
For a fixed θ_{vpp}, the final eccentricities are higher for larger initial radii, as expected. Around θ_{vpp} ≈ 0, the final eccentricities are increased relative to the values outside −90 ≲ θ_{vpp} ≲ 90. The values seem to be a simple function of θ_{vpp} and r_{init}/r_{p,peri}. But the derivation of this dependency lies beyond the scope of this paper.
4. Discussion
We have presented a method to obtain general numerical solutions for the perturbation of lowmass objects on circular Keplerian orbits around a massive object by the flyby of another massive object. As for all numerical solutions, the quality of the results depends on the chosen spacial and temporal resolution.
The temporal resolution is limited by the used integrator. Here, the LSODA integrator is a wellestablished integrator with errorcontrolled timestep size. As the maximum relative error of the integrated values, the builtin default of ≈1.5 × 10^{8} was used. Since the trajectory of each particle was integrated individually, the obtained trajectories for most particles are much more precise than, for example, those of Breslau et al. (2014). Only for some particles that approach one of the massive objects so close that very small time steps are required, the accuracy was deemed insufficient. Therefore, particles that approach one of the two massive objects closer than 0.0001·r_{p,peri} have been excluded from the shown images (see Sect. 3.1).
Other potentially critical factors for the accuracy of the results are the duration of the simulation before (t_{peri}−t_{init}) and after (t_{end}−t_{peri}) pericentre passage of the perturber. By varying the initial force influence of the perturber (see also Sect. 2.2) and hence t_{peri}−t_{init}, we investigated the influence of the starting time on the results. We found an initial relative force influence of F_{p}/F_{h} ≲ 10^{4} to be sufficient for our purpose.
The results also converge for t_{end}−t_{peri} → ∞. By integrating the particle trajectories individually until the orbital elements do not change more than 1 % during a certain time span (see also Sect. 2.2), we ensure that our results are also converged with respect to t_{end}.
For cases where the orbits of the lowmass particles are not exclusively determined by the masses of host and perturber, e.g. due to selfgravity and/or viscosity, our results can not directly be applied. Since, in these cases, the initial angular velocity of the particles no longer only depends on the distance to and the mass of the host, the particle’s VPPs cannot be determined analytically. However, as long as additional forces are small, compared to the gravitational force of the perturber around pericentre, it is still possible to make some predictions since the physical processes will be similar.
For example, for selfgravitating, nonviscous discs, it can be expected that a lowmass perturber also clears only a narrow band along its orbit, depending on the range of its force. The eccentricities of the remaining particles will also decline with increasing distance to the perturber’s interaction orbit. Similarly, even selfgravity and viscosity can only weaken but not prevent the centrifugal clearing of the upper half of our parameter space in the case of a highmass perturber (see also Sect. 3.2). Especially viscosity does not provide a strong enough additional attracting force. The repulsing effect of viscosity has also no impact here since the material is removed in a direction without other material.
In the case of viscous discs, the gradual exchange of energy and angular momentum leads to processes in the disc like accretion and discspreading, even without a perturbation. Therefore, predictions for viscous discs from our results may only be possible for a period after the perturbation that is shorter than the viscous timescale. Depending on the strength of viscosity, during this period the orbital elements of the material may deviate considerably from our results since the influence of the perturber is not yet negligible and the particles have not yet settled into stable orbits. After this short period, the viscous evolution leads to further change in the orbital elements, which makes a comparison difficult. Additionally, predictions from our results may only be possible for material that does not move through regions with higher density during and after the perturbation.
5. Application to planetary systems and nonviscous discs
The results of the restricted threebody approximation presented here are directly applicable to stars surrounded by lowmass discs, either protoplanetary or debris, or planetary systems that are disturbed by the flyby of a second star. However, this approach is only valid if viscous forces can be neglected. Furthermore, the matter has to be on circular orbits before the perturbation.
For the application to discs, we note that the presented datasets contain information what happens to disc material for all encounters from penetrating (with r_{p,peri} = 0.2 r_{disc}⇔ r_{disc}/r_{p,peri} = 5) up to very distant encounters. As already pointed out by Toomre & Toomre (1972), the results for different initial disc sizes are contained in our dataset as subsets. When investigating the effect of an encounter on a disc with a certain size, only the radial subset which corresponds to the relative size of the respective disc is considered.
If, for example, a disc with an initial outer radius of 300 AU is perturbed by a star with a periastron of 100 AU, only the parameter space up to r_{init}/r_{p,peri} = 300 AU/100 AU = 3 would be considered. For the equal mass case, this is illustrated in Fig. 7a, which shows only the inner area of Fig. 5a. Figure 7b shows the same results for a disc with an initial outer radius of 150 AU perturbed by a star with a periastron of 50 AU. Since the whole problem scales with size, the result is exactly the same as for the previous case. Figure 7c shows the results for a disc with an initial outer radius of 200 AU perturbed by a star with a periastron of 100 AU. As can be seen, the result is the inner part of Fig. 7a.
The above example also illustrates that our method of representing the results of a perturbation is independent of the initial disc size as a parameter. Previous investigations often performed simulations of perturbations of discs with different initial sizes (e.g. Pfalzner et al. 2005b). But owing to the inherent geometrical scaling of the problem (see also Sect. 2), the obtained information was redundant. Normalising the initial conditions with the perturber’s pericentre distance, as already performed by, for example, Clarke & Pringle (1993), removes this redundancy. For the simulations for a parameter study, this is a reduction by one free parameter – the initial disc size. When investigating the effect of an encounter on a certain disc, this parameter must be introduced again. But this is much faster than performing one simulation for each disc size.
Our method uses test particles to describe the effect of a flyby on objects with VPPs in a given area. How much matter resides initially in this area depends on the mass distribution within the disc. To investigate the change in a physical quantity, like mass, angular momentum, or energy, one therefore needs to attribute masses to the particles.
Furthermore, a real physical disc also has a vertical mass/particle distribution with a certain scale height, which has to be considered. The material, which initially moves on orbits that are inclined relative to the plane of the perturber orbit, is affected differently than in the coplanar case we present here. The resulting disc would then be a combination of the results for the coplanar case and the respective inclined cases. Typical scale heights of these types of discs correspond to orbital inclinations ≲10°. Tests have shown that for these inclinations the influence is of minor importance.
In contrast to the application to discs, the postencounter properties of a planet can be directly obtained from the maps by scaling the initial radius of the planetary orbit to the perturber’s periastron. For example, the final properties of a planet initially on an orbit with r = 100 AU around a star of mass M_{h} = 1M_{⊙} perturbed by another star with mass M_{p} = 1M_{⊙} and a periastron of r_{p,peri} = 200 AU can be found in the maps at r_{vpp}/r_{p,peri} = 0.5. Since, in practice, θ_{vpp} may be unknown, it might be sufficient to obtain the probability for a certain result depending on r_{vpp}/r_{p,peri}. This method was already applied by Clarke & Pringle (1993; see their Fig. 3), but our results are higher resolved and show less statistical scatter (see Fig. 8). Similarly to the data for Fig. 6, the data for Fig. 8 were obtained from simulations where the parameter space was sampled in polar coordinates.
Fig. 8 Probabilities for certain fates depending on the relative radius of a particle’s initial orbit (r_{init}/r_{p,peri}) for the massratio m = 1.0. See text for full description. 

Open with DEXTER 
6. Summary
In this work, we presented an efficient way to generalise the results of a perturbation of lowmass bodies initially on circular Keplerian orbits around a massive object by the flyby of another massive object. In contrast to previous investigations, our new parametrisation maintains not only the radial resolution for the particles, but also the angular resolution. By also exploiting the other known scaling relations of the problem, the initial conditions for the lowmass bodies for one specific perturber orbit can now be specified with just two parameters. We have created maps that relate the orbital properties of the lowmass bodies after the perturbation to this two dimensional initial condition space. These maps enable us to see the correlation between the particle’s initial position and its final properties. This provides a better understanding of the underlying dynamics.
Maps of the particle fates (whether they remain bound to the host, are captured by the perturber, or become unbound) and some final orbital elements of perturbed particle orbits (i.e. the final eccentricities and semimajor axes) were presented for three sample massratios. With these maps, a detailed comparison of the perturbations with different massratios is possible. The maps show that a lowmass perturber (M_{p}/M_{h}< 1) basically clears a narrow band along its orbit through the parameter space. The closer to the perturber orbit, the higher the final eccentricities and the lower the semimajor axes of the particles that remain bound to the host. For a highmass perturber (M_{p}/M_{h}> 1), the effect is much more violent – only material from small regions remains finally bound. In this case, the removal of particles from the host is a consequence of the centrifugal force acting on the particles during the host’s pericentre passage on its orbit around the heavier perturber. The equalmass case (M_{p}/M_{h} = 1) is a mixture of the lowmass and the highmass case.
Regions in the presented maps, where particles have the same fate, show smooth gradients in the final orbital elements. As a consequence, the presented numerical results are not only valid for exactly the integrated initial conditions, but can be generalised for the full presented parameter space by interpolation.
The results are directly applicable to astrophysical processes like the perturbation of planetary orbits or protoplanetary discs by passingby stars without the need for additional simulations. The fates of planets in an encounter and the final orbital elements can directly be obtained from the maps by scaling the initial planetary orbit with the perturber’s periastron distance. When applied to discs, our results are independent of the initial disc size, in contrast to most previous works. The result of the perturbation of a disc with finite size can be obtained from our maps by considering only the respective radial subset. For the simulations of a parameter study in this context, this reduction by one parameter reduces the numerical effort significantly.
Even in the case of selfgravitating and viscous discs, the results may be applicable to some extent. Material that is removed from the disc without moving through denser regions, as in the case of a highmass perturber, will most likely behave in a similar way. By way of contrast, in the case of crossing particle trajectories, viscosity will change the results considerably.
Indeed, this fact is used in all numerical investigations of gravitational perturbations to restrict the simulation time (e.g. Hall et al. 1996, start their investigation of the perturbation of a disc with an initial distance of at least 10 times the disc radius).
The mathematical bases can be found for example in Bate et al. (1971) or Appendix 1.5 of Kemble (2006).
As the orbital elements are calculated from the final positions of the particles relative to the centre of mass and the centre of mass moves relative to host and perturber, it is difficult to transfer this information into the VPP system, which is hostcentred. A small periapsis distance does not mean that a particle’s hyperbolic trajectory originates from close to the host. It just means that this trajectory has a small periapsis distance with the centre of mass of host and perturber. Because the orbital elements are calculated at the end of the integration of the particle’s trajectory, when the distance between host and perturber is usually large again, the position of the systems centre of mass may differ considerably from the position of the host.
Acknowledgments
We thank the anonymous referee for very useful comments that helped to improve this paper significantly. A.B. also wants to thank A. Bhandare, S. Brackertz, T. Bardenheuer, M. Engler, C. Korntreff, S. Libi, K.H. Rehren, M. Sonntag, A. Wolff, and M. XiangGruess for fruitful discussions. Additionally, A.B. is indebted to his family, without whose support this work might not have been finished.
References
 Bate, R. R., Mueller, D. D., & White, J. E. 1971, Fundamentals of astrodynamics (New York: Dover Publication) [Google Scholar]
 Breslau, A., Steinhausen, M., Vincke, K., & Pfalzner, S. 2014, A&A, 565, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Clarke, C. J., & Pringle, J. E. 1993, MNRAS, 261, 190 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, J., & Krumholz, M. R. 2013, ApJ, 769, 150 [NASA ADS] [CrossRef] [Google Scholar]
 D’Onghia, E., Vogelsberger, M., FaucherGiguere, C.A., & Hernquist, L. 2010, ApJ, 725, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Farouki, R., & Shapiro, S. L. 1981, ApJ, 243, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Fregeau, J. M., Chatterjee, S., & Rasio, F. A. 2006, ApJ, 640, 1086 [NASA ADS] [CrossRef] [Google Scholar]
 Hall, S. M. 1997, MNRAS, 287, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Hall, S. M., Clarke, C. J., & Pringle, J. E. 1996, MNRAS, 278, 303 [NASA ADS] [Google Scholar]
 Hao, W., Kouwenhoven, M. B. N., & Spurzem, R. 2013, MNRAS, 433, 867 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D. C., & Hut, P. 1993, ApJS, 85, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, C. H. 1995, ApJ, 455, 252 [NASA ADS] [CrossRef] [Google Scholar]
 Hindmarsh, A. C. 1983, IMACS Transactions on Scientific Computation, 1, 55 [Google Scholar]
 Hut, P. 1983, ApJ, 268, 342 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P., & Bahcall, J. N. 1983, ApJ, 268, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Kemble, S. 2006, Interplanetary Mission Analysis and Design (Berlin: Springer) [Google Scholar]
 Kobayashi, H., & Ida, S. 2001, Icarus, 153, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Larwood, J. D. 1997, MNRAS, 290, 490 [NASA ADS] [Google Scholar]
 Levison, H. F., Morbidelli, A., & Dones, L. 2004, AJ, 128, 2553 [NASA ADS] [CrossRef] [Google Scholar]
 Li, G., & Adams, F. C. 2015, MNRAS, 448, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Malmberg, D., Davies, M. B., & Heggie, D. C. 2011, MNRAS, 411, 859 [NASA ADS] [CrossRef] [Google Scholar]
 Marchal, C. 1990, The threebody problem (Amsterdam: Elsevier) [Google Scholar]
 Melita, M. D., Larwood, J., CollanderBrown, S., et al. 2002, in Asteroids, Comets, and Meteors: ACM 2002, ed. B. Warmbein, ESA SP, 500, 305 [Google Scholar]
 Musielak, Z. E., & Quarles, B. 2014, Rep. Prog. Phys., 77, 065901 [NASA ADS] [CrossRef] [Google Scholar]
 Olczak, C., Pfalzner, S., & Spurzem, R. 2006, ApJ, 642, 1140 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, E. C. 1994, ApJ, 424, 292 [NASA ADS] [CrossRef] [Google Scholar]
 Pfalzner, S., Umbreit, S., & Henning, T. 2005a, ApJ, 629, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Pfalzner, S., Vogel, P., Scharwächter, J., & Olczak, C. 2005b, A&A, 437, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Punzo, D., CapuzzoDolcetta, R., & Portegies Zwart, S. 2014, MNRAS, 444, 2808 [NASA ADS] [CrossRef] [Google Scholar]
 Scharwächter, J., Eckart, A., Pfalzner, S., Saviane, I., & Zuther, J. 2007, A&A, 469, 913 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Steinhausen, M., Olczak, C., & Pfalzner, S. 2012, A&A, 538, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Toomre, A., & Toomre, J. 1972, ApJ, 178, 623 [NASA ADS] [CrossRef] [Google Scholar]
 Valtonen, M., & Karttunen, H. 2006, The ThreeBody Problem (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Valtonen, M. J., & Heggie, D. C. 1979, Celest. Mech., 19, 53 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Properties of the unbound particles
In this section, we show the results for the particles that become unbound during the perturbation. Astrophysical investigations of the perturbation of protoplanetary discs or planetary systems usually focus on matter that remains bound to the host or is captured by the perturber. Therefore, these results are probably of minor importance. However, the final orbital elements of the unbound particles support our explanation for the two different processes for low and highmass perturbations.
For the unbound particles, we show the orbital elements with the centre of mass of host and perturber. The unbound particles are on hyperbolic orbits with the centre of mass and have, therefore, negative semimajor axes. Because these are difficult to interpret, we show instead the periapsis distances^{8} of those particles.
Appendix A.1: Lowmass perturbation
For the massratio of m = 0.1, the final eccentricities of the unbound particles are shown in the map in a logarithmic colour scale up to e = 10 (see Fig. A.1a). The particles from the band along the approaching branch of the perturber’s orbit finally have high eccentricities, while the particles from the band along the departing branch have eccentricities slightly above 1.
Figure A.1b shows the final periapsides of the unbound particles. The particles from the band along the approaching branch of the perturber’s orbit have finally periapsis distances with the centre of mass comparable to their initial orbital radii. In contrast, the particles from the departing branch of the perturber’s orbit have relatively small final periapsis distances. By also taking into consideration the eccentricities of the final orbits of these particles, we can say that particles with VPPs close to the approaching branch are ejected on hyperbolic orbits with relatively high periapsis distances. Particles from the departing branch become barely unbound with eccentricities just above 1 and small periapsides. For these barely hyperbolic orbits, small periapsides correspond to small semimajor axes and thus to relatively small angular momenta.
Appendix A.2: Highmass perturbation
For the massratio of m = 20, the particles that become unbound finally have remarkably low eccentricities, most have e< 2 (see Fig. A.1c). Only particles with VPPs in the upper half, close to the host (cyan region), have final eccentricities ≈3. Close to the lower border of the image there is a region from where the particles are ejected on orbits with higher eccentricities (rainbow around e.g. (0, −4.2)). Compared to the lowmass case, more particles become unbound, but on significantly less hyperbolic orbits.
The gradient in the upper half of the final periapsides of the unbound particles (see Fig. A.1d) is caused by the centrifugal clearing (see also Sect. 3.2 for explanation). Only, from the lower half do a few particles become unbound, many of them on barely hyperbolic orbits with small periapsides.
Appendix A.3: Equalmass perturbations
The similarities with the low and the highmass case, as described in Sect. 3.3, can also be seen in the final periapsides with the centre of mass of the unbound particles (see Figs. A.1e and f). In the final eccentricities and periapsides of the particles along the approaching branch of the perturber’s interaction orbit, similarities to the lowmass case can be seen. The dropshaped, centrifugally cleared region shows similarities with the highmass case. As in the lowmass case, along the approaching branch of the perturber orbit an increase of the final periapsides with increasing distance of the VPP to the host can be seen. A similar increase can be seen in the centrifugally cleared region. All other unbound particles finally have relatively small periapsides with the centre of mass. The final eccentricities are again surprisingly low; they are all below ≈3.
Fig. A.1 Maps of the final eccentricities (left column) and periapsides (right column), with the centre of mass of the unbound particles for the three massratios. 

Open with DEXTER 
All Figures
Fig. 1 Virtual pericentre positions (VPPs) for a perturbed disc. The perturber has the same mass as the host and the pericentre of the parabolic, prograde orbit is at (1, 0). In a) and b) the disc is seen face on when the perturber has nearly reached pericentre and at t = t_{init}, respectively. In c) the VPPs of the disc particles are shown. For full description see text. 

Open with DEXTER  
In the text 
Fig. 2 Maps of the final properties of the particles depending on their VPPs for the region r_{vpp}/r_{p,peri}< 5 for a massration of m = 0.1. a) and b) show whether the particles are finally bound to host or perturber, or become unbound. a) Also shows the orbit of the perturber through the Cartesian space and b) the perturber’s interaction orbit (see text for explanation). c) Shows the final eccentricities of the bound particles, and d) shows the final semimajor axes of the bound particles relative to the radii of their initial orbits. 

Open with DEXTER  
In the text 
Fig. 3 As for Fig. 2, but for a massratio of m = 20.0. 

Open with DEXTER  
In the text 
Fig. 4 Relative orbits in the case of a highmass perturber. The centre of mass is indicated by the cross. The solid line depicts the orbit of the host and the dashed line the orbit or the perturber. The grey circle denotes the lowmass material and the arrows show the rotation directions. 

Open with DEXTER  
In the text 
Fig. 5 As in Fig. 2, but for a massratio of m = 1.0. 

Open with DEXTER  
In the text 
Fig. 6 Final eccentricities of particles in the central region depending on the VPPs for the massratios m = 0.1 a), m = 1.0 b), and m = 20.0 c). 

Open with DEXTER  
In the text 
Fig. 7 Sample application of the particle fates for nonviscous, lowmass discs with different initial sizes perturbed by a perturber with a massratio of m = 1.0 and different pericentre distances. The scale is the same as in Fig. 5 for better comparison. For full description, see text. 

Open with DEXTER  
In the text 
Fig. 8 Probabilities for certain fates depending on the relative radius of a particle’s initial orbit (r_{init}/r_{p,peri}) for the massratio m = 1.0. See text for full description. 

Open with DEXTER  
In the text 
Fig. A.1 Maps of the final eccentricities (left column) and periapsides (right column), with the centre of mass of the unbound particles for the three massratios. 

Open with DEXTER  
In the text 