© ESO, 2017

1. Introduction

Circumbinary discs are accretion discs that orbit a binary system that consists for example of a binary star or a binary black hole system. A very prominent example of a circumbinary disc orbiting two stars is the system GG Tau. Due to the large star separation, the whole system (stars and disc) can be directly imaged by interferometric methods (Guilloteau et al. 1999). Later data yielded constraints on the size of the dust in the system and have indicated that the central binary may consist of multiple stellar systems (see Andrews et al. 2014; Di Folco et al. 2014, and references therein). More recent observations have pointed to possible planet formation (Dutrey et al. 2014) and streamers from the disc onto the stars (Yang et al. 2017). A summary of the properties of the GG Tau system is given in Dutrey et al. (2016).

An additional important clue for the existence and importance of circumbinary discs is given by the observed circumbinary planets. These are planetary systems where the planets do not orbit one star, but instead orbit a binary star system. These systems have been detected in recent years by the Kepler Space Mission and have inspired an intense research activity. The first such circumbinary planet, Kepler-16b, orbits its host system consisting of a K-type main-sequence star and an M-type red dwarf, in what is known as a P-type orbit with a semi-major axis of 0.7048 au and a period of 228.7 days (Doyle et al. 2011). Presently there are ten Kepler circumbinary planets known and a summary of the properties of the first five systems discovered is given in Welsh et al. (2014)¹. The observations have shown that in these systems the stars are mutually eclipsing each other, as well as the planets. The simultaneous eclipses of several objects is a clear indication of the flatness of the system because the orbital plane of the binary star has to coincide nearly exactly with the orbital plane of the planet. In this respect these systems are flatter than our own solar system. Because planets form in discs, the original protoplanetary disc in these systems had to orbit around both stars, hence it must have been a circumbinary disc that was also coplanar with respect to the orbital plane of the central binary.

In addition to their occurrence around younger stars, circumbinary discs are believed to orbit supermassive black hole binaries in the centres of galaxies (Begelman et al. 1980), where they may play an important role in the evolution of the host galaxies (Armitage & Natarajan 2002). In the early phase of the Universe there were many more close encounters between the young galaxies that already had black holes in their centres. The gravitational tidal forces between them often resulted in new merged objects that consisted of a central black hole binary surrounded by a large circumbinary disc (see e.g. Cuadra et al. 2009; Lodato et al. 2009, and references therein). The dynamical behaviour of such a system is very similar to the one described above, where we have two stars instead of black holes.

To understand the dynamics of discs around young binary stars such as GG Tau or the planet formation process in circumbinary discs or the dynamical evolution of a central black hole binary in a galaxy, it is important to understand in detail the behaviour of a disc orbiting a central binary. The strong gravitational perturbation by the binary generates spiral waves in the disc; these waves transport energy and angular momentum from the binary to the disc (Pringle 1991) and subseqently alter the evolution of the binary (Artymowicz et al. 1991). The most important impact of the angular momentum transfer to the disc is the formation of an inner cavity in the disc with a very low density whose size depends on disc parameters (such as viscosity or scale height) and on the binary properties (mass ratio, eccentricity) (Artymowicz & Lubow 1994).

The numerical simulation of these circumbinary discs is not trivial and over the years several attempts have been made. The first simulations of circumbinary discs used the smoothed particle hydrodynamics (SPH) method, a Lagrangian method where the fluid is modelled as individual particles (Artymowicz & Lubow 1996). On the one hand, this method has the advantage that the whole system can be included in the computational domain and hence the accretion onto the single binary components can be studied. On the other hand, because of the finite mass of the SPH-particles, the maximum density contrast that can be resolved is limited.

Later Rozyczka & Laughlin (1997) simulated circumbinary discs with the help of a finite difference method on a polar grid. The grid allowed for a much larger range in mass resolution than the SPH codes used to date. A polar grid is well suited for the outer parts of a disc, but suffers from the fact that there will always be an inner hole in the computational domain since the minimum radius R_min cannot be made arbitrarily small. Decreasing the minimum radius will also strongly reduce the required numerical time step, rendering long simulations unfeasible. But since circumbinary disc simulations have to be simulated for tens of thousands of binary orbits to reach a quasi-steady state, the minimum radius is often chosen such that the motion of the central stars is not covered by the computational domain, and the mass transfer from the disc to the binary cannot be studied.

Günther & Kley (2002) used a special dual-grid technique to use the best of the two worlds, the high resolution of the grid codes and the coverage of the whole domain of the particle codes. For the outer disc they used a polar grid and they overlaid the central hole with a Cartesian grid. In this way it was possible to study the complex interaction between the binary and the disc in greater detail. All these simulations indicated that despite the formation of an inner gap, material can enter the inner regions around the stars crossing the gap through stream like features to eventually become accreted onto the stars and influence their evolution (Bate et al. 2002).

Another interesting feature of circumbinary discs concerns their eccentricity. As shown in hydrodynamical simulations for planetary mass companions, discs develop a global eccentric mode even for circular binaries (Papaloizou et al. 2001). The eccentricity is confined to the inner disc region and shows a very slow precession (Kley & Dirksen 2006). These results were confirmed for equal mass binaries on circular orbits (MacFadyen & Milosavljević 2008). The disc eccentricity is excited by the 3:1 outer eccentric Lindblad resonance (Lubow 1991; Papaloizou et al. 2001), which is operative for a sufficiently cleared out gap. For circular binaries a transition, in the disc structure from more circular to eccentric is thought to occur for a mass ratio of secondary to primary above 1/25 (D’Orazio et al. 2016). However, this appears to contradict the mentioned results by Papaloizou et al. (2001) and Kley & Dirksen (2006) who show that this transition already occurs for a secondary in the planet mass regime with details depending on disc parameters such as pressure and viscosity.

Driven by the observation of planets orbiting eccentric binary stars and the possibility of studying eccentric binary black holes, a large number of numerical simulations have been performed over the last few years dealing with discs around eccentric central binaries (e.g. Pierens & Nelson 2013; Kley & Haghighipour 2014; Dunhill et al. 2015). The recent simulations of circumbinary discs by Pierens & Nelson (2013), Kley & Haghighipour (2014) and Lines et al. (2015) used grid-based numerical methods on a polar grid with an inner hole. This raises questions about the location and imposed boundary conditions at the innermost grid radius R_min. In particular, the value of R_min must be small enough to capture the development of the disc eccentricity through gravitational interaction between the binary and the disc, especially through the 3:1 Lindblad resonance (Pierens & Nelson 2013). Therefore R_min has to be chosen in a way that all important resonances lie inside the computational domain.

