Discprotoplanet interaction
Influence of circumprimary radiative discs on selfgravitating protoplanetary bodies in binary star systems
^{1} Institute for Astrophysics (IfA), University of Vienna, Türkenschanzstrasse 17, 1180 Vienna, Austria
email: markus.gyergyovits@univie.ac.at
^{2} IMCCE – Observatoire de Paris, 77 avenue DenfertRochereau, 75014 Paris, France
^{3} Planetarium Mannheim, WilhelmVarnholtAllee 1, 68165 Mannheim, Germany
^{4} Department of Physics, University of Graz, Universitätsplatz 5, 8010 Graz, Austria
Received: 7 May 2013
Accepted: 5 May 2014
Context. More than 60 planets have been discovered so far in systems that harbour two stars, some of which have binary semimajor axes as small as 20 au. It is well known that the formation of planets in such systems is strongly influenced by the stellar components, since the protoplanetary disc and the particles within are exposed to the gravitational influence of the binary. However, the question on how selfgravitating protoplanetary bodies affect the evolution of a radiative, circumprimary disc is still open.
Aims. We present our 2D hydrodynamical GPUCPU code and study the interaction of several thousands of selfgravitating particles with a viscous and radiative circumprimary disc within a binary star system. To our knowledge this program is the only one at the moment that is capable to handle this many particles and to calculate their influence on each other and on the disc.
Methods. We performed hydrodynamical simulations of a circumstellar disc assuming the binary system to be coplanar. Our gridbased staggered mesh code relies on ideas from ZEUS2D, where we implemented the FARGO algorithm and an additional energy equation for the radiative cooling according to opacity tables. To treat particle motion we used a parallelised version of the precise Bulirsch – Stoer algorithm. Four models in total where computed taking into account (i) only Nbody interaction; (ii) Nbody and disc interaction; (iii) the influence of computational parameters (especially smoothing) on Nbody interaction; and (iv) the influence of a quiet loweccentricity disc while running model (ii). The impact velocities were measured at two different time intervals and were compared.
Results. We show that the combination of disc and Nbody selfgravity can have a significant influence on the orbit evolution of roughly Moon sized protoplanets.
Conclusions. Not only gas drag can alter the orbit of particles, but the gravitational influence of the disc can accomplish this as well. The results depend strongly on the state of the disc (i.e. quiet or dynamically evolving) – according to encounterprobability distributions, planet formation can be strongly altered if there is a dynamically evolving gas disc – and also on the smoothing parameter.
Key words: accretion, accretion disks / hydrodynamics / protoplanetary disks / binaries: close
© ESO, 2014
1. Introduction
Most mainsequence stars in the solar neighbourhood are members of multiple systems (Duquennoy & Mayor 1991), and the distribution of their separations peaks at about 50 au (Eggleton 2006). This suggests a significant influence of the secondary star on the circumprimary protoplanetary disc. Planets have been discovered in midseparation binaries, with semimajor axes of about 20 au, for example, Gl86 (Santos et al. 2000; Lagrange et al. 2006), HD 41004 (Zucker et al. 2004), γ Cephei (Hatzes et al. 2003; Neuhäuser et al. 2007), and the recently discovered HD 196885 (Chauvin et al. 2011) and α Centauri B (Dumusque et al. 2012), which has led to a growing interest in understanding planetary formation processes in general. In midseparation binary star systems, the secondary and the heavily truncated and possibly distorted disc causes the formation and evolution of planets throughout several stages of the planetforming process. The protoplanetary disc is truncated by the companion star through gravitational interaction, as was shown by Artymowicz & Lubow (1994) and Savonije et al. (1994). This shortens the viscous lifetime of the disc and consequently limits the period in which gasgiants can form. During accretion phases to kmsized particles, the second star can force the shape and orientation of planetesimal orbits, which can stop planetary formation (e.g. Thebault 2011, and references therein).
Despite the considerable progress this field has seen over the past decades (Marzari & Scholl 2000; Barbieri et al. 2002; Lissauer et al. 2004; Thébault et al. 2004, 2006; Quintana et al. 2007; Haghighipour & Raymond 2007; Guedes et al. 2008; Thebault 2011), many open questions remain, especially with regard to accretion stages where kmsized particles are expected to grow into planetary embryos within a gaseous disc. This process demands fast runaway and oligarchic growth phases, which in turn require low encounter and impact velocities (dv) between protoplanetary bodies. Previous studies (Marzari & Scholl 2000; Thébault et al. 2004, 2006) have shown that coupling the protoplanetary bodies to gas via gas drag creates a sizedependent phasing of protoplanetary orbits that decreases dv for equalsized objects but increases it for all other types of impacts (Thebault 2011). Introducing small inclinations between the static, axisymmteric, circumprimary gas disc and the binary orbital plane might facilitate accretion (Xie & Zhou 2009; Xie et al. 2010). Recently, Kaib et al. (2013) showed that wide binaries can also influence the evolution of embedded planetary systems even Gyr after their formation. These studies are based on either pure Nbody simulations that neglect the influence of a gas disc (Marzari & Scholl 2000; Barbieri et al. 2002; Lissauer et al. 2004; Quintana et al. 2007; Haghighipour & Raymond 2007; Guedes et al. 2008) or only take gas drag caused by a static circular disc (Thébault et al. 2004, 2006; Thebault 2011) into account. A static circular disc, however, is not realistic for close binary star systems since the secondary exerts spiral density waves within the disc, which will presumably alter the evolution of protoplanets by either drag force or gravitational force. In an eccentric binary system even the assumption of circularity does not hold any more. This was shown by Kley & Nelson (2008) and Kley et al. (2008).
The influence of an evolving gasdisc on protoplanetary accretion has been studied first by Paardekooper et al. (2008) in two dimensions and more recently by Fragner et al. (2011) in 3D. Both studies used a locally isothermal disc as underlying model without taking the selfgravity of particles into account. Müller & Kley (2012) showed, that radiative discs have far lower eccentricities than their locally isothermal counterparts. This means that, changing the radiative properties of the protoplanetary disc can lead to a dynamically quiet environment more suitable for planet formation than the local isothermal approach.
While most of the studies described in the previous paragraph focused on embryo formation via protoplanetary accretion, little attention has been given so far to the penultimate stage in the formation scenario, where relatively large protoplanets interact gravitationally with each other as well as with a dynamically evolving disc. The lack of such investigations into this regime might be related to the huge computational requirements for simulations where selfgravitation of the protoplanets and the gravitational interaction with the disc is more important than gas drag.
We here introduce a GPUCPU 2D hydrodynamic grid code that combines hydrodynamic radiative disc computations with highly accurate Nbody simulations. The development of this new code allows us to investigate the interplay between a circumprimary gas disc and the orbital evolution of embedded selfgravitating protoplanets on timescales of one hundred binary periods. The application of our code using various models shows the differences for planetary formation, when taking into account (i)binary – disc interactions; (ii)binary – protoplanet interaction; (iii) binary – protoplanet – disc interactions. In Sect. 2 we describe important features of our code and present results for standard tests in Appendix A. The results from our numerical investigations of combined discparticle systems are presented in Sect. 3.1. A summary of our findings and concluding remarks are given in Sect. 5.
2. Code description
Our code is based on the ZEUS2D solver for hydrodynamical equations on a twodimensional polar grid (Stone & Norman 1992). It was originally developed by Korchagin & Theis (1999) for simulations of galaxies and computes interactions of up to three different fluids in one simulation (Theis & Orlova 2004). The heart of the code consists of a staggeredmesh finitedifference method using operator splitting and a firstorder integrator in time to handle the source step. The advective terms are solved by a secondorder upwind algorithm in time – for a detailed description see Stone & Norman (1992). The adaptation of the code to hydrodynamical simulations of gaseous discs in tight binary systems required the following major modifications: (i) an accelerated coordinate system that is centred on the primary star; (ii) implementing of the socalled FARGO acceleration algorithm (Masset 2000) to speed up calculations; (iii) we partly rewrote the code to run on NVIDIA graphics cards in doubleprecision accuracy to run simulations that include the full gravitational interaction of several thousand particles with the disc.
The equations of motion of the hydrodynamic part of the code in an accelerated coordinate system are given by
where v_{r}, v_{φ} are the velocity components in r and φ direction. P denotes the pressure, Σ(r,φ) is the vertically integrated density, Ψ(r) is the gravitational potential exerted by all bodies (a detailed description is given below in this section) and f_{r}, f_{φ} denote the accelerations due to the noninertial frame. Here, the secondary as well as the massive particles moving in the disc will contribute to this acceleration. The effects of viscosity are treated as discussed in Kley (1999) and Landau & Lifshitz (1959). We implemented the standard viscous stress tensor, but used the turbulent viscosity coefficient ν_{t} instead of the molecular viscosity. This leads to the equations (2)and (3)where T_{rr}, T_{φφ}, T_{rφ} are the components of the viscous stress tensor: (4)▽ ·v denotes the divergence of the velocity and η and ξ are the shear and bulk viscosity, respectively. In our simulations the bulk viscosity was set to zero. For the shear viscosity we used η = Σν_{t}. Depending on the simulation scenario, the turbulent viscosity coefficient was either set constant or according to the α prescription by Shakura & Sunyaev (1973). We found that better results – fewer postshock oscillations in the Sod shocktube test (Sod 1978) – can be achieved when we applied an artificial viscosity of von Neumann & Richtmyer type (Stone & Norman 1992).
The gravitational potential exerted on the disc by the masses (M_{i}) is given by (5)where r_{1}, r_{2} are the position vectors and M_{1}, M_{2} are the masses of the two stars (primary and secondary). The secondary is moving outside the computational grid. Furthermore, G denotes the gravitational constant, r_{i}, M_{i} are the position vectors and masses of the smaller bodies, and r is the point of interest in the disc. For bodies moving in the grid, gravity softening was used via the ϵ softening parameter (de ValBorro et al. 2006, Eq. (3)).
For closure conditions, the code offers (i) the widely known locally isothermal approach; (ii) a globally isothermal model; and (iii) a radiative disc model. To solve the required energy equation we followed Müller & Kley (2012) with one small exception: instead of mixing adiabatic and isothermal sound speed, we used only the latter one. For reasons of numerical stability we also used a lowest value for the optical depth. The vertically integrated energy equation reads (6)where e is the internal energy density, Q_{+} the heating source term, Q_{−} the cooling source term (Müller & Kley 2012). The radiative transport in the plane of the disc is treated in the fluxlimited diffusion approximation. The flux is given by (7)where c denotes the speed of light, a the radiation constant, ρ the midplane density, κ the Rosseland mean opacity, and λ is the fluxlimiter Kley (1989, Eq. (9)) and Kley & Crida (2008). The initial radial velocity v_{r} was set to zero, and v_{φ} is subKeplerian, meaning that the pressure terms as well as the viscosity are taken into account when the initial values for v_{φ} are computed from Eq. (2).
To handle the Nbody calculations accurately we implemented a Bulirsch – Stoer integrator (Bulirsch & Stoer 1964) parallelised to run on the GPU and thus gaining a speed factor of up to 40.
3. Application of the code
3.1. Models
Because the tests showed (see Appendix A) that our code provides satisfactory performance, we applied our new code to study systems with full gravitational interaction between the binary stars, the disc, and the protoplanets. In close binary star systems such as γ Cephei – which we took as our model – it can be assumed that close encounters between protoplanets and additional interaction with frequently occuring spiral waves induced by the secondary star in the gas disc lead to momentum and energy exchange of such an amount that this might increase the semimajor axis and eccentricity of the protoplanets to high values and thus hinder agglomeration and merging of particles and protoplanets.
Müller & Kley (2012) showed that in such systems the massweighted disc eccentricity and longitude of disc pericenter of the gas disc can change their dynamical behaviour from damped oscillation to a constant small one (for the first) or from rotation to oscillation (the second). Such changes in dynamical evolution of the disc might also affect movement of protoplanets and may either render planet formation impossible or provide an environment where agglomeration is not strongly affected by the disc.
Therefore we investigated:

