© ESO, 2014

1. Introduction

Most main-sequence 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 mid-separation binaries, with semi-major 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 mid-separation 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 planet-forming 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 gas-giants can form. During accretion phases to km-sized 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 km-sized 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 size-dependent phasing of protoplanetary orbits that decreases dv for equal-sized 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 N-body 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 gas-disc 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 self-gravity 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 self-gravitation of the protoplanets and the gravitational interaction with the disc is more important than gas drag.

We here introduce a GPU-CPU 2D hydrodynamic grid code that combines hydrodynamic radiative disc computations with highly accurate N-body simulations. The development of this new code allows us to investigate the interplay between a circumprimary gas disc and the orbital evolution of embedded self-gravitating protoplanets on time-scales 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 disc-particle 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 ZEUS-2D solver for hydrodynamical equations on a two-dimensional 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 staggered-mesh finite-difference method using operator splitting and a first-order integrator in time to handle the source step. The advective terms are solved by a second-order 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 so-called FARGO acceleration algorithm (Masset 2000) to speed up calculations; (iii) we partly rewrote the code to run on NVIDIA graphics cards in double-precision 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

$\begin{array}{ccc}& & {\mathit{\partial }}_{\mathit{t}}{\mathit{v}}_{\mathit{r}}\mathrm{+}{\mathit{v}}_{\mathit{r}}{\mathit{\partial }}_{\mathit{r}}{\mathit{v}}_{\mathit{r}}\mathrm{+}\frac{{\mathit{v}}_{\mathit{\phi }}}{\mathit{r}}{\mathit{\partial }}_{\mathit{\phi }}{\mathit{v}}_{\mathit{r}}\mathrm{-}\frac{{\mathit{v}}_{\mathit{\phi }}^{\mathrm{2}}}{\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{\Sigma }}{\mathit{\partial }}_{\mathit{r}}\mathit{P}\mathrm{+}\mathrm{▽}\mathrm{·}{\mathrm{T}}_{\mathit{r}}\mathrm{-}{\mathit{\partial }}_{\mathit{r}}\mathrm{\Psi }\mathrm{+}{\mathit{f}}_{\mathit{r}}\\ & & {\mathit{\partial }}_{\mathit{t}}{\mathit{v}}_{\mathit{\phi }}\mathrm{+}{\mathit{v}}_{\mathit{r}}{\mathit{\partial }}_{\mathit{r}}{\mathit{v}}_{\mathit{\phi }}\mathrm{+}\frac{{\mathit{v}}_{\mathit{\phi }}}{\mathit{r}}{\mathit{\partial }}_{\mathit{\phi }}{\mathit{v}}_{\mathit{\phi }}\mathrm{+}\frac{{\mathit{v}}_{\mathit{\phi }}{\mathit{v}}_{\mathit{r}}}{\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{\Sigma }\mathit{r}}{\mathit{\partial }}_{\mathit{\phi }}\mathit{P}\\ & & \end{array}$(1)

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 non-inertial 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 $\begin{array}{ccc}\mathrm{▽}\mathrm{·}{\mathrm{T}}_{\mathit{r}}\mathrm{=}\frac{\mathrm{1}}{\mathit{r}}\mathrm{[}{\mathit{\partial }}_{\mathit{r}}\mathrm{(}\mathit{r}{\mathit{T}}_{\mathit{rr}}\mathrm{)}\mathrm{+}{\mathit{\partial }}_{\mathit{\phi }}{\mathit{T}}_{\mathit{r\phi }}\mathrm{]}\mathrm{-}\frac{{\mathit{T}}_{\mathit{\phi \phi }}}{\mathit{r}}& & \end{array}$(2)and $\begin{array}{ccc}\mathrm{▽}\mathrm{·}{\mathrm{T}}_{\mathit{\phi }}\mathrm{=}{\mathit{\partial }}_{\mathit{r}}\mathrm{(}{\mathit{T}}_{\mathit{r\phi }}\mathrm{)}\mathrm{+}{\mathit{\partial }}_{\mathit{\phi }}\left(\frac{{\mathit{T}}_{\mathit{\phi \phi }}}{\mathit{r}}\right)\mathrm{+}\mathrm{2}\frac{{\mathit{T}}_{\mathit{r\phi }}}{\mathit{r}}\mathit{,}& & \end{array}$(3)where T_rr, T_φφ, T_rφ are the components of the viscous stress tensor: $\begin{array}{ccc}\begin{array}{c}\\ \\ \\ \\ \\ \end{array}& & \end{array}$(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 post-shock oscillations in the Sod shock-tube 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 $\begin{array}{ccc}\mathrm{\Psi }\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{1}}}{\mathrm{|}{r}\mathrm{-}{{r}}_{\mathrm{1}}\mathrm{|}}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{2}}}{\mathrm{|}{r}\mathrm{-}{{r}}_{\mathrm{2}}\mathrm{|}}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{3}}^{\mathit{n}}\frac{\mathit{G}{\mathit{M}}_{\mathit{i}}}{\mathrm{|}{r}\mathrm{-}{{r}}_{\mathit{i}}\mathrm{+}\mathit{ϵ}\mathrm{|}}\mathit{,}& & \end{array}$(5)where r₁, r₂ are the position vectors and M₁, M₂ 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 Val-Borro 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 $\begin{array}{ccc}{\mathit{\partial }}_{\mathit{t}}\mathit{e}\mathrm{+}\mathrm{▽}\mathrm{·}\mathrm{(}\mathit{e}{{v}}^{\mathrm{)}}\mathrm{=}\mathrm{-}\mathit{p}\mathrm{▽}\mathrm{·}{v}\mathrm{+}{\mathit{Q}}_{\mathrm{+}}\mathrm{-}{\mathit{Q}}_{\mathrm{-}}\mathrm{-}\mathrm{2}\mathit{H}\mathrm{▽}\mathrm{·}{F}\mathit{,}& & \end{array}$(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 flux-limited diffusion approximation. The flux is given by $\begin{array}{ccc}{F}\mathrm{=}\mathrm{-}\frac{\mathit{\lambda c}\mathrm{4}\mathit{a}{\mathit{T}}^{\mathrm{3}}}{\mathit{\rho \kappa }}\mathrm{▽}\mathit{T,}& & \end{array}$(7)where c denotes the speed of light, a the radiation constant, ρ the mid-plane density, κ the Rosseland mean opacity, and λ is the flux-limiter Kley (1989, Eq. (9)) and Kley & Crida (2008). The initial radial velocity v_r was set to zero, and v_φ is sub-Keplerian, 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 N-body 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 semi-major 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 mass-weighted 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 binary-discinteraction.

model a1: several thousand self-gravitating protoplanets distort the disc gravitationally but no back-reaction 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 two-dimensional 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 semi-major 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 self-gravitating 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 dust-to-gas-ratio 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₀Δφ, where r₀ = 1 au. Thus the last ring is close to the critical semi-major axis a = 3.920 ± 0.009 au where stable motion is still possible according to the studies by Holman & Wiegert (1999) and Pilat-Lohinger & 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.

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 time-steps. 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 self-gravity 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 zero-gradient boundary condition will create a (eccentric) hole in the vicinity of this boundary (Kley & Nelson 2008). Taking a non-reflecting 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 power-law distribution: $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}{\mathit{r}}^{\mathrm{-}\mathit{n}}\mathit{f}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,} \mathit{n}\mathrm{=}\mathrm{1}\mathit{,}& & \end{array}$(8)where Σ₀ is the value of the vertical averaged density at r = 1. To avoid unrealistically low gas-densities we applied a floor value (fl), which is given by fl = Σ₀ × 10^-7. The density was modified by the function f(r), which emulates a truncated disc. For this purpose we introduced the Fermi function $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{e}}^{\mathit{\delta }\mathrm{(}\mathit{r}\mathrm{-}{\mathit{r}}_{\mathit{t}}\mathrm{)}}}\mathit{,} \mathit{\delta }\mathrm{=}\frac{\ln \mathrm{(}{\mathrm{10}}^{\mathrm{7}}\mathit{/}\mathit{r}\mathrm{-}\mathrm{1}\mathrm{)}}{{\mathit{r}}_{\max }\mathrm{-}{\mathit{r}}_{\mathit{t}}}& & \end{array}$(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 mass-weighted eccentricity (e_mw) and argument of periastron of the gas disc (ϖ_mw) where we followed Müller & Kley (2012), and which are given by $\begin{array}{ccc}{\mathit{e}}_{\mathrm{mw}}& \mathrm{=}& \left[{\mathrm{\int }}_{{\mathit{r}}_{\min }}^{{\mathit{r}}_{\max }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{\Sigma }\mathit{er}\mathrm{d}\mathit{r}\mathrm{d}\mathit{\phi }\right]\mathrm{\times }{\left[{\mathrm{\int }}_{{\mathit{r}}_{\min }}^{{\mathit{r}}_{\max }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{\Sigma }\mathit{r}\mathrm{d}\mathit{r}\mathrm{d}\mathit{\phi }\right]}^{-1}\\ {\mathit{\varpi }}_{\mathrm{mw}}& \mathrm{=}& \left[{\mathrm{\int }}_{{\mathit{r}}_{\min }}^{{\mathit{r}}_{\max }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{\Sigma }\mathit{\varpi r}\mathrm{d}\mathit{r}\mathrm{d}\mathit{\phi }\right]\mathrm{\times }{\left[{\mathrm{\int }}_{{\mathit{r}}_{\min }}^{{\mathit{r}}_{\max }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{\Sigma }\mathit{r}\mathrm{d}\mathit{r}\mathrm{d}\mathit{\phi }\right]}^{-1}\end{array}$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

In Figs. 1–3 we show the evolution of the mass-weighted eccentricity and mass-weighted 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 N-body 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

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 particle-particle 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 gas-free disc we find that our simulation does not show the wavy forced eccentricity patterns encountered in other studies where the self-gravity 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.

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 semi-major 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 (hit-and-run, disruption, and super-catastrophic 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 hit-and-run 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 ${\mathit{v}}_{\mathrm{esc}}\mathrm{=}\sqrt{\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{p}}\mathit{/}{\mathit{R}}_{\mathrm{p}}}s\mathrm{1.957}$ km s^-1 (method m1) or for two bodies by ${\mathit{v}}_{\mathrm{esc}}\mathrm{=}\sqrt{\mathrm{4}\mathit{G}{\mathit{M}}_{\mathrm{p}}\mathit{/}{\mathit{R}}_{\mathrm{p}}}s\mathrm{2.768}$ 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 accretion-hostile environment can be explained by the action of periodically occuring and inward-moving 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 planet-formation. 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 time-scales, 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 well-known γ Cephei system (Neuhäuser et al. 2007). Introducing our hydrodynamical 2D code, which is suitable to compute several thousands of self-gravitating 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 binary-disc system with a grid centred on the primary, we compared the evolution of two distinct models of disc-protoplanet 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 sub-Keplerian initial azimuthal velocity.

As expected from the total particle-to-gas-mass ratio of 10^-2, the disc evolution is not strongly affected by the particles. But a phase drift in the disc’s mass-weighted eccentricity and argument of pericenter occurs between solutions where the protoplanets’ influence on the disc is either incorporated or neglected.

Furthermore, we showed that N-body relaxation processes occur faster in presence of an dynamically evolving and interacting gas disc, with average protoplanetary eccentricities almost twice as high as in non-interaction 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 N-body 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 gas-free disc we found that our simulation does not show the wavy forced eccentricity patterns encountered in other studies, where the self-gravity 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 orbit-crossing protoplanets was purged from the simulations, the number of potential 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 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 disc-protoplanet and protoplanet-protoplanet interaction, the type of boundaries (reflecting, outflow, non-reflecting) of the grid in the hydrodynamical part, different flux-limiter 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 P-20216-N16 and S11608-N16. M. Gyergyovits is very grateful to W. Kley and T. Müller for valuable discussions.

References

Appendix A: standard test for the code

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: $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{(}{\mathit{r}}^{\mathrm{)}}\mathrm{=}\frac{\mathit{C}}{\mathit{\tau }{\mathit{x}}^{\mathrm{1}\mathit{/}\mathrm{4}}}\exp \left(\mathrm{-}\frac{\mathrm{1}\mathrm{+}{\mathit{x}}^{\mathrm{2}}}{\mathit{\tau }}\right){\mathit{I}}_{\mathrm{1}\mathit{/}\mathrm{4}}\left(\frac{\mathrm{2}\mathit{x}}{\mathit{\tau }}\right)\mathrm{·}& & \end{array}$(A.1)Here, x denotes the normalised radius x = r/r₀ where the initial ring was located at r₀, $\mathit{C}\mathrm{=}\mathrm{1}\mathit{/}\mathrm{\left(}\mathit{\pi }{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\right)}$ 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 ${\mathit{t}}_{\mathit{v}}\mathrm{=}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{12}\mathit{\nu }\mathrm{\right)}$. 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₀ = 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 non-viscous disc, as can be found in de Val-Borro et al. (2006). We also applied an additional artificial viscosity of von Neumann & Richtmyer type (Stone & Norman 1992), where C² = 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 $\begin{array}{ccc}\frac{\mathrm{d}\mathit{x}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{0}}}{\mathit{\tau }}\mathit{R}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}& & \end{array}$(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 Val-Borro et al. 2006, Eq. (10)). The semi-major axis of the planet was set to a = 1 au and the planet-to-star-mass ratio μ was 0.001. The initial surface density is constant: $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{=}\mathrm{0.002}\frac{{\mathit{M}}_{\mathrm{s}}}{\mathit{\pi }{\mathit{a}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(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 $\begin{array}{ccc}\mathit{M}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{M}}_{\mathrm{p}}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\pi t}}{\mathrm{10}{\mathit{P}}_{\mathrm{p}}}\right)\mathit{,} \mathit{t}\mathrm{=}\mathrm{0}\mathrm{to}\mathrm{5}{\mathit{P}}_{\mathrm{p}}\mathit{,}& & \end{array}$(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 non-viscous case and viscous case and also the evolution of the total torque for the whole computation time of 200 orbits shown in de Val-Borro 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 Val-Borro et al. (2006), our results are similar to the outcome of RH2D, a code written by Kley (1989).

Appendix A.3: Radiative disc

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 two-dimensional, non-self-gravitating gas that can be described by the Navier-Stokes 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 two-dimensional (r − φ) computational domain consisting of a ring-shaped 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 semi-major axis of Jupiter. The total disc mass was set to M_d = 0.01 M_⊙ and the constant kinematic viscosity ν = 10¹⁵ cm²/ s. Following a power law, the initial density distribution of the gas is given by Σ(r) = Σ₀r^{− 1 / 2}, the temperature T(r) ∝ r^-1 and the rotation is sub-Keplerian. A simulation pre-run 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 κ = κ₀T². 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 κ = κ₀T² 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 ${\mathit{c}}_{\mathrm{s}}\mathrm{\propto }\sqrt{\mathit{T}}$ and H/r = c_s/v_K.

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 sub-Keplerian azimuthal velocities for the disc and damped towards these values instead of purely Keplerian ones.

All Tables

All Figures