In addition to the position of the inner boundary, there is no common agreement on which numerical boundary condition better describes the disc. In simulations concerning circumbinary planets mainly two types of boundary conditions are used: closed inner boundaries (Pierens & Nelson 2013; Lines et al. 2015) that do not allow for mass flow onto the central binary, and outflow boundaries (Kley & Haghighipour 2014, 2015) that do allow for accretion onto the binary. While Lines et al. (2015) do not find a significant difference between these two cases, Kley & Haghighipour (2014) see a clear impact on the surface density profile, and they were also able to construct discs with (on average) constant mass flow through the disc, which is not possible for closed boundaries. Driven by these discussions in the literature we decided to perform dedicated numerical studies to test the impact of the location of R_min, the chosen boundary condition, and other numerical aspects in more detail. During our work we became aware of other simulations that were tackling similar problems. In October 2016 alone three publications appeared on astro-ph that described numerical simulations of circumbinary discs (Fleming & Quinn 2017; Miranda et al. 2017; Mutter et al. 2017). While the first paper considered SPH-simulations, in the latter two papers grid-codes were used and in both the boundary condition at R_min is discussed.

In our paper we start out with numerical considerations and investigate the necessary conditions at the inner boundary in detail, and we discuss aspects of the other recent results in the presentation of our findings below. Additionally, we present a detailed parameter study of the dynamical behaviour of circumbinary discs as a function of binary and disc properties.

We organised this paper in the following way. In Sect. 2 we describe the numerical and physical setup of our circumbinary simulations. In Sect. 3 we briefly describe the numerical methods of the different codes we used for our simulations. In Sect. 4 we examine the disc structure by varying the inner boundary condition and its location. The influence of different binary parameters are examined in Sect. 5. Different disc parameters and their influence are studied in Sect. 6. Finally we summarise and discuss our results in Sect. 7.

2. Model setup

To study the evolution of circumbinary discs we perform locally isothermal hydrodynamical simulations. As a reference system we consider a binary star having the properties of Kepler-16, whose dynamical parameters are presented in Table 1. This system has a typical mass ratio and a moderate eccentricity of the orbit, and it has been studied frequently in the literature. In our parameter study we start from this reference system and vary the different aspects systematically.

2.1. Physics and equations

Inspired by the flatness of the observed circumbinary planetary systems, for example in the Kepler-16 system the motion takes place in single plane to within 0.5° (Doyle et al. 2011), we make the following two assumptions:

• 1.

  the vertical thickness H of the disc is small compared to the distance from the centre R;

• 2.

  there is no vertical motion.

It is therefore acceptable to reduce the problem to two dimensions by vertically averaging the hydrodynamical equations. In this case it is meaningful to work in cylindrical coordinates (R,ϕ,z)² and the averaging process is done in the z-direction. Furthermore, we choose the centre of mass of the binary as the origin of our coordinate system with the vertical axis aligned with the rotation axis of the binary. The averaged hydrodynamical equations are then given by where is the surface density, the vertically integrated pressure, and u = (u_R,u_ϕ)^T the two-dimensional velocity vector. We close this system of equations with a locally isothermal equation of state (3)with the local sound speed c_s(R).

The gravitational potential Φ of both stars is given by (4)Here M_k denotes the mass of the kth star, G is the gravitational constant, and R − R_k a vector from a point in the disc to the primary or secondary star. From a numerical point of view, the smoothing factor εH is not necessary if the orbit of the binary is not inside the computational domain. We use this smoothing factor to account for the vertically extended three-dimensional disc in our two-dimensional case (Müller et al. 2012). For all simulations we use a value of ε = 0.6.

In our simulations with no bulk viscosity the viscous stress tensor Π is given in tensor notation by (5)where η = Σν is the vertically integrated dynamical viscosity coefficient. To model the viscosity in the disc we use the α-disc model by Shakura & Sunyaev (1973), where the kinematic viscosity is given by ν = αc_sH. The parameter α is less than one, and for our reference model we use α = 0.01.

To calculate the aspect ratio h = ^H/R we assume a vertical hydrostatic equilibrium (6)with the density ϱ and pressure p. To solve this equation we use the full three-dimensional potential of the binary and the isothermal equation of state (we also assume an isothermal disc in the z-direction). Integration over z yields (7)with the disc scale height (Günther & Kley 2002) (8)and the midplane density (9)which can be calculated by using the definition of the surface density.

If we concentrate the mass of the binary M_bin = M_A + M_B in its centre of mass (R_k → 0), this reduces to (10)with the Keplerian velocity . Test calculations with the simpler disc height (10)showed no significant difference compared to calculations with the more sophisticated disc height (8). Therefore, we use the simpler disc height in all our calculation to save some computation costs.

For our locally isothermal simulations we use a temperature profile of T ∝ R^-1, which corresponds to a disc with constant aspect ratio. The local sound speed is then given by c_s(R) = hu_K ∝ R^{− 1/ 2}. If not stated otherwise we use a disc aspect ratio of h = 0.05.

2.2. Initial disc parameters

For the initial disc setup we follow Lines et al. (2015). The initial surface density in all our models is given by (11)with the reference surface density Σ_ref = 10^-4M_⊙ au^-2. The initial slope is α_Σ = 1.5 and the gap function, which models the expected cavity created by the binary-disc interaction, is given by (12)(Günther & Kley 2002), with the transition width ΔR = 0.1 R_gap and the estimated size of the gap R_gap = 2.5 a_bin (Artymowicz & Lubow 1994). The initial radial velocity is set to zero u_R(t = 0) = 0, and for the initial azimuthal velocity we choose the local Keplerian velocity u_ϕ(t = 0) = u_K.

2.3. Initial binary parameters

For models run with Pluto or Fargo we start the binary at periastron at time t = 0 (upper signs in Eqs. (13)and (14)). Models carried out with Rh2d start the binary at apastron (lower signs in Eqs. (13)and (14)); with K₁ = M_B/M_bin, K₂ = M_A/M_bin, and the period of the binary .

2.4. Numerics

In our simulations we use a logarithmically increasing grid in the R-direction and a uniform grid in the ϕ-direction. The physical and numerical parameter for our reference system used in our extensive parameter studies in Sect. 5 and Sect. 6 are quoted in Table 2. In the radial direction the computational domain ranges from R_min = 0.25 au to R_max = 15.4 au and in the azimuthal direction from 0 to 2π, with a resolution of 762 × 582 grid cells. For our numerical studies we also use also different configurations as explained below.

Because of the development of an inner cavity where the surface density drops significantly, we use a density floor Σ_floor = 10^-9 (in code units) to avoid numerical difficulties with too low densities. Test simulations using lower floor densities did not show any differences.

At the outer radial boundary we implement a closed boundary condition where the azimuthal velocity is set to the local Keplerian velocity. At the inner radial boundary the standard condition is the zero-gradient outflow condition as in Kley & Haghighipour (2014) or Miranda et al. (2017), but we also test different possibilities such as closed boundaries or the viscous outflow condition (Kley & Nelson 2008; Mutter et al. 2017).

The open boundary is implemented in such a way that gas can leave the computational domain but cannot reenter it. This is done by using zero-gradient boundary conditions (∂/∂R = 0) for Σ and negative u_R. For positive u_R we use a reflecting boundary to prevent mass in-flow. Since there is no well-defined Keplerian velocity at the inner boundary, due to the strong binary-disc interaction, we also use a zero-gradient boundary condition for the angular velocity Ω_ϕ = u_ϕ/R. By using the zero-gradient condition for the angular velocity instead of the azimuthal velocity we ensure a zero-torque boundary. In the ϕ-direction we use periodic boundary conditions.