whether the combined gravitational influence (stars,protoplantes and disc) dampens or increases the orbital forcing of the protoplanets; and

to which degree the choice of the smoothing parameter affects these simulations.
For this purpose we ran the following simulations:
reference model: here we included only the binarydiscinteraction.
model a1: several thousand selfgravitating protoplanets distort the disc gravitationally but no backreaction from the disc on the particles is considered. At the same time, the disc and the bodies move under the gravitational influence of the binary star.
model a2: the same initial conditions as model a1 are used except for a higher smoothing parameter when calculating the particle – particle forces.
model b1: the full gravitational interaction between particles, the disc and the binary star is taken into account.
model b2: we recalculated the reference model for 50 000 yr, reset disc mass to its initial value, and then inserted particles in the disc taking into account the full gravitational interaction between particles, the disc and the binary star.
The results of models a1, a2, b1 and b2 were compared with the disc evolution of our reference model and with each other.3.2. Numerical setup
For the simulations we used a twodimensional logarithmic polar grid ranging from 0.5 au to 8 au, and the initial conditions given in Table 1. The masses of the primary and secondary were 1.4 M_{⊙} and 0.4 M_{⊙}. Their semimajor axis was 20 au and the chosen eccentricity was 0.4. These values are close to those of γ Cephei provided by Neuhäuser et al. (2007) and were also used in the study of Müller & Kley (2012). Our initial disc mass was 0.01 M_{⊙} and we used an alpha viscosity prescription, where α has a value of 5 × 10^{3}. The adiabatic index γ was set to 7/5. The initial density and temperature profiles of the gas disc follow a ∝ r^{1} distribution. Our initial disc aspect ratio (H/r) was 0.05, which is a commonly used value. Moreover, we chose different grid resolutions in φ and r directions to allow for locally square cells where Δr ≈ rΔφ. We ran simulations with 2048 selfgravitating protoplanetary bodies and the two stellar components of the binary, where we fixed the masses to values that would follow from a complete coagulation of dust – the dusttogasratio being 10^{2}. This gives particles of mass ~0.016 M_{⊕}. As a first approximation we removed the term connected to the radiative flux in the disc in the energy equation, since its contribution to the whole equation is several orders of magnitude lower than that of the other terms. The total integration time for our reference model as well as the simulation runs of models a1, a2, b1, and b2 was 100 binary periods, corresponding to 6670 years. Initially, the protoplanets were distributed on 32 circular orbits centred on the binary’s primary star. On each of the 32 orbits, 64 particles were evenly distributed so that all particles are started on linearly spaced circular orbits ranging from 1 to ~4 au in radial direction and chosen in such a way that Δr ≈ r_{0}Δφ, where r_{0} = 1 au. Thus the last ring is close to the critical semimajor axis a = 3.920 ± 0.009 au where stable motion is still possible according to the studies by Holman & Wiegert (1999) and PilatLohinger & Dvorak (2002). We chose this quite artificial approach of an ordered state because we did not wish to introduce additional perturbation in the system. We assumed that due to the given gravitational influences the system will adjust itself during the simulation. Of course, one could have tested a variety of different initial conditions and their influence on the outcome of the simulations, but this would be a task of its own and will be presented in a forthcoming article.
Initial conditions for the simulations.
Since these particles are moving within the hydrodynamical grid, we used a gravity softening of ϵ = 0.6H_{p}, a value commonly found in literature. We did not include the more sophisticated treatment suggested by Müller et al. (2012). For the chosen particle masses the corresponding softening sphere is slightly larger than the particles’ Roche lobe with respect to the primary. Using gravity softening for the particle – particle interactions as well (ϵ_{p} = 10^{5} au), we decided not to take protoplanet merging into account. For protoplanets with mean lunar density, the associated softening distance is of the order of their respective radii. Such a weak softening is not expected to influence the computations significantly, but it prevents very small computational timesteps. For this reason processes such as gravitational focusing of the bodies or coagulation of particles are probably not inhibited during the simulation. Given the current choices of initial and boundary conditions, the gas drag exerted by the disc on the protoplanets would be two to three orders of magnitude lower than the selfgravity of the two protoplanets at opposing outer edges of the truncated disc. We therefore neglected possible influences of gas drag and concentrated only on gravitational particle – particle and particle – disc interactions. To compare all our results we kept the secondary moving on a fixed ellipse according to its initial conditions. Reflecting boundary conditions at the inner boundary and zero gradient boundary condition at the outer boundary were applied for all simulations (for a short description on boundary condition details see also Müller & Kley 2012). Modelling the inner boundary is quite a delicate task since we do not have any information about this region. A zerogradient boundary condition will create a (eccentric) hole in the vicinity of this boundary (Kley & Nelson 2008). Taking a nonreflecting boundary condition (Godon 1996), one also has to make assumptions about disc behaviour inside the boundary (which means outside the computational domain). We have chosen this rigid boundary since this initial density distribution needs to adjust itself to ongoing perturbation caused by the secondary and therefore will move inward at first. Of course, this boundary will reflect every wave into the computational regime and cause artificial superpositions of waves. But since the density is very high at the inner regions we expect these reflected waves to be damped very fast. The initial density of the gas disc was modelled according to the following powerlaw distribution: (8)where Σ_{0} is the value of the vertical averaged density at r = 1. To avoid unrealistically low gasdensities we applied a floor value (fl), which is given by fl = Σ_{0} × 10^{7}. The density was modified by the function f(r), which emulates a truncated disc. For this purpose we introduced the Fermi function (9)with r_{t} = 5 au. The radial velocity (v_{r}) of the disc was set to zero and the azimuthal velocity (v_{φ}) was adjusted taking the pressure gradient and the viscosity of the disc into account. To analyse the disc behaviour and compare our work with others we also computed the massweighted eccentricity (e_{mw}) and argument of periastron of the gas disc (ϖ_{mw}) where we followed Müller & Kley (2012), and which are given by and e and ϖ are the eccentricity and argument of pericenter computed for each cell separately. The integrals where evaluated by summation over all grid cells. This definition favours the inner parts of the disc, which are more dense, and thus represents the evolution of the region – at least in these parameters – in which we are interested in.
4. Results and discussion
4.1. Evolution of the disc
Fig. 1 Time evolution of the massweighted eccentricity for a disc without particles (reference model) and model a2 (upper plot), and models a1 and b1 (lower plot). Time is given in terms of binary orbits, where one orbit corresponds to 66.7 yr. See Sect. 3.1 for details. 