In all our simulations we use dimensionless units. The unit of length is R₀ = 1 au, the unit of mass is the sum of the primary and secondary mass M₀ = M_A + M_B and the unit of time is so that the gravitational constant G is equal to one. The unit density is then given by .

2.5. Monitored parameters

Since our goal is to study the evolution of the disc under the influence of the binary we calculated the disc eccentricity e_disc and the argument of the disc periastron ϖ_disc ten times per binary orbit. To calculate these quantities we treat each cell as a particle with the cells mass and velocity on an orbit around the centre of mass of the binary. Thus, the eccentricity vector of a cell is given by (15)with the specific angular momentum j = R × u of that cell. We have a flat system, thus the angular momentum vector only has a z-component. The eccentricity e_cell and longitude of periastron ϖ_cell of the cell’s orbit are therefore given by The global disc values are then calculated through a mass-weighted average of each cell’s eccentricity and longitude of periastron (Kley & Dirksen 2006) The integrals are simply evaluated by summing over all grid cells. The lower bound is always R₁ = R_min. For the disc eccentricity we integrate over the whole disc (R₂ = R_max) if not stated otherwise, whereas for the disc’s longitude of periastron it is suitable to integrate only over the inner disc (R₂ = 1.0 au) to obtain a well-defined value of ϖ_disc since animations clearly show a precession of the inner disc (see e.g. Kley & Haghighipour 2015). The radial eccentricity distribution of the disc is given by (20)

3. Hydrodynamic codes

Since the system under analysis is, in particular near the central binary, very dynamical we decided to compare the results from three different hydrodynamic codes to make sure that the observed features are physical and not numerical artefacts. We use codes with very different numerical approaches. Pluto solves the hydrodynamical equations in conservation form with a Godunov-type shock-capturing scheme, whereas Rh2d and Fargo are second-order upwind methods on a staggered mesh. In the following sections we describe each code and its features briefly.

3.1. Pluto

We use an in-house developed GPU version of the Pluto 4.2 code (Mignone et al. 2007). Pluto solves the hydrodynamic equations using the finite-volume method which evolves volume averages in time. To evolve the solution by one time step, three substeps are required. First, the cell averages are interpolated to the cell interfaces, and then in the second step a Riemann problem is solved at each interface. In the last step the averages are evolved in time using the interface fluxes.

For all three substeps Pluto offers many different numerical options. We found that the circumbinary disc model is very sensitive to the combination of these options. Third-order interpolation and time evolution methods lead very quickly to a negative density from which the code cannot recover, even though we set a density floor. We therefore use a second-order reconstruction of states and a second-order Runge-Kutta scheme for the time evolution. Another important parameter is the limiter, which is used during the reconstruction step to avoid strong oscillations. For the most diffusive limiter, the minmod limiter, no convergence is reached for higher grid resolutions. For the least diffusive limiter, the mc limiter, the code again produces negative densities and aborts the calculation. Thus, we use the van Leer limiter which, in kind of diffusion, lies between the minmod and the mc limiter.

For the binary position we solve Kepler’s equation using the Newton–Raphson method at each Runge-Kutta substep.

3.2. Rh2d

The Rh2d code is a two-dimensional radiation hydrodynamics code originally designed to study boundary layers in accretion discs (Kley 1989), but later extended to perform flat disc simulations with embedded objects (Kley 1999). It is based on the second-order upwind method described in Hawley et al. (1984) and Rozyczka (1985). It uses a staggered grid with second-order spatial derivatives and through operator-splitting the time integration is semi-second order. Viscosity can be treated explicitly or implicitly, artificial viscosity can be applied, and the Fargo algorithm has also been added. The motion of the binary stars is integrated using a fourth-order Runge-Kutta algorithm.

3.3. Fargo

We adopted the adsg version of Fargo (Masset 2000; Baruteau & Masset 2008) updated by Müller & Kley (2012). This code uses a staggered mesh finite difference method to solve the hydrodynamic equations. Conceptually, the Fargo code uses the same methods as Rh2d or the Zeus code, but employs the special Fargo algorithm that avoids the time step limitations that are due to the rotating shear flow (Masset 2000). In our situation the application of the Fargo algorithm is not always beneficial because of the larger deviations from pure Keplerian flow near the binary. The position of the binary stars is calculated by a fifth-order Runge-Kutta algorithm.

4. Numerical considerations

Before describing our results on the disc dynamics, we present two important numerical issues that can have a dramatic influence on the outcome of the simulations: the inner boundary condition (open or closed) and the location of the inner radius of the disc. We show that “unfortunate” choices can lead to incorrect results. In this section we use a numerical setup different from that in Table 2. Specifically, the base model has a radial extent from R_min = 0.25 au to R_max = 4.0 au, which is covered with 448 × 512 gridcells for simulations with Rh2d and FARGO, and 512 × 580 with Pluto (see Appendix A for an explanation of why we use different resolutions for different codes). For simulations with varying inner radii, the number of gridcells in the radial direction is adjusted to always give the same resolution in the overlapping domain.

4.1. Inner boundary condition

All simulations using a polar-coordinate grid expericence the same problem: there is a hole in the computational domain because R_min cannot be zero. Usually this hole exceeds the binary orbit. Therefore, the area where gas flows from the disc onto the binary and where circumstellar discs around the binary components form is not part of the simulation. This complex gas flow around the binary has been shown by e.g. Günther & Kley (2002) with a special dual-grid technique that covers the whole inner cavity. Their code is no longer in use and modern efficient codes (Fargo and Pluto) that run in parallel do not have the option of a dual-grid.

To reproduce nevertheless some results of Günther & Kley (2002), we carried out a simulation on a polar grid with an inner radius of R_min = 0.02 au so that both orbits of the primary and secondary lie well inside the computational domain. A snapshot of this simulation is plotted in Fig. 1. The logarithm of the surface density is colour-coded, while the orbits of the primary and secondary stars are shown in white and green. The red cross marks the centre of mass of the binary, which lies outside the computational domain. Figure 1 shows circumstellar discs around the binary components, as well as a complex gas flow through the gap onto the binary. Although including the binary orbit inside the computational would be desirable, it is at the moment not feasible for long-term simulations because of the strong time step restriction resulting from the small cell size at the inner boundary. This means that for simulations where the binary orbit is not included in the computational domain, the boundary condition at the inner radius has to allow for flow into the inner cavity, at least approximately. Two boundary conditions are usually used in circumbinary disc simulations: closed boundaries (Lines et al. 2015; Pierens & Nelson 2013) and outflow boundaries (Kley & Haghighipour 2014, 2015; Miranda et al. 2017).

Given the importance of the boundary condition at R_min, we carried out dedicated simulations for a model adapted from Lines et al. (2015) and used an inner radius of R_min = 0.345 au, but otherwise it was identical to our standard model. We applied closed boundaries and open ones in order to examine their influence on the disc structure. We also varied different numerical parameters (resolution, radial grid spacing, integrator) for this simulation series.

The top panel of Fig. 2 shows the azimuthally averaged density profiles for different grid resolutions and spacings after 16 000 binary orbits. All these simulations were performed with Pluto and a closed inner boundary. The strong dependence of the density distribution and inner gap size on the numerical setup stands out. This dependence on numerics is very surprising and not at all expected since the physical setup was identical in all these simulations and therefore the density profiles should all be similar and converge upon increasing resolution. Not only does the gap size show this strong dependence, but also the radial eccentricity distribution (bottom panel of Fig. 2). To examine this dependence further, which could in principle be the result of a bug in our code, we reran the identical physical setup with a different code, Rh2d. The results of these Rh2d simulations show the same strong numerical dependence as the Pluto results. In Fig. 3 the total disc eccentricity time evolution for the Rh2d simulations is plotted. Again, one would expect that different numerical parameters should produce approximately the same disc eccentricity, but our simulations show a radical change in the disc eccentricity if the numerical methods are slightly different.

We note that using a viscous outflow condition where the radial velocity at R_min is fixed to − β3ν/ (2R) with β = 5 (as suggested by Mutter et al. 2017) resulted in the same dynamical behaviour as the closed boundary condition. The reason for this lies in the fact that the viscous speed is very low in comparison to the radial velocity induced by the perturbations of the central binary and hence there is little difference between a closed and a viscous boundary.

In contrast to the closed inner boundary simulations, simulations with a zero-gradient outflow inner boundary reach a quasi-steady state. For the open boundaries Pluto simulations with different numerical parameters produce very similar surface density and radial disc eccentricity profiles (Fig. 4) for different resolutions and grid spacings. Only the values for the uniform radial grid with a resolution of 395 × 512 (purple line in Fig. 4) deviate. Here the resolution in the inner computational domain is not high enough, because a simulation with a higher resolution uniform grid (790 × 512, orange line in Fig. 4) produces approximately the same results as the simulations with a logarithmic radial grid spacing.

We do not have a full explanation for this strong variability and non-convergence of the flow when using the closed inner boundary, but this result seems to imply that this problem is ill-posed. A closed inner boundary creates a closed cavity, which apparently implies in this case that the details of the flow depend sensitively on numerical diffusion as introduced by different spatial resolutions, time integrators, or codes. As a consequence, for circumbinary disc simulations a closed inner boundary is not recommended for two reasons. Firstly, it produces the described numerical instabilities and secondly, for physical reasons the inner boundary should be open because otherwise no material can flow into the inner region and be accreted by the stars. This is also indicated in Fig. 1 which shows mass flow through the inner gap, which cannot be modelled with a closed inner boundary.

4.2. Location of the inner radius

After determining that a zero gradient outflow condition is necessary, we use it now in all the following simulations and investigate in a second step the optimal location of R_min. It should be placed in such a way that it simultaneously insures reliable results and computational efficiency. For the excitation of the disc eccentricity through the binary-disc interaction, non-linear mode coupling, and the 3:1 Lindblad resonance are important (Pierens & Nelson 2013). Therefore, the location of the inner boundary is an important parameter and R_min should be chosen such that all major mean-motion resonances between the disc and the binary lie inside the computational domain. To investigate the influence of the inner boundary position, we set up simulations with an inner radius from R_min = 0.12 au to R_min = 0.5 au for the standard system parameter as noted in Table 2. Figure 5 shows the results of these simulations. Displayed are the azimuthally averaged surface density and the radial disc eccentricity distribution after 16 000 binary orbits. In agreement with Kley & Haghighipour (2014), we find that the radial disc eccentricity does not depend much on the location of the inner radius and remains low in the outer disc (R>2 au). Only for a large inner radius of R_min = 0.5 au do we observe almost no disc eccentricity growth (brown line in Fig. 5). In this case the 3:1 Lindblad resonance (R_3:1 = 0.46 au) lies outside the computational domain, confirming the conclusion from Pierens & Nelson (2013) about the importance of the 3:1 Lindblad resonance for the disc eccentricity growth. The variation of the radial disc eccentricity for inner radii smaller than 0.3 au occurs because the precessing discs are in different phases after 16 000 binary orbits. Contrary to the eccentricity distribution, the azimuthally averaged surface density shows a stronger dependence on the inner radius. As seen in the top panel of Fig. 5, the maximum of the surface density increases monotonically as the location of the inner boundary decreases because a smaller inner radius means also that the “area” where material can leave the computational domain decreases. For inner radii smaller than R_min≤0.25 au the change in the maximum surface density becomes very low upon further reduction of R_min. This implies that for too large R_min there is too much mass leaving the domain and it has to chosen small enough to capture all dynamics. Although the maximum value of the surface density increases, the position of the maximum does not depend on the inner radius and remains at approximately R = 1.1 au as long as the 3:1 Lindblad resonance lies inside the computational domain. Furthermore, we find that the slope of the surface density profile at the gap’s edge does not depend on the location of the inner boundary. These results are in contrast to Mutter et al. (2017) who find a stronger dependence of the density profile on the location of the inner radius. However, their results may be a consequence of the viscous outflow condition used.

As is discussed below, the disc becomes eccentric with a slow precession that depends on the binary eccentricity, e_bin. To explore how the location of the inner radius affects the disc dynamics, we ran additional simulations with varying R_min for different e_bin. We chose e_bin = 0.08 and e_bin = 0.32 in addition to the Kepler-16 value of e_bin = 0.16, which lies near the bifurcation point, e_crit. Figure 6 shows the precession period of the inner gap for varying inner radii and three different binary eccentricities, where the blue curve refers to the model shown above with e_bin = 0.16. For a higher binary eccentricity of e_bin = 0.32 (green curve in Fig. 6 and on the upper branch of the bifurcation diagram) the disc dynamics seems to be captured well even for higher values of R_min. Only for values lower than R_min = 0.22 au are deviations seen. The simulation with R_min = 0.20 au was more unstable, probably because the secondary moved in and out of the computational domain on its orbit. On the lower branch of the bifurcation diagram the convergence with decreasing R_min is slower as indicated by the case e_bin = 0.08 (red curve in Fig. 6). Here, an inner radius of R_bin = a_bin or even slightly lower may be needed. One explanation for this behaviour is that on the lower branch (for low e_bin) the inner gap is smaller but nevertheless quite eccentric such that the disc is influenced more by the location of the inner boundary. The convergence of the results for smaller inner radii is also visible in the azimuthally averaged surface density and radial eccentricity profiles for e_bin = 0.08 and e_bin = 0.32 displayed in Fig. 7.

In summary, from the results shown in Figs. 5–7 we can conclude that to model the circumbinary disc properly, the inner radius should correspond to R_min ≈ a_bin. This should be used in combination with a zero-gradient outflow (hereafter just outflow) condition. Since a smaller inner boundary increases the number of radial cells and imposes a stricter condition on the time step, long-term simulations with a small inner radius are very expensive. As a compromise to make long-term simulations possible, we have chosen in all simulations below an inner radius of R_min = 0.25 au.