Open with DEXTER 
Fig. 2 Time evolution of the massweighted argument of the pericenter for a disc without particles (reference model) and model a2 (upper plot), and models a1 and b1 (lower plot). Time is given in terms of binary orbits, where one orbit corresponds to 66.7 yr. 

Open with DEXTER 
Fig. 3 Time evolution of massweighted argument of the pericenter for a disc without particles (upper plot) and of massweighted eccentricity (lower plot) for 50 000 yr after applying a running window average. Time is given in terms of binary orbits, where each orbit corresponds to 66.7 yr. We find transition from circulation to oscillation within 200 binary orbits for ϖ_{mw} and a subsequent oscillation around ≈0.6 rad. A damped oscillation is visible for e_{mw}, reaching a value around ≈0.0275. 

Open with DEXTER 
In Figs. 1–3 we show the evolution of the massweighted eccentricity and massweighted argument of the pericenter for all models. To eliminate small oscillation caused by the periodical approach of the secondary (orbital period ~66.7 yr) we applied a running window filter. For all of the models the eccentricity of the disc, which initially was zero, was increased to values of nearly 0.06–0.07 and reached ~0.03 (model a1, b1) or ~0.035 (reference model and model a2) at the end of the simulation (Fig. 1). The argument of the pericenter is rotating retrograde at the beginning of the simulation, but as time proceeds, the circulation turns into a damped oscillation around 0.5 radians for models a1 and b1 (Fig. 2, lower plot). For the reference model and model a2 (Fig. 2, upper plot) this effect seems to be present as well, but not as pronounced as in the first case. The results for the reference model agree well with those obtained by Müller & Kley (2012).
The extended integration of our reference model (Fig. 3) shows a similar behaviour as models a1 and b1. We find transition from circulation to oscillation within 200 binary orbits for ϖ_{mw} and a subsequent oscillation around ≈0.6 rad. A damped oscillation is visible for e_{mw}, reaching a value of about ≈0.0275. A similar transitional behaviour was also found in Müller & Kley (2012).
A comparison of Figs. 1 and 2 indicates that a protoplanetary disc seems to damp the evolution of e_{mw} and ϖ_{mw} such that e_{mw} for models a1 and b1 is lower than e_{mw} for the reference model and the amplitude of oscillation of ϖ_{mw} is also decreased. On the other hand, we see that an increasing protoplanetary gravitational softening parameter for Nbody interaction increases the amplitudes for e_{mw} as well as for ϖ_{mw}, indicating that this parameter has a significant influence on the computation.
Using an initially truncated disc not only depicts the response of the disc due to a perturber in a better way, it also diminishes the mass loss. This is because using a simple r^{1} density distribution places much material outside the Hill sphere of the primary where it is strongly affected by the perturbations of the secondary. This material is quickly removed from the disc within a few revolutions of the second star. This behaviour was found by the authors during test runs and can also be seen in Fig. 12 of Müller & Kley (2012) (small jump within the first 1000 years). The mass loss during our simulations with respect to the initial disc mass was 0.68% for the reference model, 0.52% for model a1, 0.75% for model a2, 0.42% for model b1, and 0.08% for model b2.
Since we are using a rigid boundary condition at the inner boundary we can not provide data about mass loss at this boundary. For all models we use an alpha value of α = 0.005.
4.2. Evolution of particles
Fig. 4 Time evolution of the particle semimajor axis influenced by the disc (model b1). We find a transition from ordered motion close to the initial positions on the grid (left lower corner) to a distribution of semimajor axes that spreads across the whole stable region of the disc within 900 yr. 

Open with DEXTER 
Fig. 5 Evolution of protoplanet mean (upper plot) and root mean square eccentricity (lower plot) for all four models. Highest values in both plots are reached by model b, whereas model b1 shows a similar behaviour as model a, and the lowest values are reached by model a1. 

Open with DEXTER 
Because periodical perturbations of the secondary alter the disc, the embedded particle orbits are affected by the changes in the gas disc and the direct perturbations of the secondary. The strong perturbations experienced by the outer disc propagate inwards through particleparticle interactions that cause strong deviations from the protoplanets’ initial orbits. This process takes ~900 years, as shown for model b1 in Fig. 4. We find a similar behaviour for model a1 with an increased evolution time of ~2200 years, for model a2 where the time amounts to ~3500, and for model b2 it takes ~1100 years (all not shown). The eccentricity of the particles undergoes strong variations in model b1, reaching values of up to ~0.9 in the first 1000 years of the simulation. For model a1 the evolution is proceeding more slowly, but on average stays below 0.4. To compare the outcome for the four models quantitatively we show the time evolution of the protoplanets’ mean eccentricity ⟨e⟩ as well as their root mean square (rms) eccentricity in Fig. 5. For models a1 and a2 ⟨e⟩ increases almost linearly from 2000 years onwards, whereas for model b1 ⟨e⟩ rises quickly to 0.3 and oscillates finally between 0.2 and 0.25. For model b2 from 10 orbital periods on we find a linear increase in ⟨e⟩ with superposed oscillations. From 70 to 100 orbits the evolution of ⟨e⟩ is lower than that of model a1. To give an appropriate picture of the state of the system we also computed the rms eccentricity (lower plot of Fig. 5), which is of the order of ⟨e⟩, reflecting the fact that due to the gravitational perturbation of the secondary we find a strong variation of the particle eccentricity within the computational domain. The RMS values of the eccentricities for models a1 and a2 increase linearly from ~15 orbits onwards, reaching values between 0.15 and 0.1 (model a1) and 0.8 and 0.9 (model a2) in the last ten orbits of the simulation. Model b1, in contrast, shows a general increase in the protoplanets’ eccentricities which reaches a value of 0.16 within the first ten orbits but slightly decreases to 0.135 at the end of the run. In model b2 the RMS rises to 0.13 within 20 orbits and decreases to values below 0.1. It again is lower than values of model a1 in the range of 70 to 100 orbits. This shows that the influence of a dynamically evolving gas disc onto the protoplanets causes an increase in particle eccentricities, whereas a quiet disc damps particle eccentricity even when we solely take disc gravity into account.
For the case of a gasfree disc we find that our simulation does not show the wavy forced eccentricity patterns encountered in other studies where the selfgravity of the protoplanets has been neglected (Marzari & Scholl 2000; Thébault et al. 2006; Paardekooper et al. 2008). Such patterns are quickly removed by orbital crossings that start at the outer edge of the system and propagate inside. Because none of the orbit crossing protoplanets was purged from the simulations, the number of potentially high velocity collisions is most likely overestimated and the vanishing intermittent phasing as well the lack of the wavy forced eccentricity patterns may be nothing else than a numerical effect of our applied conditions.
Fig. 6 Collision probabilities for 30–40 orbital periods for model a1 (upper left), model b1 (upper middle), and model b2 (upper right) as well as for 50–60 orbits for model a1 (lower left), model b1 (lower right), and model b2 (lower right). We show the plot radius (in au) versus the probability for an event (in percent) for disruption and merging for methods m1 and m2 (explained in the text). 