4.3. Location of the outer radius

While investigating the numerical conditions at the inner radius, we used an outer radius of R_max = 4.0 au and assumed that this value is high enough, so that the outer boundary would not interfere with the dynamical behaviour of the inner disc. During our study of different binary and disc parameters, we found that in some cases this assumption is not true. Especially for binary eccentricities greater than 0.32, reflections from the outer boundary interfered with the complex inner disc structure. Therefore, we increased the outer radius of all the following simulations to R_max = 15.4 au = 70 a_bin, a value used by Miranda et al. (2017). In order to keep the same spatial resolution as before we also increased the resolution to 762 × 582 cells. This ensures that in the radial direction we still have 512 grid cells between 0.25 au and 4.0 au. To save computational time it is also possible to use a smaller outer radius (R_max ≈ 40 a_bin) with a damping zone where the density and radial velocity are relaxed to their initial value. A detailed description of the implementation of such a damping zone can be found in de Val-Borro et al. (2006).

Hence, for our subsequent simulations we use from now on the numerical setup quoted in Table 2, unless otherwise stated. To study the dependence on different physical parameter we start from the standard values of the Kepler-16 system in Table 2 and vary individual parameter separately.

5. Variation of binary parameters

In this section we study the influence of the central binary parameter specifically its orbital eccentricity and mass ratio on the disc structure. Throughout this section we use a disc aspect ratio of h = 0.05 and a viscosity α = 0.01. For the inner radius and the inner boundary condition we use the results from the previous section and apply the outflow condition at R_min (see also Table 2). All the results in this section were obtained with Pluto, but comparison simulations using Rh2d show very similar behaviour (see Fig. A.4).

5.1. Dynamics of the inner cavity

Before discussing in detail the impact of varying e_bin and q_bin, we comment first on the general dynamical behaviour of the inner disc. Figure 8 shows snapshots of the inner disc after 16 000 binary orbits for different binary eccentricities. As seen in several other studies (as mentioned in the introduction) we find that the inner disc becomes eccentric and shows a coherent slow precession. In the figure we overplot ellipses (white dashed lines) that are approximate fits to the central inner cavity, where the semi-major axis (a_gap) and eccentricity (e_gap) are indicated. To calculate these parameters we assumed first that the focus of the gap-ellipse is the centre of mass of the binary. We then calculated from our data the coordinates (R_Σmax,ϕ_Σmax) of the cell with the highest surface density value, which defines the direction to the gap’s apocentre. We defined the apastron of the gap R_apa as the minimum radius along the line (R,ϕ_Σmax) which fulfils the condition (21)The periastron of the gap R_peri is defined analogously as the minimum radius along the line (R,ϕ_Σmax + π) in the opposite direction which fulfils (22)Using the apastron and periastron of the gap as defined above, the eccentricity and semi-major axis of the gap are given by As shown in Fig. 8 this purely geometrical construction matches the shape of the central cavity very well. Even though the overall disc behaviour shows a rather smooth slow precession, the dynamical action of the central binary is visible as spiral waves near the gap’s edges, most clearly seen in the first and last panel. We prepared some online videos to visualise the dynamical behaviour of the inner disc.

An alternative way to characterise the disc gap dynamically is by using the orbital elements (e_max,a_max) of the cell with the maximum surface density. These orbital elements can be calculated with Eq. (15)and the vis-viva equation (25)The time evolution of these gap characteristics are plotted in Fig. 9. The gap’s size and eccentricity show in phase oscillatory behaviour with a larger amplitude in the eccentricity variations. The period of the oscillation is identical to the precession period of the disc. Clearly, the extension of the gap always lies inside of the location of maximum density a_gap<a_max, but for the eccentricities the ordering is not so clear. The radial variations of the disc eccentricity as shown in Figs. 4 and 5 indicate that the inner regions of the disc typically have a higher eccentricity. In Fig. 9e_disc is lowest because it is weighted with the density which is very low in the inner disc regions. The eccentricity for the gap e_gap stems from a geometric fit to the very inner disc regions and is the highest. Inside the maximum density and even slightly beyond that radius the disc precesses coherently in the sense that the pericentres estimated at different radii are aligned. The data show further that when the disc eccentricity is highest the disc is fully aligned with the orbit of the secondary star, i.e. the pericentre of disc and binary lie in the same direction (see also Appendix A).

5.2. Binary eccentricity

For simulations in this section we fixed q_bin = 0.29 and varied e_bin from 0.0 to 0.64. The binary eccentricity strongly influences the size of the gap as well as the disc precession period of the inner disc.

In Fig. 10 the azimuthally averaged surface density profiles are shown for the various e_bin at 16 000 binary orbits. In all cases there is a pronounced density maximum visible which is the strongest for the circular binary with e_bin = 0 (red curve). The position of the peak varies systematically with binary eccentricity. Increasing e_bin from zero to higher values the peak shifts inward. Not only does the gap size decrease, but the maximum surface density also drops with increasing e_bin until a minimum at e_bin = 0.18 is reached. Increasing e_bin further causes the density peak to move outward again, whereas the maximum surface density stays roughly constant for higher binary eccentricities.

As mentioned above, the eccentric disc precesses slowly and in Fig. 11 the longitude of the inner disc’s pericentre is shown versus time over 16 000 binary orbits. The disc displays a slow prograde precession with inferred period of several thousand binary orbits. For the measurement of the precession period we discarded the first 6000 binary orbits since plots of the longitude of periastron show that during this time span the precession period is not constant yet (see Fig. 11). The period is then calculated by averaging over at least two full periods. To be able to average over two periods, models with a very long T_prec (e.g. e_bin = 0.64) were simulated for more than 16 000 binary orbits.

Figure 12 summarises the results from simulations with varying binary eccentricities. The top panel of Fig. 12 shows different measures for the gap size for varying binary eccentricities. In addition to the radial position where the azimuthally averaged surface density reaches its maximum (R_peak), we plot the positions where the density drops to 50 and 10 percent of its maximum value (R_0.5 and R_0.1). The value R_0.5 was used by Kley & Haghighipour (2014) as a measure for the gap size, whereas Miranda et al. (2017) used R_0.1 and Mutter et al. (2017)R_peak. Starting from a non-eccentric binary the curves for R_peak and R_0.5 decrease with increasing binary eccentricity. For e_bin ≈ 0.18 the gap size reaches a minimum and then increases again for higher binary eccentricities. In agreement with Miranda et al. (2017) we see an almost monotonic increase of R_0.1 for increasing binary eccentricities. The gap size, a_gap, correlates well with R_0.5 and is always about 14% smaller.

The middle panel of Fig. 12 shows the eccentricity of the inner cavity, calculated with the method described in Sect. 5.1. The disc eccentricity also changes systematically with e_bin. For circular binaries it reaches e_gap ≈ 0.44, and then it drops down to about 0.25 for the turning point e_bin = 0.18, and increases again reaching e_gap ≈ 0.4 for the highest e_bin = 0.64. Hence, for nearly circular binaries e_gap can be higher than for more eccentric binaries. A similar variation of the disc’s eccentricity and precession rate with binary eccentricity has been noticed by Fleming & Quinn (2017) who attribute this to the possible excitation of higher order resonances for higher disc eccentricities.

The bottom panel of Fig. 12 shows the precession period of the inner disc (R<1 au) for different binary eccentricities. The curve for the precession period shows a similar behaviour to the upper curves for the gap size. To investigate further the correlation of the gap size (here represented by R_0.5) and the precession period, we show in Fig. 13 the two quantities plotted against each other. Two branches can be seen as indicated by the dashed curves. One starts at e_bin = 0.0 and goes to e_bin = 0.18 where the gap size and precession period decrease with increasing binary eccentricities, and the other branch starts at the minimum at e_bin = 0.18 where both gap properties increase again with increasing binary eccentricities. These two branches may indicate two different physical processes that are responsible for the creation of the eccentric inner gap and show the direct correlation between precession period and gap size.

In our simulations all discs became eccentric and showed a prograde precession, and we did not find any indication for “stand still” discs that are eccentric but without any precession. In contrast, Miranda et al. (2017) find for their disc setup (h = 0.1,α = 0.1) that the disc does not precess for binary eccentricities between 0.2 and 0.4. We do not find this behaviour for our disc models, but note that we use a disc with a lower aspect ratio and a lower α-value. To investigate this further we set up a simulation with the same numerical parameters as Miranda et al. (2017) and used a physical setup where they did not observe a precession of the disc (q_bin = 1.0, e_bin = 0.4, h = 0.1, and α = 0.1). We carried out two simulations, one with an inner radius of R_min = (1 + e_bin)a_bin and one with an inner radius of R_min = 1.136 a_bin. The first radius was used by Miranda et al. (2017) and in this case we also observe no precession of the disc. The second radius corresponds to R_min = 0.25 au, which is the radius we established in Sect. 4.2 as the optimum location of the inner boundary. In the case with the smaller inner radius we observe a clear precession of the inner disc (Fig. 14) with a relatively short precession period T_prec as expected for high viscosity and high h (see below). This is another indication of the necessity of choosing R_min sufficiently small to capture all effects properly.

5.3. Binary mass ratio

In this section we study the influence of the binary mass ratio. We carried out simulations with a mass ratio from q_bin = 0.1 to q_bin = 1.0; all other parameters were set according to Table 2. Figure 15 shows the azimuthally averaged density profiles for varying binary mass ratios. The density profiles do not show such a strong variation as in simulations with varying binary eccentricities. Only the density profiles of the two lowest simulated mass ratios q_bin = 0.1 and q_bin = 0.2 differ significantly from the other profiles. For these models the position of the peak surface density is closer to the binary and the maximum surface density roughly 20 percent lower. The different gap size measures, introduced in the previous section, are shown in the top panel of Fig. 16. In general, the variations of all three gap-size indicators with mass ratio q_bin are relatively weak. All three gap size measurements remain nearly constant for mass ratios greater than 0.3. Since the density profiles do not show a significant variation this is not surprising. At the same time the size and eccentricity of the gap do not vary significantly and lie in the range a_gap ≈ 4 a_bin and e_gap ≈ 0.27. Overall, compared to simulations with varying binary eccentricity the binary mass ratio does not have such a strong influence on the inner disc structure for q_bin ≳ 0.3.

In contrast, the precession period of the inner disc shows a far stronger dependence on the binary mass ratio than the gap size (bottom panel of Fig. 16). Discs around binaries with a higher mass ratio have a lower precession period than discs around binaries with a low mass ratio. This dependence can be understood in terms of free particle orbits around a binary star that also display a reduction of the precession period with increasing q, as we discuss in more detail in Appendix B.

6. Variation of disc parameters

In this section we explore the influence of discs parameters, namely the aspect ratio and the alpha viscosity, on the inner disc structure. We use the Kepler-16 values for the binary and our standard numerical setup (Table 2).

While performing the simulations for this paper we observed that models with low pressure (low H/R) and high viscosity (high α) seem to be very challenging for the Riemann-Code Pluto. As a consequence we also present in this section results using the Upwind-Code Rh2d. As discussed in Appendix A, the two codes produce results in very good agreement with each other for our standard model. For other parameters they do not necessarily produce results with such good agreement, but the simulation results from both codes always show the same trend. Therefore, we decided to mix simulation results from Pluto and Rh2d to cover a broader parameter range.

6.1. Disc aspect ratio

First we varied the disc aspect ratio h from 0.03 to 0.1. Our simulation results show a decreasing gap size for higher aspect ratios (Fig. 17). The gaps precession period is again directly correlated to the gap size, an observation we have already seen for different binary eccentricities. In general, a drop in T_prec with higher h is expected because an increase in pressure will tend to reduce the gap size, which will in turn lead to a faster precession. This trend is indeed observed in our simulations for higher h. For simulations with h< 0.05 we observed a decrease in gap size and precession period. The drop in T_prec with lower values of h may be partly due to the lack of numerical resolution and partly to the lack of pressure support, which may no longer allows for coherent disc precession.

6.2. Alpha viscosity

In this section we discuss the influence of the magnitude of viscosity on the disc structure and therefore set up simulations with different viscosity coefficients ranging from α = 0.001 to α = 0.1. We then analysed the structure and behaviour of the inner cavity as before. The results are summarised in Fig. 18. A clear trend is visible for the gap size and for the precession period of the gap. For higher viscosities the gap size shrinks and the precession period decreases. Again, we see the direct correlation of gap size and precession rate, as already observed for aspect ratio and binary eccentricity variations. The decrease in the disc’s gap size can be explained: for higher α values the viscous spreading of the disc increases. This viscous spreading counteracts the gravitational torques of the binary, which are responsible for the gap creation.

7. Summary and discussion

Using two-dimensional hydrodynamical simulations, we studied the structure of circumbinary discs for different numerical and physical parameters. Since these simulations are, from a numerical point of view, not trivial and often it is not easy to distinguish physical features and numerical artefacts, we checked our results with three different numerical codes. In the following we summarise the most important results from our simulations with different numerical and physical parameters.

Concerning the numerical treatment the most crucial issue with simulations using grid codes in cylindrical coordinates is the treatment of the inner boundary condition. As there has been some discussion in the literature about this issue (Pierens & Nelson 2013; Kley & Haghighipour 2014; Lines et al. 2015; Miranda et al. 2017; Mutter et al. 2017), we decided to perform a careful study of the location of the inner boundary and the conditions on the hydrodynamical variables at R_min. Our results can be summarised as follows:

• Inner boundary condition. Since there was no agreement onwhich type of inner boundary condition should be used for cir-cumbinary disc simulations, especially if closed or open bound-ary conditions should be used, we investigated the influence of theinner boundary condition through dedicated simulations withtwo different codes. We observed that closed boundaries can leadto numerical instabilities for all codes used in our comparison andno convergence to a unique solution could be found. The use of aviscous outflow condition where the velocity is set to the meanviscous disc speed produced very similar results to the closedboundary because the viscous speed is very low in comparison tothe velocities induced by the dynamical action of the binary stars.Open inner boundaries, on the other hand, lead to numericallystable results and are also, from a physical point of view, morelogical because material can flow onto the binary and be accretedby one of the stars. Hence, open boundaries with free outflow arethe preferred conditions at R_min.