Open with DEXTER 
In a second step we performed simulations where we computed differential velocities (dv) for each close encounter of embryos. A close encounter occurs when an embryo enters the Hill sphere of another. For numerical reasons we computed this in two time intervals: (1) from the 30th to 40th binary orbit (I (30–40)); and (2) from the 50th to 60th binary orbit (I (50–60)), and recorded all encounter events together with the velocities and positions of the involved bodies. The two intervals were chosen for two reasons: (i) as mentioned before, it needs time for perturbations induced by the secondary to move inward and distort the innermost ring (in the worst case this needs ~2000 yr, corresponding 30 orbital periods of the binary) and (ii) in all simulations we found periastra alignment for all bodies. For models a1 and a2 this alignment is destroyed very quickly, while for models b1 and b2 in the region of semimajor axes between 1.8 to 3 au this pattern remains for 60 orbits. This periastra alignment might be indicate an enhanced probability of planet formation since even for highly eccentric bodies such an alignment would lead to low differential velocities and consequently to less violent collisions.
We furthermore computed probability distributions where we divided the radial range (0.5–8 au) into 100 rings and counted the number of close encounters within each ring. We distinguished two cases: merging and disruption. A merger occurs if the difference velocity of the two bodies is lower than v_{esc}, otherwise it is assumed as a disruption. This is a very simplified criterion that – compared with a more sophisticated model (Leinhardt & Stewart 2012; Stewart & Leinhardt 2012) – overestimates the disruption case. We did not take into account the more elaborate major collision regimes (hitandrun, disruption, and supercatastrophic disruption) since this is beyond the aim of our study. This means in our case, for instance, that what we call a disruptive event does not necessarily lead to destruction of the bodies, they might also just undergo a hitandrun scenario, which we cannot properly resolve. The application of these regimes will be carried out in a forthcoming article.
For now we wish to show that different types of disc behaviour (dynamically evolving or quiet) can influence protoplanetary orbits such that the probability of further growth may be reduced. In our study the escape velocity of one body is given by km s^{1} (method m1) or for two bodies by km s^{1} (method m2), where M_{p} and R_{p} are the mass and radius of a body. Method m2 is definitely the upper boundary of merging since we assumed that all material has agglomerated into one body. We considered m1 to be a realistic intermediate method as a second indicator to monitor what the common probability for merging or disruption is. Applying method m1, we monitored if one of the bodies could leave the gravitational sphere of influence of the second one, while in m2 we are interested to see whether the body can leave the sphere of influence of the sum of two masses. From this it follows that a disruptive impact occurs when the differential velocity is higher than v_{esc}.
The probability distributions where computed such that the numbers of mergers and disruptive impacts within each bin were divided by the total number of collisions for each time interval I (30–40) or I (50–60). The result is shown in Fig. 6, where plots for model a1 (upper/lower left), model b1 (upper/lower middle), and model b2 (upper/lower right) are shown. The upper panels display the probability distributions for I (30–40) and the lower plots for I (50–60). For model a1 the probability for a merger in I (30–40) is in both cases (m1 and m2) always higher than for disruption and peaks around 3 to 3.5 au. In I (50–60) disruption overwhelms merging in case of m1, but for m2 the distributions for merging and disruption are almost equal and peak around 2.5 to 3 au. For model b1 the situation is different. There the probabilities for disruptions are two to three times higher than for a merger for both methods (m1, m2) and both computational intervals, and they peak around 2 au for disruption and 2.5 au for merger. The distributions for model b2 are almost similar to those for a1. In the first time interval I (30–40) for case of m1 merger and disruption probabilities are similar, where disruption slightly outbalances merger. For m2 plots for merger and disruption are clearly distinct and merger dominates disruption. The distributions peak around 2.5 to 3 au. In the second interval I (50–60) for m1 the distributions for merger and disruption are clearly separated and disruption dominates by a factor of two; for m2 they are similar and the disruption distribution is slightly higher than that of merger. The distributions peak around 2.5 au.
For models a1 and b1 the distributions show different peaks. For a disruption event chances are higher in the inner part of the disc (1–2.5 au), while for a merger the outer part (2–3 au) is more advantageous. Comparing models a1 and b1, one clearly sees that a dynamically evolving gas disc could alter the planet evolution tremendously. This accretionhostile environment can be explained by the action of periodically occuring and inwardmoving spiral waves in the disc that are exited by the secondary together with periodically varying disc eccentricity. This leads to excitations of the embryos and thus higher encounter velocities. But as mentioned before for model b1, we also find that the peak probabilities for merger and disruption are not aligned, giving a small chance for an embryo to grow in the region between 2–3 au. Comparing models a1 and b2, we see almost similar distributions but with a peak shifted towards the star for model b2. For b2 the influence of the disc is visible through slightly elevated disruption distributions, whereas for a merger they did not change. It seems that a cold disc does not hinder planet growth strongly.
For models a1, a2, b1, and b2 we implicitly assumed that the planetesimal accretion phase was successful and only investigated the later stages of planetformation. This condition is far from being granted. Furthermore one has to keep in mind that the given initial distribution of protoplanetary masses and orbital parameters is rather artificial and was chosen so that the gravitational interplay between the disc and the protoplanets could be studied on short timescales, and therefore a detailed study that varies the initial conditions of the particles needs to be conducted.
5. Summary and conclusions
We investigated the evolution of particles embedded in a protostellar disc by using orbital parameters similar to the wellknown γ Cephei system (Neuhäuser et al. 2007). Introducing our hydrodynamical 2D code, which is suitable to compute several thousands of selfgravitating particles that interact gravitationally with a viscous and radiative circumprimary disc within a binary star system, we presented a series of tests that validated our program (see the Appendix A). To our knowledge, this program is the only one at the moment that can handle this many particles and calculates their influence on each other and on the disc. By simulating a coplanar binarydisc system with a grid centred on the primary, we compared the evolution of two distinct models of discprotoplanet interactions with a purely hydrodynamical reference model. Unlike in similar attempts, such as Müller & Kley (2012), the mass loss of the disc during our simulations is small. This is due to the more carefully chosen initial density truncation as well as the subKeplerian initial azimuthal velocity.
As expected from the total particletogasmass ratio of 10^{2}, the disc evolution is not strongly affected by the particles. But a phase drift in the disc’s massweighted eccentricity and argument of pericenter occurs between solutions where the protoplanets’ influence on the disc is either incorporated or neglected.
Furthermore, we showed that Nbody relaxation processes occur faster in presence of an dynamically evolving and interacting gas disc, with average protoplanetary eccentricities almost twice as high as in noninteraction models. We also found that the relaxation time depends on the smoothing parameter, which introduces a numerical viscosity and keeps eccentricities on low values. The growth from embryos to planets within a dynamically evolving gas disc is strongly altered by a dynamically evolving disc, which leads to a decreased probability for planet evolution a least in the inner parts of the gas disc. For a cold disc we found distributions almost similar to a solely Nbody case, but with slightly elevated disruption distributions, which favours them for planet formation. Reasons for these differences between the models are numerous. The influence of the gas and chaotic energy transfer due to close encounters certainly play a role, but more studies are required to identify all contributors. For a gasfree disc we found that our simulation does not show the wavy forced eccentricity patterns encountered in other studies, where the selfgravity of the protoplanets has been neglected (Marzari & Scholl 2000; Thébault et al. 2006; Paardekooper et al. 2008). Such patterns are quickly removed by orbital crossings starting at the outer edge of the system and propagating inside. Since none of the orbitcrossing protoplanets was purged from the simulations, the number of potential highvelocity collisions is most likely overestimated and the vanishing intermittent phasing as well the lack of the wavy forced eccentricity patterns may be nothing but a numerical effect of our applied conditions. There are, of course, several other computational parameters that affect the outcome of the simulations as well, such as the smoothing parameter ϵ regarding the discprotoplanet and protoplanetprotoplanet interaction, the type of boundaries (reflecting, outflow, nonreflecting) of the grid in the hydrodynamical part, different fluxlimiter functions in the advection part of the code, and the orbit evolution of the secondary star. Future studies will investigate their respective influence and shed more light on the problems of planet formation in binary star systems.
Acknowledgments
This work was supported by the FWF Projects P20216N16 and S11608N16. M. Gyergyovits is very grateful to W. Kley and T. Müller for valuable discussions.
References
 Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Barbieri, M., Marzari, F., & Scholl, H. 2002, A&A, 396, 219 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bulirsch, R., & Stoer, J. 1964, Num. Math., 6, 413 [CrossRef] [Google Scholar]
 Chauvin, G., Beust, H., Lagrange, A.M., & Eggenberger, A. 2011, A&A, 528, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Pepe, F., Lovis, C., et al. 2012, Nature, 491, 207 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 [NASA ADS] [Google Scholar]
 Eggleton, P. 2006, Evolutionary Processes in Binary and Multiple Stars (Cambridge: Cambridge University Press) [Google Scholar]
 Fragner, M. M., Nelson, R. P., & Kley, W. 2011, A&A, 528, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Godon, P. 1996, MNRAS, 282, 1107 [NASA ADS] [CrossRef] [Google Scholar]
 Guedes, J. M., Rivera, E. J., Davis, E., et al. 2008, ApJ, 679, 1582 [NASA ADS] [CrossRef] [Google Scholar]
 Haghighipour, N., & Raymond, S. N. 2007, ApJ, 666, 436 [NASA ADS] [CrossRef] [Google Scholar]
 Hatzes, A. P., Cochran, W. D., Endl, M., et al. 2003, ApJ, 599, 1383 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Kaib, N. A., Raymond, S. N., & Duncan, M. 2013, Nature, 493, 381 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
 Kley, W. 1999, MNRAS, 303, 696 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., & Crida, A. 2008, A&A, 487, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., & Nelson, R. P. 2008, A&A, 486, 617 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., Papaloizou, J. C. B., & Ogilvie, G. I. 2008, A&A, 487, 671 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Korchagin, V., & Theis, C. 1999, A&A, 347, 442 [NASA ADS] [Google Scholar]
 Lagrange, A.M., Beust, H., Udry, S., Chauvin, G., & Mayor, M. 2006, A&A, 459, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1959, Fluid Mech. (Pergamon press, Reading Mass, Addison – Weley Pub. Co.) [Google Scholar]
 Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, 745, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Quintana, E. V., Chambers, J. E., Duncan, M. J., & Adams, F. C. 2004, in Rev. Mex. Astron. Astrofis. Conf. Ser. 22, eds. G. GarciaSegura, G. TenorioTagle, J. Franco, & H. W. Yorke, 99 [Google Scholar]
 Marzari, F., & Scholl, H. 2000, ApJ, 543, 328 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Müller, T. W. A., & Kley, W. 2012, A&A, 539, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neuhäuser, R., Mugrauer, M., Fukagawa, M., Torres, G., & Schmidt, T. 2007, A&A, 462, 777 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paardekooper, S.J., Thébault, P., & Mellema, G. 2008, MNRAS, 386, 973 [NASA ADS] [CrossRef] [Google Scholar]
 PilatLohinger, E., & Dvorak, R. 2002, Celest. Mech. Dyn. Astron., 82, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Quintana, E. V., Adams, F. C., Lissauer, J. J., & Chambers, J. E. 2007, ApJ, 660, 807 [NASA ADS] [CrossRef] [Google Scholar]
 Santos, N. C., Mayor, M., Naef, D., et al. 2000, A&A, 361, 265 [NASA ADS] [Google Scholar]
 Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Sod, G. A. 1978, J. Comput. Phys., 27, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Stewart, S. T., & Leinhardt, Z. M. 2012, ApJ, 751, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Thebault, P. 2011, Celest. Mech. Dyn. Astron., 111, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Thébault, P., Marzari, F., Scholl, H., Turrini, D., & Barbieri, M. 2004, A&A, 427, 1097 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thébault, P., Marzari, F., & Scholl, H. 2006, Icarus, 183, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Theis, C., & Orlova, N. 2004, A&A, 418, 959 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Xie, J.W., & Zhou, J.L. 2009, ApJ, 698, 2066 [NASA ADS] [CrossRef] [Google Scholar]
 Xie, J.W., Zhou, J.L., & Ge, J. 2010, ApJ, 708, 1566 [NASA ADS] [CrossRef] [Google Scholar]
 Zucker, S., Mazeh, T., Santos, N. C., Udry, S., & Mayor, M. 2004, A&A, 426, 695 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
Appendix A: standard test for the code
Fig. A.1 Normalised surface density for τ = 0.016, τ = 0.021, τ = 0.058, τ = 0.091. The solid lines indicate the analytical solution. The initial profile was set according to the analytical solution at τ = 0.016. For details see Sect. A.1. 

Open with DEXTER 
To validate our code we ran several test cases proposed in literature. The following subsections briefly introduce the tests performed to monitor the respective influence of the viscosity tensor (Sect. A.1), a gravitational perturber (Sect. A.2) as well as energy loss via irradiation (Sect. A.3) onto the evolution of the circumprimary disc.
Appendix A.1: Viscosity tensor
Because our code applies a viscosity tensor, we used the method reported by Kley (1999) to test it. We computed the axisymmetric problem of an expanding ring whose initial surface density profile corresponds to a δdistribution spreading under the influence of a low viscosity and negligible pressure. Assuming Keplerian rotation, the analytical solution of the surface density Σ(r) can be given in terms of Bessel functions: (A.1)Here, x denotes the normalised radius x = r/r_{0} where the initial ring was located at r_{0}, and I_{1 / 4} is the modified Bessel function. The dimensionless time is denoted by τ = t/t_{v} given in units of the viscous spreading time . We computed the full equations of motion on a 128 × 128 grid using an isothermal approach where the sound speed is set to zero. Figure A.1 shows the evolution of the ring for different dimensionless times. Initially, the ring is located at r_{0} = 1. The starting time τ = 0.016 was chosen to be consistent with the initial conditions of Kley (1999), so that the initial distribution no longer resembles a perfect delta distribution. Instead, the ring has already spread somewhat. Our constant dimensionless viscosity was chosen to be ν = 10^{6} and no artificial viscosity was introduced. The numerical solution (points) agrees well with the analytical results (lines).
Appendix A.2: Disc–planet interaction
The disc–planet interaction was tested via recomputing the case of Jupiter in a viscous and a nonviscous disc, as can be found in de ValBorro et al. (2006). We also applied an additional artificial viscosity of von Neumann & Richtmyer type (Stone & Norman 1992), where C^{2} = 2. According to the test conditions given in this article, we used n_{r} × n_{φ} = (128,384) cells in radial and azimuthal direction, reflecting boundaries, and applied a damping function (A.2)where x represents the surface density and velocity components, τ is the orbital period at the corresponding boundary, and R(r) is a parabolic function (de ValBorro et al. 2006, Eq. (10)). The semimajor axis of the planet was set to a = 1 au and the planettostarmass ratio μ was 0.001. The initial surface density is constant: (A.3)This entails a total disc mass of M_{d} = 0.012 M_{⊙}. M_{s} is the mass of the star. We used a locally isothermal approach and H/r = 0.05, which results in c_{s} = 0.05v_{K} where c_{s} is the isothermal sound speed and v_{K} is the local Keplerian velocity. The grid extends from 0.4 to 2.5 au, and we applied the artificial damping of Eq. (A.2) between 0.4 and 0.5 as well as between 2.1 and 2.5 au. The planet mass was slowly increased to its final value during its first five orbits according to (A.4)where P_{p} is the period of the planet. The kinematic viscosity was chosen to be ν = 1 × 10^{5} in dimensionless units, where a = 1, P_{p} = 2π and M_{s} = 1 − μ. We can reproduce the normalised azimuthally averaged surface density after 100 planetary orbits for the nonviscous case and viscous case and also the evolution of the total torque for the whole computation time of 200 orbits shown in de ValBorro et al. (2006). Averaging the torques between 175–200 periods in dimensionless units produces values of − 7.73 × 10^{5} for the viscous case and − 1.75 × 10^{5} for the inviscid case. The time is given in orbital periods of the planet. Compared with the values of Tables 5 and 6 in de ValBorro et al. (2006), our results are similar to the outcome of RH2D, a code written by Kley (1989).
Appendix A.3: Radiative disc
Fig. A.2 Comparison of analytical and simulated disc temperatures for cases I and II. The radius is normalised r = a/ (5.2 au) and the temperature is given in Kelvin. 