• Location of inner boundary. The location of the inner boundary is also an important numerical parameter. Since an inner hole in the computational domain cannot be avoided in polar coordinates, the inner radius has to be sufficiently small to capture all relevant physical effects. In particular, all important mean motion resonances, which may be responsible for the development of the eccentric inner cavity, should lie inside the computational domain. Through a parameter study we were able to determine that the radius of the inner boundary R_min has to be of the order of the binary separation a_bin to capture all physical effects in a proper way.

After having determined the best values for the numerical issues we performed, in the second part of the paper a careful study to determine the physical aspects that determine the dynamics of circumbinary discs. Starting from a reference model, where we chose the binary parameter of Kepler-16, we varied individual parameters of the binary and the disc and studied their impact on the disc dynamics. First of all, for all choices of our parameters the models produce an eccentric inner disc that shows a coherent prograde precession. However, the size of the gap, its eccentricity, and its precession rate depend on the physical parameters of binary and disc. Our results can be summarised as follows:

• Binary parameters. To study the influence of the binary star on thedisc we systematically varied the eccentricity(e_bin) and the mass ratio (q_bin) of the binary. The parameter study showed that e_bin has a strong influence on the gap size and on the precession period of the gap. We found that two regimes exist where the disc behaves differently to an increasing binary eccentricity (see Fig. 8). The two branches bifurcate at a critical binary eccentricity, e_crit = 0.18 from each other. From e_bin = 0.0 to e_crit the gap size and precession period decrease, and from e_bin = 0.18 onward both gap parameters become larger again, as displayed in Fig. 13. The bifurcation of the two branches near e_crit = 0.18 strongly suggests that different physical mechanisms, responsible for the creation of the eccentric inner cavity, operate in the two regimes. For the lower branch the excitation at the 3:1 outer Lindblad resonance may be responsible, which is supported by a simulation where the inner boundary was outside the 3:1 radius and no disc eccentricity was found. For the upper branch (high e_bin) non-linear effects may be present, as suggested by Fleming & Quinn (2017) who found similar behaviour. As second binary parameter we varied the mass ratio of the secondary to the primary star, q_bin, between [ 0.1,1.0 ] and studied its impact on the gap size and precession period. Overall the variation of R_0.5 with q_bin is weak. For low q_bin the gap size increases until it becomes nearly constant for q_bin ≳ 0.3 with R_0.5 ≈ 4.6 a_bin, as shown in Fig. 16. The precession period decreases on average with increasing T_prec and for large q_bin the behaviour is equivalent to a single particle at a separation of 4.9 au, as we discuss in more detail in Appendix B.

• Disc parameters. We also varied different disc parameters, namely the pressure (through the aspect ratio h = H/R) and the viscosity (through α). The results are displayed in Figs. 17 and 18. For changes in the aspect ratio no clear trend was visible; both gap size and precession period first increase with increasing h until they reach a maximum at h = 0.05; from there both gap properties decrease again with increasing aspect ratio. For high h the behaviour can be understood in term of the gap closing tendency of higher disc pressure. On the other hand, for very low pressure it may be more difficult to sustain a coherent precession of a large inner hole in the disc due to the reduced sound speed. Additionally, the damping action of the viscosity becomes more important for discs with lower sound speed. A clear monotonic trend was visible for the viscosity variation, where higher α leads to smaller gaps and shorter gap precession periods. This can be directly attributed to the gap closing tendency of viscosity.

To allow for an alternative view of the impact of different binary and disc parameters on the dynamical structure of the disc we display in Fig. 19 the eccentricity (top) and the precession period (bottom) of the gap plotted against the semi-major axis of the gap for all our simulations where different colours stand for different model series. Pluto simulations are represented by dots whereas squares stand for Rh2d simulations. The reference model is marked with a star. Interestingly the majority of the models lie on a main-sequence branch, which has the trend of increasing eccentricity and precession period with increasing gap width. The smaller branch that bifurcates near the reference model represents the low e_bin sequence, which has higher e_gap but smaller T_prec. By coincidence, our reference model with the Kepler-16 parameters lies near the bifurcation point.

How the location of the bifurcation point depends of the disc and binary parameter needs to be explored in subsequent numerical studies. A first simulation series with q_bin = 0.6 suggest that this critical eccentricity does not depend on the mass ratio of the binary. As outlined in the Appendix, for parameters very near to the critical e_bin the models can have the tendency to switch between the two states during a single simulation. The physical cause of the bifurcation needs to be analysed in subsequent simulations.

In addition, there are a more possibilities to further improve our model. First, the use of a non-isothermal equation should be considered to study the effects of viscous heating and radiative cooling. Especially for simulations, which also evolve planets inside the disc, Kley & Haghighipour (2014) showed that radiative effects should been taken into account. Self-gravity effects have been studied by Lines et al. (2015) and Mutter et al. (2017) and are important for galactic discs around a central binary black hole. The backreaction of the disc onto the binary could also be considered in more detail. Angular momentum exchange with the binary occurs through gravitational torques from the disc and direct accretion of angular momentum from the material that enters the central cavity. The first part always leads a decrease in the binary semi-major axis and an increase in eccentricity (Kley & Haghighipour 2015), while Miranda et al. (2017) pointed out that angular momentum advection might even be dominant leading to a binary separation.

Another interesting question is the evolution of planets inside the circumbinary disc and their influence on the disc structure. Since an in situ formation of the observed planets seems very unlikely, the migration and final parking position of planets is an interesting topic that requires further investigations. In view of our extensive parameter studies there is the indication that in several past studies (Pierens & Nelson 2013; Kley & Haghighipour 2014, 2015) an inner radius was considered that was too large and that did not allow for full disc dynamics.

Through fully three-dimensional simulations it will be possible to study the dynamical evolution of inclined discs.

---

Acknowledgments

Daniel Thun was funded by grant KL 650/26 of the German Research Foundation (DFG), and Giovanni Picogna acknowledges DFG support grant KL 650/21 within the collaborative research program “The first 10 Million Years of the Solar System”. Most of the numerical simulations were performed on the bwForCLuster BinAC, supported by the state of Baden-Württemberg through bwHPC, and the German Research Foundation (DFG) through grant INST 39/963-1 FUGG. All plots in this paper were made with the Python library matplotlib (Hunter 2007).

References

Appendix A: Convergence study

The strong gravitational influence of the binary onto the disc, and the complex dynamics of the inner disc presents a big challenge for numerical schemes of grid codes. We therefore study in this appendix the dependence of the disc structure on numerical parameters, such as the numerical scheme or the grid resolution. First, we simulated our reference system, using an open inner boundary and an inner radius of R_min = 0.25 au, with the different codes described in Sect. 3.