Open with DEXTER 
To test our implementation of the energy equation we compared our code with the fully radiative reference model by Kley & Crida (2008) (see also Eq. (6) of this article). The protoplanetary disc is treated as a twodimensional, nonselfgravitating gas that can be described by the NavierStokes equations. The embedded planet is modelled as a point mass that orbits the central star on a fixed, circular orbit. For the planetary potential, we used a smoothing length of = 0.6H, where H is the vertical scale height of the disc. To calculate the gravitational torques acting on the planet, we excluded the inner parts of the Hill sphere of the planet. The test was carried out twice.
Case I: in this case we implemented the energyequation given in Kley & Crida (2008) and alsoused the adiabatic sound speed.
Case II: for the second case the method suggested in Müller & Kley (2012) was implemented using only the isothermal sound speed.
The numerical parameters for both cases were chosen as follows: a twodimensional (r − φ) computational domain consisting of a ringshaped protoplanetary disc, which is centred on the star. This domain extends from r_{min} = 0.4 a_{Jup} to r_{max} = 2.5a_{Jup}, where a_{Jup} = 5.2 au is the semimajor axis of Jupiter. The total disc mass was set to M_{d} = 0.01 M_{⊙} and the constant kinematic viscosity ν = 10^{15} cm^{2}/ s. Following a power law, the initial density distribution of the gas is given by Σ(r) = Σ_{0}r^{− 1 / 2}, the temperature T(r) ∝ r^{1} and the rotation is subKeplerian. A simulation prerun of 200 orbits of Jupiter was performed to ensure that the disc is in equilibrium for the actual test calculations. The disc boundary conditions were chosen to be reflecting, and dampening was applied in the regions (0.4a_{Jup}, 0.5a_{Jup}) and (2.1a_{Jup}, 2.5a_{Jup}).The analytical temperature profile ranging between 0 K and 170 K was derived from the equilibrium of heating and cooling (Q_{+} = Q_{−}) as well as from fact that the Rosseland mean opacity is κ = κ_{0}T^{2}. We furthermore assumed that the radial velocity of the disc is zero. This gives an equation of sixth order in T, which can be reduced to third order and finally solved to arrive at one real solution. In Fig. A.2 we show a comparison of simulated and analytical disc temperatures for cases I and II. The analytical and computed solutions match very well except for high temperatures, where the assumption κ = κ_{0}T^{2} is no longer valid. Comparing cases I and II, the temperature profile in the latter case is lower and the point where analytical and numerical solutions branch moves inward. At this point κ starts to follow a different power law in T. A similar behaviour can be seen in the aspect ratio (Fig. A.3). In both cases the aspect ratio does not increase steadily until the inner boundary is reached. It reaches a highest value of 0.059 at 0.77 a_{Jup} for case I and 0.0526 at 0.6 a_{Jup} for case II. The values in case I are higher than in case II and the maximum of the numerical solution is shifted to smaller radii for case II compared with case I. This is expected since and H/r = c_{s}/v_{K}.
Fig. A.3 Same as Fig. A.2, but for aspect ratios H/r. 

Open with DEXTER 
Fig. A.4 Time evolution of the specific total torque T ((a^{2}Ω^{2})_{Jup}) exerted by the disc on a planet of 20 M_{Earth} for the adiabatic, local cool/heat, and fully radiative case. 

Open with DEXTER 
In a second test we recomputed the specific torque the disc exerts on a 20 M_{Earth} planet using the implementations of case I for three modes, neglecting different terms in the energy equation (Kley & Crida 2008, Fig. 2). Omitting all terms except the first one on the right side in Eq. (6) gives an adiabatic model, whereas in the cool/heat case we just neglected the last term. In the fully radiative case we included all terms in the equation. In Fig. A.4 we show that we can find an evolution of the specific torque the disc exerts on the planet, which is similar to the outcome of Kley & Crida (2008, Fig. 2), especially during the time between 200 and 300 planetary orbits. The local cool/heat case matches very well, whereas the specific torque in the fully
radiative case is slightly higher than in Kley & Crida (2008). The adiabatic case shows stronger oscillations in the beginning and approaches a higher equilibrium value than presented in Kley & Crida (2008). These small differences might arise because we used subKeplerian azimuthal velocities for the disc and damped towards these values instead of purely Keplerian ones.
All Tables
All Figures
Fig. 1 Time evolution of the massweighted eccentricity for a disc without particles (reference model) and model a2 (upper plot), and models a1 and b1 (lower plot). Time is given in terms of binary orbits, where one orbit corresponds to 66.7 yr. See Sect. 3.1 for details. 

Open with DEXTER  
In the text 
Fig. 2 Time evolution of the massweighted argument of the pericenter for a disc without particles (reference model) and model a2 (upper plot), and models a1 and b1 (lower plot). Time is given in terms of binary orbits, where one orbit corresponds to 66.7 yr. 

Open with DEXTER  
In the text 
Fig. 3 Time evolution of massweighted argument of the pericenter for a disc without particles (upper plot) and of massweighted eccentricity (lower plot) for 50 000 yr after applying a running window average. Time is given in terms of binary orbits, where each orbit corresponds to 66.7 yr. We find transition from circulation to oscillation within 200 binary orbits for ϖ_{mw} and a subsequent oscillation around ≈0.6 rad. A damped oscillation is visible for e_{mw}, reaching a value around ≈0.0275. 

Open with DEXTER  
In the text 
Fig. 4 Time evolution of the particle semimajor axis influenced by the disc (model b1). We find a transition from ordered motion close to the initial positions on the grid (left lower corner) to a distribution of semimajor axes that spreads across the whole stable region of the disc within 900 yr. 

Open with DEXTER  
In the text 
Fig. 5 Evolution of protoplanet mean (upper plot) and root mean square eccentricity (lower plot) for all four models. Highest values in both plots are reached by model b, whereas model b1 shows a similar behaviour as model a, and the lowest values are reached by model a1. 

Open with DEXTER  
In the text 
Fig. 6 Collision probabilities for 30–40 orbital periods for model a1 (upper left), model b1 (upper middle), and model b2 (upper right) as well as for 50–60 orbits for model a1 (lower left), model b1 (lower right), and model b2 (lower right). We show the plot radius (in au) versus the probability for an event (in percent) for disruption and merging for methods m1 and m2 (explained in the text). 

Open with DEXTER  
In the text 
Fig. A.1 Normalised surface density for τ = 0.016, τ = 0.021, τ = 0.058, τ = 0.091. The solid lines indicate the analytical solution. The initial profile was set according to the analytical solution at τ = 0.016. For details see Sect. A.1. 

Open with DEXTER  
In the text 
Fig. A.2 Comparison of analytical and simulated disc temperatures for cases I and II. The radius is normalised r = a/ (5.2 au) and the temperature is given in Kelvin. 

Open with DEXTER  
In the text 
Fig. A.3 Same as Fig. A.2, but for aspect ratios H/r. 

Open with DEXTER  
In the text 
Fig. A.4 Time evolution of the specific total torque T ((a^{2}Ω^{2})_{Jup}) exerted by the disc on a planet of 20 M_{Earth} for the adiabatic, local cool/heat, and fully radiative case. 

Open with DEXTER  
In the text 