In Fig. A.1 the long-time evolution of the disc eccentricity (top) and of the longitude of periastron (bottom) are shown. All three codes produce comparable results, although we should note that to get this agreement simulations performed with Pluto need a slightly higher resolution. Pluto simulations with a resolution of 448 × 512 grid cells produce a slightly smaller precession period of the inner disc (red line in Fig. A.3, see also Table A.1). One possible explanation for this could be the size of the numerical stencil. Fargo and Rh2d have a very compact stencil because they use a staggered grid. In contrast, Pluto has a wider stencil because a collocated grid is used, where all variables are defined at the cell centres. In particular the viscosity routine of Pluto has a very wide stencil to achieve the correct centring of the viscous stress tensor components. This could explain why for Pluto a slightly higher resolution is needed.

The π-shift of the longitude of periastron calculated with Rh2d compared to the values calculated with Pluto and Fargo (bottom panel of Fig. A.1) can be explained in the following way. Simulations performed with Pluto and Fargo started the binary at periastron with the pericentre of the secondary at ϕ = π (see red ellipses in Fig. 8), whereas Rh2d started the binary at apastron with the pericentre of the secondaries orbit at ϕ = 0. The data show that when the disc eccentricity is at its maximum the disc is aligned with the orbit of the secondary star. Since the pericentre of the secondary in Rh2d simulations is shifted by π we can also see this π-shift in the disc’s longitude of periastron.

Figure A.2 shows the azimuthally averaged density profiles at t = 16 000 T_bin. Even after such a long time span all three codes agree very closely. All codes produce roughly the same density maximum and the same density slopes.

To check the influence of the numerical resolution on the disc structure we run Pluto simulations with different resolutions, summarised in Table A.1. The results of this resolution test are shown in Fig. A.3. The disc eccentricity converges to the same value, whereas it seems that the disc precession rates increases with higher resolution. A closer look at the data shows that after around 6000 orbits the precession rate for the higher resolution converged. Table A.1 summarises the measured precession rate for different resolutions. Since there is no large variation between 512 × 580 and 600 × 1360 grid cells, we used a resolution of 512 × 580 grid cells for our Pluto simulations.

As already discussed in Sect. 4.2 in some cases (for example high binary eccentricities) an outer radius of R_max = 4.0 au can lead to wave reflections which interfere with the inner disc. Therefore we increased the outer boundary to R_max = 15.4 au = 70 a_bin (and to ensure the same resolution between 0.25 au and 4.0 au we also increased the number of grid points to 762 × 582). Figure A.4 compares the effect of different outer radii for the case of varying binary eccentricities. For binary eccentricities greater than 0.32, a larger outer radius makes a difference. Also, around the bifurcation point, at e_bin = 0.18, a larger outer radius changes the gap size and precession rate of the disc. We also added results from Rh2d simulations with an outer radius of R_max = 18.18 a_bin. These Rh2d results follow the same trend as the Pluto simulations, although for most binary eccentricities we could not reproduce the very good agreement of the e_bin = 0.16 case. Around e_bin = 0.18, which is the critical binary eccentricity where the two branches bifurcate, simulations tend to switch between two states. This jump happens after roughly 3000 binary orbits and does not depend on physical parameters. Changes in numerical parameters, like the Courant-Friedrichs-Lewy condition or R_max, can trigger the jump. This can be seen in Fig. A.4 for e_bin = 0.16 where simulations with R_max = 18.18 a_bin (red curve) and R_max = 70.0 a_bin (blue curve) produce different values for the gap size and the precession period. An outer radius of R_max = 18.18 a_bin is in the case of e_bin = 0.16 not too small, since a simulation with R_max = 100 a_bin again produced the same values for the gap size and the precession period. For simulations with e_bin far away from e_crit, we did not observe these jumps between states. So far we have only observed these jumps in simulations carried out with Pluto.

Appendix B: Comparison with massless test particle

In our circumbinary disc evolution we found that the inner disc becomes eccentric with a constant precession rate. In this section we compare this to test particle trajectories around the binary and investigate how a massless particle with orbital elements similar to the disc gap (Fig. 19) would behave under the influence of the binary potential. In particular, we would like to know at what distance from the binary a test particle has to be positioned such that its precession period agrees approximately with that of the disc. Using secular perturbation theory for the coplanar motion of a massless test particle around an eccentric binary star Moriwaki & Nakagawa (2004) derived the following formula for the particle’s precession period for low binary eccentricities (B.1)where a_p is the semi-major axis of the particle. The same relation was quoted by Miranda et al. (2017) based on results by Liu et al. (2015). A similar relation was given by Leung & Lee (2013) without the e_bin term.

In order to verify the applicability of expression (B.1) when varying of binary and planet parameters, we performed series of three-body simulations with very low planet mass (10^-6M_⊙), where we varied individually q_bin, e_bin, and a_p. In addition we varied the planet eccentricity e_p. For the simulations we used the parameter of the Kepler-16 systems as a reference (see Table 1) and integrated the system for several 10 000 binary orbits. Our results of these three-body simulations are shown in Fig. B.1. We find that the precession rate scales with q_bin, a_p, and e_bin exactly as expected from relation (B.1). The fourth panel of Fig. B.1 indicates that the precession period also depends on the particle eccentricity as , where we plotted the average eccentricity of the orbit which is equivalent to the particle’s free eccentricity. The agreement holds up to about for the used a_p. Clearly, for higher values of e_p the particle’s orbit will be unstable as this leads to close encounters with the binary. The value of e_p where this happens depends on the distance from the binary star. For lower a_p the range will be more limited. In summary, we find for the precession period of a particle around a binary star (B.2)This relation is plotted in Fig. B.1 as the solid black line. In addition to the scaling behaviour we find that the numerical results are about 3 − 5% lower than the theoretical estimate. The agreement of the numerical results with the theoretical prediction becomes better for larger a_p because Eq. (B.2) is derived from an approximation for large a_p/a_bin. The full relation (B.2) can be inferred directly from Eq. (11) in Georgakarakos & Eggl (2015).

Now we compare the particle precession rate to our disc simulations. The best option for this comparison is to check the scaling of the precession rate for models where q_bin has been varied because for binary mass ratio greater than q_bin = 0.3 the gap size and eccentricity do not change much (blue points in Fig. 19). This is supported by the small variation of the gap radii in Fig. 16. For higher mass ratios the gap size is approximately a_gap = 4.0 a_bin and e_gap = 0.27. As seen from Eq. (B.2), a test particle with these orbit elements would have a precession period that is far too short. However, if we assume that the test particle has a semi-major axis of a_p = 4.9 a_bin (roughly 20 percent more than a_gap), the precession period of the gap and the particle match very well for different q_bin (Fig. B.2), at least for the higher mass ratios. In general, however, due to the strong sensitivity of the precession period with a_p, an exact agreement will be difficult to obtain.

Movies

Movie of Fig. 1 Access here

Movie 1 of Fig. 8 Access here

Movie 2 of Fig. 8 Access here

Movie 3 of Fig. 8 Access here

Movie 4 of Fig. 8 Access here

All Tables

All Figures