Circumbinary discs: Numerical and physical behaviour
AAstronomy & Astrophysics manuscript no. cb © ESO 2018June 14, 2018
Circumbinary discs: Numerical and physical behaviour
Daniel Thun , Wilhelm Kley , and Giovanni Picogna Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, D-72076 Tübingen, Germanye-mail: {daniel.thun@, wilhelm.kley@}uni-tuebingen.de Universitäts-Sternwarte, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 München, Germanye-mail: [email protected]
ABSTRACT
Aims.
Discs around a central binary system play an important role in star and planet formation and in the evolution of galactic discs.These circumbinary discs are strongly disturbed by the time varying potential of the binary system and display a complex dynamicalevolution that is not well understood. Our goal is to investigate the impact of disc and binary parameters on the dynamical aspects ofthe disc.
Methods.
We study the evolution of circumbinary discs under the gravitational influence of the binary using two-dimensional hydro-dynamical simulations. To distinguish between physical and numerical e ff ects we apply three hydrodynamical codes. First we analysein detail numerical issues concerning the conditions at the boundaries and grid resolution. We then perform a series of simulationswith di ff erent binary parameters (eccentricity, mass ratio) and disc parameters (viscosity, aspect ratio) starting from a reference modelwith Kepler-16 parameters. Results.
Concerning the numerical aspects we find that the length of the inner grid radius and the binary semi-major axis must becomparable, with free outflow conditions applied such that mass can flow onto the central binary. A closed inner boundary leads tounstable evolutions.We find that the inner disc turns eccentric and precesses for all investigated physical parameters. The precession rate is slow withperiods ( T prec ) starting at around 500 binary orbits ( T bin ) for high viscosity and a high aspect ratio H / R where the inner hole is smallerand more circular. Reducing α and H / R increases the gap size and T prec reaches 2500 T bin . For varying binary mass ratios q bin the gapsize remains constant, whereas T prec decreases with increasing q bin .For varying binary eccentricities e bin we find two separate branches in the gap size and eccentricity diagram. The bifurcation occurs ataround e crit ≈ .
18 where the gap is smallest with the shortest T prec . For e bin lower and higher than e crit , the gap size and T prec increase.Circular binaries create the most eccentric discs. Key words.
Hydrodynamics – Methods: numerical – Planets and satellites: formation – Protoplanetary discs – Binaries: close
1. Introduction
Circumbinary discs are accretion discs that orbit a binary sys-tem that consists for example of a binary star or a binary blackhole system. A very prominent example of a circumbinary discorbiting two stars is the system GG Tau. Due to the large starseparation, the whole system (stars and disc) can be directly im-aged by interferometric methods (Guilloteau et al. 1999). Laterdata yielded constraints on the size of the dust in the system andhave indicated that the central binary may consist of multiplestellar systems (see Andrews et al. (2014), Di Folco et al. (2014)and references therein). More recent observations have pointedto possible planet formation (Dutrey et al. 2014) and streamersfrom the disc onto the stars (Yang et al. 2017). A summary of theproperties of the GG Tau system is given in Dutrey et al. (2016).An additional important clue for the existence and impor-tance of circumbinary discs is given by the observed circumbi-nary planets. These are planetary systems where the planets donot orbit one star, but instead orbit a binary star system. Thesesystems have been detected in recent years by the
Kepler SpaceMission and have inspired an intense research activity. The firstsuch circumbinary planet, Kepler-16b, orbits its host systemconsisting of a K-type main-sequence star and an M-type reddwarf, in what is known as a P-type orbit with a semi-major axis of 0 . . . The observations have shownthat in these systems the stars are mutually eclipsing each other,as well as the planets. The simultaneous eclipses of several ob-jects is a clear indication of the flatness of the system becausethe orbital plane of the binary star has to coincide nearly exactlywith the orbital plane of the planet. In this respect these systemsare flatter than our own solar system. Because planets form indiscs, the original protoplanetary disc in these systems had toorbit around both stars, hence it must have been a circumbinarydisc that was also coplanar with respect to the orbital plane ofthe central binary.In addition to their occurrence around younger stars, cir-cumbinary discs are believed to orbit supermassive black holebinaries in the centres of galaxies (Begelman et al. 1980), wherethey may play an important role in the evolution of the host Additional known main-sequence binaries with circumbinary planetsbesides Kepler-16 are Kepler-34 and -35 (Welsh et al. 2012), Kepler-38 (Orosz et al. 2012a), Kepler-47 (Orosz et al. 2012b), Kepler-64(Schwamb et al. 2013), Kepler-413 (Kostov et al. 2014), Kepler-453(Welsh et al. 2015), and Kepler-1647 (Kostov et al. 2016)Article number, page 1 of 20 a r X i v : . [ a s t r o - ph . E P ] J un & A proofs: manuscript no. cb galaxies (Armitage & Natarajan 2002). In the early phase ofthe Universe there were many more close encounters betweenthe young galaxies that alread had black holes in their centres.The gravitational tidal forces between them often resulted innew merged objects that consisted of a central black hole binarysurrounded by a large circumbinary disc (see e.g. Cuadra et al.2009; Lodato et al. 2009, and references therein). The dynamicalbehaviour of such a system is very similar to the one describedabove, where we have two stars instead of black holes.To understand the dynamics of discs around young binarystars such as GG Tau or the planet formation process in cir-cumbinary discs or the dynamical evolution of a central blackhole binary in a galaxy, it is important to understand in detail thebehaviour of a disc orbiting a central binary. The strong gravi-tational perturbation by the binary generates spiral waves in thedisc; these waves transport energy and angular momentum fromthe binary to the disc (Pringle 1991) and subseqently alter theevolution of the binary (Artymowicz et al. 1991). The most im-portant impact of the angular momentum transfer to the disc isthe formation of an inner cavity in the disc with a very low den-sity whose size depends on disc parameters (such as viscosity orscale height) and on the binary properties (mass ratio, eccentric-ity) (Artymowicz & Lubow 1994).The numerical simulation of these circumbinary discs is nottrivial and over the years several attempts have been made. Thefirst simulations of circumbinary discs used the smoothed parti-cle hydrodynamics (SPH) method, a Lagrangian method wherethe fluid is modelled as individual particles (Artymowicz &Lubow 1996). On the one hand, this method has the advantagethat the whole system can be included in the computational do-main and hence the accretion onto the single binary componentscan be studied. On the other hand, because of the finite massof the SPH-particles, the maximum density contrast that can beresolved is limited.Later Rozyczka & Laughlin (1997) simulated circumbinarydiscs with the help of a finite di ff erence method on a polar grid.The grid allowed for a much larger range in mass resolution thanthe SPH codes used to date. A polar grid is well suited for theouter parts of a disc, but su ff ers from the fact that there will al-ways be an inner hole in the computational domain since theminimum radius R min cannot be made arbitrarily small. Decreas-ing the minimum radius will also strongly reduce the requirednumerical time step, rendering long simulations unfeasible. Butsince circumbinary disc simulations have to be simulated for tensof thousands of binary orbits to reach a quasi-steady state, theminimum radius is often chosen such that the motion of the cen-tral stars is not covered by the computational domain, and themass transfer from the disc to the binary cannot be studied.Günther & Kley (2002) used a special dual-grid techniqueto use the best of the two worlds, the high resolution of thegrid codes and the coverage of the whole domain of the particlecodes. For the outer disc they used a polar grid and they overlaidthe central hole with a Cartesian grid. In this way it was possibleto study the complex interaction between the binary and the discin greater detail. All these simulations indicated that despite theformation of an inner gap, material can enter the inner regionsaround the stars crossing the gap through stream like featuresto eventually become accreted onto the stars and influence theirevolution (Bate et al. 2002).Another interesting feature of circumbinary discs concernstheir eccentricity. As shown in hydrodynamical simulations forplanetary mass companions, discs develop a global eccentricmode even for circular binaries (Papaloizou et al. 2001). Theeccentricity is confined to the inner disc region and shows a very slow precession (Kley & Dirksen 2006). These results were con-firmed for equal mass binaries on circular orbits (MacFadyen &Milosavljevi´c 2008). The disc eccentricity is excited by the 1:3outer eccentric Lindblad resonance (Lubow 1991; Papaloizouet al. 2001), which is operative for a su ffi ciently cleared out gap.For circular binaries a transition, in the disc structure from morecircular to eccentric is thought to occur for a mass ratio of sec-ondary to primary above 1 /
25 (D’Orazio et al. 2016). However,this appears to contradict the mentioned results by Papaloizouet al. (2001) and Kley & Dirksen (2006) who show that this tran-sition already occurs for a secondary in the planet mass regimewith details depending on disc parameters such as pressure andviscosity.Driven by the observation of planets orbiting eccentric bi-nary stars and the possibility of studying eccentric binary blackholes, a large number of numerical simulations have been per-formed over the last few years dealing with discs around ec-centric central binaries (e.g. Pierens & Nelson 2013; Kley &Haghighipour 2014; Dunhill et al. 2015). The recent simula-tions of circumbinary discs by Pierens & Nelson (2013), Kley& Haghighipour (2014) and Lines et al. (2015) used grid-basednumerical methods on a polar grid with an inner hole. This raisesquestions about the location and imposed boundary conditionsat the innermost grid radius R min . In particular, the value of R min must be small enough to capture the development of thedisc eccentricity through gravitational interaction between thebinary and the disc, especially through the 1:3 Lindblad reso-nance (Pierens & Nelson 2013). Therefore R min has to be chosenin a way that all important resonances lie inside the computa-tional domain.In addition to the position of the inner boundary, there is nocommon agreement on which numerical boundary condition bet-ter describes the disc. In simulations concerning circumbinaryplanets mainly two types of boundary conditions are used: closedinner boundaries (Pierens & Nelson 2013; Lines et al. 2015) thatdo not allow for mass flow onto the central binary, and outflowboundaries (Kley & Haghighipour 2014, 2015) that do allow foraccretion onto the binary. While Lines et al. (2015) do not find asignificant di ff erence between these two cases, Kley & Haghigh-ipour (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 forclosed boundaries. Driven by these discussions in the literaturewe decided to perform dedicated numerical studies to test theimpact of the location of R min , the chosen boundary condition,and other numerical aspects in more detail. During our workwe became aware of other simulations that were tackling similarproblems. In October 2016 alone three publications appeared on astro-ph that described numerical simulations of circumbinarydiscs (Fleming & Quinn 2017; Miranda et al. 2017; Mutter et al.2017). While the first paper considered SPH-simulations, in thelatter two papers grid-codes were used and in both the boundarycondition at R min is discussed.In our paper we start out with numerical considerations andinvestigate the necessary conditions at the inner boundary in de-tail, and we discuss aspects of the other recent results in the pre-sentation of our findings below. Additionally, we present a de-tailed parameter study of the dynamical behaviour of circumbi-nary discs as a function of binary and disc properties.We organised this paper in the following way. In Sect. 2 wedescribe the numerical and physical setup of our circumbinarysimulations. In Sect. 3 we briefly describe the numerical methodsof the di ff erent codes we used for our simulations. In Sect. 4 weexamine the disc structure by varying the inner boundary condi- Article number, page 2 of 20aniel Thun et al.: Circumbinary discs: Numerical and physical behaviour
Table 1.
Orbital parameters of the Kepler-16 system. M A [ M (cid:12) ] M B [ M (cid:12) ] q bin a bin [au] e bin T bin [d]0.69 0.20 0.29 0.22 0.16 41.08 Notes.
Values taken from Doyle et al. (2011) (rounded to two decimals).The mass ratio is defined as q bin = M B / M A . tion and its location. The influence of di ff erent binary parametersare examined in Sect. 5. Di ff erent disc parameters and their in-fluence are studied in Sect. 6. Finally we summarise and discussour results in Sect. 7.
2. Model setup
To study the evolution of circumbinary discs we perform locallyisothermal hydrodynamical simulations. As a reference systemwe consider a binary star having the properties of Kepler-16,whose dynamical parameters are presented in Table 1. This sys-tem has a typical mass ratio and a moderate eccentricity of theorbit, and it has been studied frequently in the literature. In ourparameter study we start from this reference system and vary thedi ff erent aspects systematically. Inspired by the flatness of the observed circumbinary planetarysystems, for example in the Kepler-16 system the motion takesplace in single plane to within 0 . ◦ (Doyle et al. 2011), we makethe following two assumptions:1. The vertical thickness H of the disc is small compared to thedistance from the centre R ;2. There is no vertical motion.It is therefore acceptable to reduce the problem to two dimen-sions by vertically averaging the hydrodynamical equations. Inthis 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 theorigin of our coordinate system with the vertical axis alignedwith the rotation axis of the binary. The averaged hydrodynami-cal equations are then given by ∂∂ t Σ + ∇ · ( Σ u ) = , (1) ∂∂ t ( Σ u ) + ∇ · ( Σ u ⊗ u − Π ) = − ∇ P − Σ ∇ Φ , (2)where Σ = (cid:82) ∞−∞ (cid:37) d z is the surface density, P = (cid:82) ∞−∞ p d z the verti-cally integrated pressure, and u = ( u R , u ϕ ) T the two-dimensionalvelocity vector. We close this system of equations with a locallyisothermal equation of state P = c ( R ) Σ , (3)with the local sound speed c s ( R ).The gravitational potential Φ of both stars is given by Φ = − (cid:88) k = GM k (cid:104) ( R − R k ) + ( ε H ) (cid:105) / . (4) With r we denote the three-dimensional positional vector r = R ˆe R + z ˆe z and with R the two-dimensional positional vector in the R − ϕ -plane R = R ˆe R Here M k denotes the mass of the k -th star, G is the gravitationalconstant, and R − R k a vector from a point in the disc to theprimary or secondary star. From a numerical point of view, thesmoothing factor ε H is not necessary if the orbit of the binaryis not inside the computational domain. We use this smoothingfactor to account for the vertically extended three-dimensionaldisc in our two-dimensional case (Müller et al. 2012). For allsimulations we use a value of ε = . Π is given in tensor notation by Π i j = η (cid:34) (cid:16) u j ; i + u i ; j (cid:17) − δ i j ( ∇ · u ) (cid:35) , (5)where η = Σ ν is the vertically integrated dynamical viscositycoe ffi cient. To model the viscosity in the disc we use the α -discmodel by Shakura & Sunyaev (1973), where the kinematic vis-cosity is given by ν = α c s H . The parameter α is less than one,and for our reference model we use α = . h = H / R we assume a verticalhydrostatic equilibrium1 (cid:37) ∂ p ∂ z = − ∂ Φ ∂ z , (6)with the density (cid:37) and pressure p . To solve this equation we usethe full three-dimensional potential of the binary Φ = − (cid:80) k GM k | r − r k | and the isothermal equation of state p = c ( R ) (cid:37) (we also assumean isothermal disc in the z -direction). Integration over z yields (cid:37) = (cid:37) exp (cid:40) − z H (cid:41) , (7)with the disc scale height (Günther & Kley 2002) H = (cid:88) k c GM k | R − R k | − (8)and the midplane density (cid:37) = Σ √ π H , (9)which can be calculated by using the definition of the surfacedensity.If we concentrate the mass of the binary M bin = M A + M B inits centre of mass ( R k → H = c s u K R , (10)with the Keplerian velocity u K = √ GM bin / R . Test calculationswith the simpler disc height (10) showed no significant di ff er-ence compared to calculations with the more sophisticated discheight (8). Therefore, we use the simpler disc height in all ourcalculation to save some computation costs.For our locally isothermal simulations we use a tempera-ture profile of T ∝ R − , which corresponds to a disc withconstant aspect ratio. The local sound speed is then given by c s ( R ) = hu K ∝ R − / . If not stated otherwise we use a disc aspectratio of h = . Article number, page 3 of 20 & A proofs: manuscript no. cb
For the initial disc setup we follow Lines et al. (2015). The initialsurface density in all our models is given by Σ ( t = = f gap Σ ref R − α Σ , (11)with the reference surface density Σ ref = − M (cid:12) au − . The ini-tial slope is α Σ = . f gap = (cid:34) + exp (cid:32) − R − R gap ∆ R (cid:33)(cid:35) − (12)(Günther & Kley 2002), with the transition width ∆ R = . R gap and the estimated size of the gap R gap = . a bin (Artymowicz& Lubow 1994). The initial radial velocity is set to zero u R ( t = =
0, and for the initial azimuthal velocity we choose the localKeplerian velocity u ϕ ( t = = u K . For models run with P luto or F argo we start the binary at pe-riastron at time t = h d start the binary at apastron (lowersigns in equation (13) and (14)). R A = (cid:32) K a bin (1 ∓ e bin )0 (cid:33) , v A = K π T bin a bin (cid:113) ± e bin ∓ e bin , (13) R B = (cid:32) − K a bin (1 ∓ e bin )0 (cid:33) , v B = − K π T bin a bin (cid:113) ± e bin ∓ e bin , (14)with K = M B / M bin , K = M A / M bin , and the period of the binary T bin = π (cid:104) a / ( GM bin ) (cid:105) / . In our simulations we use a logarithmically increasing grid inthe R -direction and a uniform grid in the ϕ -direction. The phys-ical and numerical parameter for our reference system used inour extensive parameter studies in Sect. 5 and Sect. 6 are quotedin Table 2 below. In the radial direction the computational do-main ranges from R min = .
25 au to R max = . π , with a resolution of 762 × ff erentconfigurations as explained below.Because of the development of an inner cavity where the sur-face density drops significantly, we use a density floor Σ floor = − (in code units) to avoid numerical di ffi culties with too lowdensities. Test simulations using lower floor densities did notshow any di ff erences.At the outer radial boundary we implement a closed bound-ary condition where the azimuthal velocity is set to the localKeplerian velocity. At the inner radial boundary the standardcondition is the zero-gradient outflow condition as in Kley &Haghighipour (2014) or Miranda et al. (2017), but we also testdi ff erent possibilities such as closed boundaries or the viscousoutflow condition (Kley & Nelson 2008; Mutter et al. 2017).The open boundary is implemented in such a way that gascan leave the computational domain but cannot reenter it. This isdone 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 Keplerianvelocity at the inner boundary, due to the strong binary-disc in-teraction, we also use a zero-gradient boundary condition for theangular velocity Ω ϕ = u ϕ / R . By using the zero-gradient condi-tion for the angular velocity instead of the azimuthal velocity weensure a zero-torque boundary. In the ϕ -direction we use peri-odic boundary conditions.In all our simulations we use dimensionless units. The unitof length is R = M = M A + M B and the unit of time is t = (cid:113) R / ( GM ) so that the gravitational constant G is equal toone. The unit density is then given by Σ = M / R . Since our goal is to study the evolution of the disc under theinfluence of the binary we calculated the disc eccentricity e disc and the argument of the disc periastron (cid:36) disc ten times per binaryorbit. To calculate these quantities we treat each cell as a particlewith the cells mass and velocity on an orbit around the centre ofmass of the binary. Thus, the eccentricity vector of a cell is givenby e cell = u × j GM bin − R | R | (15)with the specific angular momentum j = R × u of that cell. Wehave a flat system, thus the angular momentum vector only hasa z -component. The eccentricity e cell and longitude of periastron (cid:36) cell of the cell’s orbit are therefore given by e cell = | e cell | , (16) (cid:36) cell = atan2( e y , e x ) . (17)The global disc values are then calculated through a mass-weighted average of each cell’s eccentricity and longitude of pe-riastron (Kley & Dirksen 2006) e disc = (cid:34)(cid:90) R R (cid:90) π Σ e cell R d ϕ d R (cid:35) / (cid:34)(cid:90) R R (cid:90) π Σ R d ϕ d R (cid:35) , (18) (cid:36) disc = (cid:34)(cid:90) R R (cid:90) π Σ (cid:36) cell R d ϕ d R (cid:35) / (cid:34)(cid:90) R R (cid:90) π Σ R d ϕ d R (cid:35) . (19)The integrals are simply evaluated by summing over all gridcells. The lower bound is always R = R min . For the disc ec-centricity we integrate over the whole disc ( R = R max ) if notstated otherwise, whereas for the disc’s longitude of periastronit is suitable to integrate only over the inner disc ( R = . (cid:36) disc since animations clearlyshow a precession of the inner disc (see e.g. Kley & Haghigh-ipour 2015). The radial eccentricity distribution of the disc isgiven by e ring ( R ) = (cid:34)(cid:90) R + d RR (cid:90) π Σ e cell R (cid:48) d ϕ d R (cid:48) (cid:35) / (cid:34)(cid:90) R + d RR (cid:90) π Σ R (cid:48) d ϕ d R (cid:48) (cid:35) . (20)
3. Hydrodynamic codes
Since the system under analysis is, in particular near the cen-tral binary, very dynamical we decided to compare the results
Article number, page 4 of 20aniel Thun et al.: Circumbinary discs: Numerical and physical behaviour from three di ff erent hydrodynamic codes to make sure thatthe observed features are physical and not numerical artefacts.We use codes with very di ff erent numerical approaches. P luto solves the hydrodynamical equations in conservation form witha Godunov-type shock-capturing scheme, whereas R h d andF argo are second-order upwind methods on a staggered mesh.In the following sections we describe each code and its featuresbriefly. P luto We use an in-house developed GPU version of the P luto luto solves the hydrodynamicequations using the finite-volume method which evolves volumeaverages in time. To evolve the solution by one time step, threesubsteps are required. First, the cell averages are interpolated tothe cell interfaces, and then in the second step a Riemann prob-lem is solved at each interface. In the last step the averages areevolved in time using the interface fluxes.For all three substeps P luto o ff ers many di ff erent numeri-cal options. We found that the circumbinary disc model is verysensitive to the combination of these options. Third-order inter-polation and time evolution methods lead very quickly to a neg-ative density from which the code cannot recover, even thoughwe set a density floor. We therefore use a second-order recon-struction of states and a second-order Runge-Kutta scheme forthe time evolution. Another important parameter is the limiter,which is used during the reconstruction step to avoid strong os-cillations. For the most di ff usive limiter, the minmod limiter, noconvergence is reached for higher grid resolutions. For the leastdi ff usive limiter, the mc limiter, the code again produces negativedensities and aborts the calculation. Thus, we use the van Leerlimiter which, in kind of di ff usion, lies between the minmod andthe mc limiter.For the binary position we solve Kepler’s equation using theNewton–Raphson method at each Runge-Kutta substep. R h d The R h d code is a two-dimensional radiation hydrodynamicscode originally designed to study boundary layers in accretiondiscs (Kley 1989), but later extended to perform flat disc sim-ulations with embedded objects (Kley 1999). It is based on thesecond-order upwind method described in Hawley et al. (1984)and Rozyczka (1985). It uses a staggered grid with second-orderspatial derivatives and through operator-splitting the time inte-gration is semi-second order. Viscosity can be treated explicitlyor implicitly, artificial viscosity can be applied, and the Fargo al-gorithm has also been added. The motion of the binary stars isintegrated using a fourth-order Runge-Kutta algorithm. F argo We adopted the adsg version of F argo (Masset 2000; Baruteau& Masset 2008) updated by Müller & Kley (2012). This codeuses a staggered mesh finite di ff erence method to solve the hy-drodynamic equations. Conceptually, the F argo code uses thesame methods as R h d or the Z eus code, but employs the spe-cial Fargo algorithm that avoids the time step limitations that aredue to the rotating shear flow (Masset 2000). In our situation theapplication of the Fargo algorithm is not always beneficial be-cause of the larger deviations from pure Keplerian flow near the Fig. 1.
Two-dimensional plot of the inner disc of one of our locallyisothermal Kepler-16 simulations where both orbits of the binary lieinside the computational domain. The logarithm of the surface densityis colour-coded. The orbits of the primary and secondary are shown inwhite and green. The white inner region lies outside the computationaldomain; the red cross marks the centre of mass of the binary. A videoof the simulation can be found online. 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 presenttwo important numerical issues that can have a dramatic influ-ence on the outcome of the simulations: the inner boundary con-dition (open or closed) and the location of the inner radius of thedisc. We show that “unfortunate” choices can lead to incorrectresults. In this section we use a numerical setup di ff erent fromthat in Table 2. Specifically, the base model has a radial extentfrom R min = .
25 au to R max = . ×
512 gridcells for simulations with R h d and FARGO, and512 ×
580 with P luto (see Appendix A for an explanation ofwhy we use di ff erent resolutions for di ff erent codes). For sim-ulations with varying inner radii, the number of gridcells in theradial direction is adjusted to always give the same resolution inthe overlapping domain. All simulations using a polar-coordinate grid expericence thesame problem: there is a hole in the computational domain be-cause R min cannot be zero. Usually this hole exceeds the binaryorbit. Therefore, the area where gas flows from the disc onto thebinary and where circumstellar discs around the binary compo-nents form is not part of the simulation. This complex gas flowaround the binary has been shown by e.g. Günther & Kley (2002)with a special dual-grid technique that covers the whole innercavity. Their code is no longer in use and modern e ffi cient codes(F argo and P luto ) that run in parallel do not have the option ofa dual-grid.To reproduce nevertheless some results of Günther & Kley(2002), we carried out a simulation on a polar grid with an in-ner radius of R min = .
02 au so that both orbits of the primaryand secondary lie well inside the computational domain. A snap-shot of this simulation is plotted in Fig. 1. The logarithm of
Article number, page 5 of 20 & A proofs: manuscript no. cb h Σ i ϕ (cid:2) g c m − (cid:3) log , × , × , × , × , × t = 0) . . . . . . . . R [au] . . . . . . e r i n g Fig. 2.
Azimuthally averaged surface density profiles (top) and ra-dial disc eccentricity (bottom) for our Kepler-16 reference setup (us-ing R min = .
345 au) with a closed inner boundary condition at t =
16 000 T bin . Coloured solid lines show the results from simulationswith di ff erent resolutions or di ff erent grid-spacing in the R -direction(logarithmic and uniform spacing). The initial density profile is alsoshown as a dashed black line. All simulations were performed withP luto . the surface density is colour-coded, while the orbits of the pri-mary and secondary stars are shown in white and green. The redcross marks the centre of mass of the binary, which lies outsidethe computational domain. Figure 1 shows circumstellar discsaround the binary components, as well as a complex gas flowthrough the gap onto the binary. Although including the binaryorbit inside the computational would be desirable, it is at the mo-ment not feasible for long-term simulations because of the strongtime step restriction resulting from the small cell size at the in-ner boundary. This means that for simulations where the binaryorbit is not included in the computational domain, the boundarycondition at the inner radius has to allow for flow into the innercavity, at least approximately. Two boundary conditions are usu-ally used in circumbinary disc simulations: closed boundaries(Lines et al. 2015; Pierens & Nelson 2013) and outflow bound-aries (Kley & Haghighipour 2014, 2015; Miranda et al. 2017).Given the importance of the boundary condition at R min , wecarried out dedicated simulations for a model adapted from Lineset al. (2015) and used an inner radius of R min = .
345 au, but oth-erwise it was identical to our standard model. We applied closedboundaries and open ones in order to examine their influence onthe disc structure. We also varied di ff erent numerical parameters(resolution, radial grid spacing, integrator) for this simulation se-ries. t [ T bin ] . . . . . e d i s c stdunsplit rk2uni squared t [yr] Fig. 3.
Time evolution of the total disc eccentricity for simulations of thestandard model with closed inner boundary and R min = .
345 au. Dif-ferent numerical parameters such as grid spacing, time integrator, andoperator splitting have been used. The legend refers to: std standard in-tegrator (operator and directional splitting), unsplit no directional split-ting, R k uni uniform grid, squared squared grid cells. All simulations were performed with R h d . The top panel of Fig. 2 shows the azimuthally averageddensity profiles for di ff erent grid resolutions and spacings after16 000 binary orbits. All these simulations were performed withP luto and a closed inner boundary. The strong dependence ofthe density distribution and inner gap size on the numerical setupstands out. This dependence on numerics is very surprising andnot at all expected since the physical setup was identical in allthese simulations and therefore the density profiles should all besimilar and converge upon increasing resolution. Not only doesthe gap size show this strong dependence, but also the radial ec-centricity distribution (bottom panel of Fig. 2). To examine thisdependence further, which could in principle be the result of abug in our code, we reran the identical physical setup with a dif-ferent code, R h d . The results of these R h d simulations showthe same strong numerical dependence as the P luto results. InFig. 3 the total disc eccentricity time evolution for the R h d sim-ulations is plotted. Again, one would expect that di ff erent numer-ical parameters should produce approximately the same disc ec-centricity, but our simulations show a radical change in the disceccentricity if the numerical methods are slightly di ff erent.We note that using a viscous outflow condition where theradial velocity at R min is fixed to − β ν/ (2 R ) with β = ff erence between a closed and aviscous boundary.In contrast to the closed inner boundary simulations, simula-tions with a zero-gradient outflow inner boundary reach a quasi-steady state. For the open boundaries P luto simulations with dif-ferent numerical parameters produce very similar surface densityand radial disc eccentricity profiles (Fig. 4) for di ff erent resolu-tions and grid spacings. Only the values for the uniform radialgrid with a resolution of 395 ×
512 (purple line in Fig. 4) de-viate. Here the resolution in the inner computational domain isnot high enough, because a simulation with a higher resolutionuniform grid (790 × Article number, page 6 of 20aniel Thun et al.: Circumbinary discs: Numerical and physical behaviour h Σ i ϕ (cid:2) g c m − (cid:3) log , × , × , × , × , × t = 0) . . . . . . . . R [au] . . . . . . e r i n g Fig. 4.
Azimuthally averaged surface density profiles (top) and radialdisc eccentricity (bottom) for our reference setup ( R min = .
345 au) witha zero-gradient outflow inner boundary condition after 16 000 binaryorbits. Coloured solid lines show results from simulations with di ff erentresolutions or di ff erent grid-spacing in the R -direction (logarithmic anduniform spacing). The initial density profile is also shown as a dashedblack line. All simulations were performed with P luto . imately the same results as the simulations with a logarithmicradial grid spacing.We do not have a full explanation for this strong variabilityand non-convergence of the flow when using the closed innerboundary, but this result seems to imply that this problem is ill-posed. A closed inner boundary creates a closed cavity, whichapparently implies in this case that the details of the flow dependsensitively on numerical di ff usion as introduced by di ff erent spa-tial resolutions, time integrators, or codes. As a consequence,for circumbinary disc simulations a closed inner boundary is notrecommended for two reasons. Firstly, it produces the describednumerical instabilities and secondly, for physical reasons the in-ner boundary should be open because otherwise no material canflow into the inner region and be accreted by the stars. This isalso indicated in Fig. 1 which shows mass flow through the in-ner gap, which cannot be modelled with a closed inner boundary. After determining that a zero gradient outflow condition is nec-essary, we use it now in all the following simulations and inves-tigate in a second step the optimal location of R min . It should beplaced in such a way that it simultaneously insures reliable re-sults and computational e ffi ciency. For the excitation of the disc h Σ i ϕ (cid:2) g c m − (cid:3) . . . . . . . . R [au] . . . . . e r i n g R min = 0 .
12 au R min = 0 .
20 au R min = 0 .
25 au R min = 0 .
30 au R min = 0 .
345 au R min = 0 .
40 au R min = 0 .
50 au
Fig. 5.
Influence of di ff erent inner radii on the disc structure for sim-ulations with a zero-gradient outflow inner boundary condition. Top:
Azimuthally averaged surface density profiles after 16 000 binary or-bits.
Bottom:
Radial disc eccentricity distribution at the same time. eccentricity through the binary-disc interaction, non-linear modecoupling, and the 3:1 Lindblad resonance are important (Pierens& Nelson 2013). Therefore, the location of the inner boundaryis an important parameter and R min should be chosen such thatall major mean-motion resonances between the disc and the bi-nary lie inside the computational domain. To investigate the in-fluence of the inner boundary position, we set up simulationswith an inner radius from R min = .
12 au to R min = . R > R min = . R = .
46 au) lies outside the computational domain, confirm-ing the conclusion from Pierens & Nelson (2013) about the im-portance of the 3:1 Lindblad resonance for the disc eccentric-ity growth. The variation of the radial disc eccentricity for innerradii smaller than 0 . ff erent phases after 16 000 binary orbits. Contrary to the ec-centricity distribution, the azimuthally averaged surface densityshows a stronger dependence on the inner radius. As seen in thetop panel of Fig. 5, the maximum of the surface density increasesmonotonically as the location of the inner boundary decreases Article number, page 7 of 20 & A proofs: manuscript no. cb .
20 0 .
22 0 .
24 0 .
26 0 .
28 0 .
30 0 .
32 0 . R min [au] T p r ec [ T b i n ] e bin = 0 . e bin = 0 . e bin = 0 . Fig. 6.
Variation of precession period of the inner gap of the disc forvarying inner radii, R min , and di ff erent binary eccentricities, e bin . because a smaller inner radius means also that the “area” wherematerial can leave the computational domain decreases. For in-ner radii smaller than R min ≤ .
25 au the change in the maximumsurface density becomes very low upon further reduction of R min .This implies that for too large R min there is too much mass leav-ing the domain and it has to chosen small enough to capture alldynamics. Although the maximum value of the surface densityincreases, the position of the maximum does not depend on theinner radius and remains at approximately R = . e bin . Toexplore how the location of the inner radius a ff ects the disc dy-namics, we ran additional simulations with varying R min for dif-ferent e bin . We chose e bin = .
08 and e bin = .
32 in addition tothe Kepler-16 value of e bin = .
16, which lies near the bifurca-tion point, e crit . Fig. 6 shows the precession period of the innergap for varying inner radii and three di ff erent binary eccentrici-ties, where the blue curve refers to the model shown above with e bin = .
16. For a higher binary eccentricity of e bin = .
32 (greencurve in Fig. 6 and on the upper branch of the bifurcation di-agram) the disc dynamics seems to be captured well even forhigher values of R min . Only for values lower than R min = .
22 auare deviations seen. The simulation with R min = .
20 au wasmore unstable, probably because the secondary moved in and outof the computational domain on its orbit. On the lower branch ofthe bifurcation diagram the convergence with decreasing R min isslower as indicated by the case e bin = .
08 (red curve in Fig. 6).Here, an inner radius of R bin = a bin or even slightly lower may beneeded. One explanation for this behaviour is that on the lowerbranch (for low e bin ) the inner gap is smaller but neverthelessquite eccentric such that the disc is influenced more by the lo-cation of the inner boundary. The convergence of the results forsmaller inner radii is also visible in the azimuthally averaged sur-face density and radial eccentricity profiles for e bin = .
08 and e bin = .
32 displayed in Fig. 7.
Table 2.
Setup for our reference model.
Numerical parameters R min .
25 au R max . × a bin .
22 au e bin . q bin . h . α . R min ≈ a bin . This should beused in combination with a zero-gradient outflow (hereafter justoutflow) condition. Since a smaller inner boundary increases thenumber of radial cells and imposes a stricter condition on thetime step, long-term simulations with a small inner radius arevery expensive. As a compromise to make long-term simulationspossible, we have chosen in all simulations below an inner radiusof R min = .
25 au.
While investigating the numerical conditions at the inner radius,we used an outer radius of R max = . ff erent binary and disc parameters, we found thatin some cases this assumption is not true. Especially for binaryeccentricities greater than 0 .
32, reflections from the outer bound-ary interfered with the complex inner disc structure. Therefore,we increased the outer radius of all the following simulations to R max = . = a bin , a value used by Miranda et al. (2017).In order to keep the same spatial resolution as before we alsoincreased the resolution to 762 ×
582 cells. This ensures that inthe radial direction we still have 512 grid cells between 0 .
25 auand 4 . R max ≈ a bin ) with a damping zone wherethe density and radial velocity are relaxed to their initial value.A detailed description of the implementation of such a dampingzone can be found in de Val-Borro et al. (2006).Hence, for our subsequent simulations we use from now onthe numerical setup quoted in Table 2, unless otherwise stated.To study the dependence on di ff erent physical parameter we startfrom the standard values of the Kepler-16 system in Table 2 andvary individual parameter separately.
5. Variation of binary parameters
In this section we study the influence of the central binary pa-rameter specifically its orbital eccentricity and mass ratio on thedisc structure. Throughout this section we use a disc aspect ratioof h = .
05 and a viscosity α = .
01. For the inner radius and
Article number, page 8 of 20aniel Thun et al.: Circumbinary discs: Numerical and physical behaviour h Σ i ϕ (cid:2) g c m − (cid:3) e bin = 0 . R min = 0 .
20 au R min = 0 .
22 au R min = 0 .
25 au R min = 0 .
30 au R min = 0 .
35 au e bin = 0 . R [au] . . . . . e r i n g R [au] Fig. 7.
Influence of di ff erent inner radii and binary eccentricities on the disc structure for simulations with a zero-gradient outflow inner boundary condition. Top:
Azimuthally averaged surface density profiles after 16 000 binary orbits.
Bottom:
Radial disc eccentricity distribution at the sametime.
Left: e bin = . Right: e bin = . the inner boundary condition we use the results from the previ-ous section and apply the outflow condition at R min (see also Ta-ble 2). All the results in this section were obtained with P luto ,but comparison simulations using R h d show very similar be-haviour (see Fig. A.4). Before discussing in detail the impact of varying e bin and q bin , wecomment first on the general dynamical behaviour of the innerdisc. Fig. 8 shows snapshots of the inner disc after 16 000 binaryorbits for di ff erent binary eccentricities. As seen in several otherstudies (as mentioned in the introduction) we find that the innerdisc becomes eccentric and shows a coherent slow precession.In the figure we overplot ellipses (white dashed lines) that areapproximate fits to the central inner cavity, where the semi-majoraxis ( a gap ) and eccentricity ( e gap ) are indicated. To calculate theseparameters we assumed first that the focus of the gap-ellipse isthe centre of mass of the binary. We then calculated from our datathe coordinates ( R Σ max , ϕ Σ max ) of the cell with the highest surfacedensity value, which defines the direction to the gap’s apocentre.We defined the apastron of the gap R apa as the minimum radiusalong the line ( R , ϕ Σ max ) which fulfils the condition Σ ( R , ϕ Σ max ) ≥ . · Σ ( R Σ max , ϕ Σ max ) . (21)The periastron of the gap R peri is defined analogously as the min-imum radius along the line ( R , ϕ Σ max + π ) in the opposite direction which fulfils Σ ( R , ϕ Σ max + π ) ≥ . · Σ ( R Σ max , ϕ Σ max ) . (22)Using the apastron and periastron of the gap as defined above,the eccentricity and semi-major axis of the gap are given by a gap = . (cid:16) R apa + R peri (cid:17) , (23) e gap = R apa / a gap − . (24)As shown in Fig. 8 this purely geometrical construction matchesthe shape of the central cavity very well. Even though the over-all disc behaviour shows a rather smooth slow precession, thedynamical action of the central binary is visible as spiral wavesnear the gap’s edges, most clearly seen in the first and last panel.We prepared some online videos to visualise the dynamical be-haviour of the inner disc.An alternative way to characterise the disc gap dynamicallyis by using the orbital elements ( e max , a max ) of the cell with themaximum surface density. These orbital elements can be calcu-lated with equation (15) and the vis-viva equation a = (cid:32) R − u GM bin (cid:33) − . (25)The time evolution of these gap characteristics are plotted inFig. 9. The gap’s size and eccentricity show in phase oscillatorybehaviour with a larger amplitude in the eccentricity variations. Article number, page 9 of 20 & A proofs: manuscript no. cb
Fig. 8.
Structure of the inner disc for di ff erent binary eccentricities after 16 000 binary orbits. The surface density is colour-coded in cgs-units.The orbits of the primary and secondary around the centre of mass (central point) are shown in blue and red. The white dashed lines representapproximate ellipses fitted to match the extension of the inner gap (see explanation in text). Videos of these simulations can be found online. The period of the oscillation is identical to the precession periodof the disc. Clearly, the extension of the gap always lies insideof the location of maximum density a gap < a max , but for theeccentricities the ordering is not so clear. The radial variationsof the disc eccentricity as shown in Figs. 4 and 5 indicate thatthe inner regions of the disc typically have a higher eccentricity.In Fig. 9 e disc is lowest because it is weighted with the densitywhich is very low in the inner disc regions. The eccentricity forthe gap e gap stems from a geometric fit to the very inner disc re-gions and is the highest. Inside the maximum density and evenslightly beyond that radius the disc precesses coherently in thesense that the pericentres estimated at di ff erent radii are aligned.The data show further that when the disc eccentricity is highestthe disc is fully aligned with the orbit of the secondary star, i.e.the pericentre of disc and binary lie in the same direction (seealso Appendix A). For simulations in this section we fixed q bin = .
29 and varied e bin from 0 . .
64. The binary eccentricity strongly influencesthe size of the gap as well as the disc precession period of theinner disc.In Fig. 10 the azimuthally averaged surface density profilesare shown for the various e bin at 16 000 binary orbits. In allcases there is a pronounced density maximum visible which isthe strongest for the circular binary with e bin = e bin from zero to higher values the peak shiftsinward. Not only does the gap size decrease, but the maximumsurface density also drops with increasing e bin until a minimumat e bin = .
18 is reached. Increasing e bin further causes the den-sity peak to move outward again, whereas the maximum surfacedensity stays roughly constant for higher binary eccentricities. Article number, page 10 of 20aniel Thun et al.: Circumbinary discs: Numerical and physical behaviour . . . . . . e e gap e max e disc t [ T bin ] . . . . . . a [ a b i n ] a gap a max t [yr] Fig. 9.
Evolution of gap eccentricity (top) and semi-major axis (bottom)for a binary with e bin = .
30. Shown are di ff erent measurements forthe disc eccentricity. The red lines ( e gap and a gap ) show the eccentricityand semi-major axis of the geometrically fitted ellipses, see Fig. 8. Theblue lines ( e max and a max ) show the eccentricity and semi-major axisof the cell with the maximum surface density. The green curve showsthe averaged eccentricity of the inner disc. Since the gap size varieswith time we use R = R Σ max in eq. (18) to calculate the inner disceccentricity. R [ a bin ] h Σ i ϕ (cid:2) g c m − (cid:3) e bin = 0 . e bin = 0 . e bin = 0 . e bin = 0 . e bin = 0 . e bin = 0 . e bin = 0 . e bin = 0 . e bin = 0 . e bin = 0 . Fig. 10.
Azimuthally averaged density profiles for varying binary ec-centricities after 16 000 binary orbits. For clarity we only show selectedresults from the simulated parameter space. ππ π π $ d i s c e bin = 0 . , T prec = 2119 . T bin ππ π π $ d i s c e bin = 0 . , T prec = 1915 . T bin ππ π π $ d i s c e bin = 0 . , T prec = 2948 . T bin t [ T bin ] ππ π π $ d i s c e bin = 0 . , T prec = 4027 . T bin t [yr] Fig. 11.
Disc longitude of periastron for di ff erent binary eccentricitiesfor the inner disc R < . T prec , is indicatedin the legend. As mentioned above, the eccentric disc precesses slowly andin Fig. 11 the longitude of the inner disc’s pericentre is shownversus time over 16 000 binary orbits. The disc displays a slowprograde precession with inferred period of several thousand bi-nary orbits. For the measurement of the precession period wediscarded the first 6000 binary orbits since plots of the longi-tude of periastron show that during this time span the precessionperiod is not constant yet (see Fig. 11). The period is then cal-culated by averaging over at least two full periods. To be ableto average over two periods, models with a very long T prec (e.g. e bin = .
64) were simulated for more than 16 000 binary orbits.Fig. 12 summarises the results from simulations with vary-ing binary eccentricities. The top panel of Fig. 12 shows di ff er-ent measures for the gap size for varying binary eccentricities.In addition to the radial position where the azimuthally averagedsurface density reaches its maximum ( R peak ), we plot the posi-tions where the density drops to 50 and 10 percent of its max-imum value ( R . and R . ). The value R . was used by Kley& Haghighipour (2014) as a measure for the gap size, whereasMiranda et al. (2017) used R . and Mutter et al. (2017) R peak .Starting from a non-eccentric binary the curves for R peak and R . decrease with increasing binary eccentricity. For e bin ≈ .
18 thegap size reaches a minimum and then increases again for higherbinary eccentricities. In agreement with Miranda et al. (2017)we see an almost monotonic increase of R . for increasing bi-nary eccentricities. The gap size, a gap , correlates well with R . and is always about 14% smaller. Article number, page 11 of 20 & A proofs: manuscript no. cb R [ a b i n ] R peak R . R . . . . . . . . e ga p . . . . . . . . e bin T p r ec [ T b i n ] Fig. 12.
Influence of di ff erent binary eccentricities on gap size, gap ec-centricty, and precession rate. Top: Di ff erent measures of the gap sizeare shown, averaged over several thousand binary orbits. The value R peak refers to the location of maximum of the averaged surface density (seeFig. 10), while R . and R . refer to density drops of 50 and 10 per-cent of the maximum value. Middle : Eccentricity of the gap.
Bottom:
Precession period of the disc gap.
The middle panel of Fig. 12 shows the eccentricity of the in-ner cavity, calculated with the method described in Sec. 5.1. Thedisc eccentricity also changes systematically with e bin . For cir-cular binaries it reaches e gap ≈ .
44, and then it drops down toabout 0.25 for the turning point e bin = .
18, and increases againreaching e gap ≈ . e bin = .
64. Hence, for nearlycircular binaries e gap can be higher than for more eccentric bi-naries. A similar variation of the disc’s eccentricity and preces-sion rate with binary eccentricity has been noticed by Fleming . . . . . . . . . R . [ a bin ] T p r ec [ T b i n ] . . . . .
18 0 .
20 0 .
24 0 .
28 0 . .
32 0 .
40 0 .
48 0 .
56 0 . Fig. 13.
Precession rate of the gap plotted against the radius where theazimuthally averaged surface density drops to 50 percent of its maxi-mum value ( R . ), averaged over several thousand binary orbits. Di ff er-ent dots correspond to di ff erent binary eccentricities. The dashed linesare not fit curves but are just drawn to guide the eye. & Quinn (2017) who attribute this to the possible excitation ofhigher order resonances for higher disc eccentricities.The bottom panel of Fig. 12 shows the precession period ofthe inner disc ( R < ff erent binary eccentricities. Thecurve for the precession period shows a similar behaviour to theupper curves for the gap size. To investigate further the correla-tion of the gap size (here represented by R . ) and the precessionperiod, we show in Fig. 13 the two quantities plotted againsteach other. Two branches can be seen as indicated by the dashedcurves. One starts at e bin = . e bin = .
18 wherethe gap size and precession period decrease with increasing bi-nary eccentricities, and the other branch starts at the minimumat e bin = .
18 where both gap properties increase again with in-creasing binary eccentricities. These two branches may indicatetwo di ff erent physical processes that are responsible for the cre-ation of the eccentric inner gap and show the direct correlationbetween precession period and gap size.In our simulations all discs became eccentric and showeda prograde precession, and we did not find any indication for“stand still” discs that are eccentric but without any preces-sion. In contrast, Miranda et al. (2017) find for their disc setup( h = . , α = .
1) that the disc does not precess for binary ec-centricities between 0 . .
4. We do not find this behaviourfor our disc models, but note that we use a disc with a lower as-pect ratio and a lower α -value. To investigate this further we setup a simulation with the same numerical parameters as Mirandaet al. (2017) and used a physical setup where they did not ob-serve a precession of the disc ( q bin = . e bin = . h = . α = . R min = (1 + e bin ) a bin and one with an inner radius of R min = . a bin . The first radius was used by Miranda et al.(2017) and in this case we also observe no precession of thedisc. The second radius corresponds to R min = .
25 au, whichis the radius we established in Sect. 4.2 as the optimum locationof the inner boundary. In the case with the smaller inner radiuswe observe a clear precession of the inner disc (Fig. 14) with arelatively short precession period T prec as expected for high vis-cosity and high h (see below). This is another indication of the Article number, page 12 of 20aniel Thun et al.: Circumbinary discs: Numerical and physical behaviour t [ T bin ] . ππ π π $ d i s c R min = 1 . a bin R min = 1 . a bin Fig. 14.
Disc longitude of periastron for di ff erent inner radii of the com-putational domain. The red line shows results from a simulation with theMiranda et al. (2017) setup where no precession of the disc is observ-able. The blue line shows results from the same physical setup but witha smaller inner disc radius. In this case a clear precession of the innerdisc is visible. R [ a bin ] h Σ i ϕ (cid:2) g c m − (cid:3) q bin = 0 . q bin = 0 . q bin = 0 . q bin = 0 . q bin = 0 . q bin = 0 . q bin = 0 . q bin = 0 . q bin = 0 . q bin = 1 . Fig. 15.
Azimuthally averaged density profiles for varying binary massratio after 16 000 binary orbits. necessity of choosing R min su ffi ciently small to capture all e ff ectsproperly. In this section we study the influence of the binary mass ratio.We carried out simulations with a mass ratio from q bin = . q bin = .
0; all other parameters were set according to Table 2.Fig. 15 shows the azimuthally averaged density profiles for vary-ing binary mass ratios. The density profiles do not show such astrong variation as in simulations with varying binary eccentric-ities. Only the density profiles of the two lowest simulated massratios q bin = . q bin = . ff er significantly from theother profiles. For these models the position of the peak surfacedensity is closer to the binary and the maximum surface densityroughly 20 percent lower. The di ff erent gap size measures, in-troduced in the previous section, are shown in the top panel ofFig. 16. In general, the variations of all three gap-size indicators R [ a b i n ] R peak R . R . . . . . . . . . . . q bin T p r ec [ T b i n ] Fig. 16.
Influence of di ff erent binary mass ratios on gap size and preces-sion rate. Top: Di ff erent measures of the gap size averaged over severalthousand binary orbits. Bottom:
Precession period of the disc gap. with mass ratio q bin are relatively weak. All three gap size mea-surements remain nearly constant for mass ratios greater than0 .
3. Since the density profiles do not show a significant varia-tion this is not surprising. At the same time the size and eccen-tricity of the gap do not vary significantly and lie in the range a gap ≈ a bin and e gap ≈ .
27. Overall, compared to simula-tions with varying binary eccentricity the binary mass ratio doesnot have such a strong influence on the inner disc structure for q bin (cid:38) . q , aswe 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 discstructure. We use the Kepler-16 values for the binary and ourstandard numerical setup (Table 2).While performing the simulations for this paper we observedthat models with low pressure (low H / R ) and high viscosity(high α ) seem to be very challenging for the Riemann-CodeP luto . As a consequence we also present in this section results Article number, page 13 of 20 & A proofs: manuscript no. cb . . . . R [ a b i n ] PlutoPluto (lower state)RH2D .
03 0 .
04 0 .
05 0 .
06 0 .
07 0 .
08 0 .
09 0 . h T p r ec [ T b i n ] Fig. 17.
Influence of di ff erent disc aspect ratios on gap size and pre-cession rate. Top:
Radius where the density drops to 50 percent of itsmaximum value averaged over several thousand binary orbits. Since re-sults from two di ff erent codes are shown, we omit the two other mea-surements for the gap size presented in previous sections. Bottom:
Pre-cession period of the disc gap. For e bin ≈ .
16 we observed that P luto simulations can switch between two states, an upper state (red) and alower state (green). For RH2D simulations we only observed the lowerstate (blue). This is discussed in more detail in Appendix A. using the Upwind-Code R h d . As discussed in Appendix A, thetwo codes produce results in very good agreement with eachother for our standard model. For other parameters they do notnecessarily produce results with such good agreement, but thesimulation results from both codes always show the same trend.Therefore, we decided to mix simulation results from P luto andR h d to cover a broader parameter range. First we varied the disc aspect ratio h from 0 .
03 to 0 .
1. Our simu-lation results show a decreasing gap size for higher aspect ratios(Fig. 17). The gaps precession period is again directly correlatedto the gap size, an observation we have already seen for di ff er-ent binary eccentricities. In general, a drop in T prec with higher h is expected because an increase in pressure will tend to reducethe gap size, which will in turn lead to a faster precession. Thistrend is indeed observed in our simulations for higher h . For sim-ulations with h < .
05 we observed a decrease in gap size andprecession period. The drop in T prec with lower values of h maybe partly due to the lack of numerical resolution and partly tothe lack of pressure support, which may no longer allows forcoherent disc precession. R [ a b i n ] R peak R . R . .
02 0 .
04 0 .
06 0 .
08 0 . α T p r ec [ T b i n ] Fig. 18.
Influence of di ff erent disc viscosities on gap size and preces-sion rate. Top: Di ff erent measures of the gap size averaged over severalthousand binary orbits. Bottom:
Precession period of the disc gap.
In this section we discuss the influence of the magnitude ofviscosity on the disc structure and therefore set up simulationswith di ff erent viscosity coe ffi cients ranging from α = .
001 to α = .
1. We then analysed the structure and behaviour of theinner cavity as before. The results are summarised in Fig. 18. Aclear trend is visible for the gap size and for the precession periodof the gap. For higher viscosities the gap size shrinks and the pre-cession period decreases. Again, we see the direct correlation ofgap size and precession rate, as already observed for aspect ratioand binary eccentricity variations. The decrease in the disc’s gapsize can be explained: for higher α values the viscous spread-ing of the disc increases. This viscous spreading counteracts thegravitational torques of the binary, which are responsible for thegap creation.
7. Summary and discussion
Using two-dimensional hydrodynamical simulations, we studiedthe structure of circumbinary discs for di ff erent numerical andphysical parameters. Since these simulations are, from a numer-ical point of view, not trivial and often it is not easy to distin-guish physical features and numerical artefacts, we checked ourresults with three di ff erent numerical codes. In the following wesummarise the most important results from our simulations withdi ff erent numerical and physical parameters.Concerning the numerical treatment the most crucial issuewith simulations using grid codes in cylindrical coordinates is Article number, page 14 of 20aniel Thun et al.: Circumbinary discs: Numerical and physical behaviour the treatment of the inner boundary condition. As there has beensome discussion in the literature about this issue (Pierens & Nel-son 2013; Kley & Haghighipour 2014; Lines et al. 2015; Mi-randa et al. 2017; Mutter et al. 2017), we decided to perform acareful study of the location of the inner boundary and the con-ditions on the hydrodynamical variables at R min . Our results canbe summarised as follows: – Inner boundary condition . Since there was no agreementon which type of inner boundary condition should be usedfor circumbinary disc simulations, especially if closed oropen boundary conditions should be used, we investigatedthe influence of the inner boundary condition through ded-icated simulations with two di ff erent codes. We observedthat closed boundaries can lead to numerical instabilities forall codes used in our comparison and no convergence to aunique solution could be found. The use of a viscous out-flow condition where the velocity is set to the mean viscousdisc speed produced very similar results to the closed bound-ary because the viscous speed is very low in comparison tothe velocities induced by the dynamical action of the binarystars. Open inner boundaries, on the other hand, lead to nu-merically stable results and are also, from a physical point ofview, more logical because material can flow onto the binaryand be accreted by one of the stars. Hence, open boundarieswith free outflow are the preferred conditions at R min . – Location of inner boundary . The location of the innerboundary is also an important numerical parameter. Since aninner hole in the computational domain cannot be avoided inpolar coordinates, the inner radius has to be su ffi ciently smallto capture all relevant physical e ff ects. In particular, all im-portant mean motion resonances, which may be responsiblefor the development of the eccentric inner cavity, should lieinside the computational domain. Through a parameter studywe were able to determine that the radius of the inner bound-ary R min has to be of the order of the binary separation a bin to capture all physical e ff ects in a proper way.After having determined the best values for the numerical issueswe performed, in the second part of the paper a careful study todetermine the physical aspects that determine the dynamics ofcircumbinary discs. Starting from a reference model, where wechose the binary parameter of Kepler-16, we varied individualparameters of the binary and the disc and studied their impact onthe disc dynamics. First of all, for all choices of our parametersthe models produce an eccentric inner disc that shows a coherentprograde precession. However, the size of the gap, its eccentric-ity, and its precession rate depend on the physical parameters ofbinary and disc. Our results can be summarised as follows: – Binary parameters To study the influence of the binary staron the disc we systematically varied the eccentricity ( e bin )and the mass ratio ( q bin ) of the binary. The parameter studyshowed that e bin has a strong influence on the gap size and onthe precession period of the gap. We found that two regimesexist where the disc behaves di ff erently to an increasing bi-nary eccentricity (see Fig. 8). The two branches bifurcate ata critical binary eccentricity, e crit = .
18 from each other.From e bin = . e crit the gap size and precession perioddecrease, and from e bin = .
18 onward both gap parametersbecome larger again, as displayed in Fig. 13. The bifurcationof the two branches near e crit = .
18 strongly suggests thatdi ff erent physical mechanisms, responsible for the creationof the eccentric inner cavity, operate in the two regimes. For the lower branch the excitation at the 3:1 outer Lindblad res-onance may be responsible, which is supported by a simu-lation where the inner boundary was outside the 3:1 radiusand no disc eccentricity was found. For the upper branch(high e bin ) non-linear e ff ects may be present, as suggestedby Fleming & Quinn (2017) who found similar behaviour.As second binary parameter we varied the mass ratio of thesecondary to the primary star, q bin , between [0 . , .
0] andstudied its impact on the gap size and precession period.Overall the variation of R . with q bin is weak. For low q bin the gap size increases until it becomes nearly constant for q bin (cid:38) . R . ≈ . a bin , as shown in Fig. 16. Theprecession period decreases on average with increasing T prec and for large q bin the behaviour is equivalent to a single par-ticle at a separation of 4 . – Disc parameters We also varied di ff erent disc parameters,namely the pressure (through the aspect ratio h = H / R )and the viscosity (through α ). The results are displayed inFigs. 17 and Fig. 18. For changes in the aspect ratio no cleartrend was visible; Both gap size and precession period firstincrease with increasing h until they reach a maximum at h = .
05; from there both gap properties decrease again withincreasing aspect ratio. For high h the behaviour can be un-derstood in term of the gap closing tendency of higher discpressure. On the other hand, for very low pressure it maybe more di ffi cult to sustain a coherent precession of a largeinner hole in the disc due to the reduced sound speed. Addi-tionally, the damping action of the viscosity becomes moreimportant for discs with lower sound speed.A clear monotonic trend was visible for the viscosity vari-ation, where higher α leads to smaller gaps and shorter gapprecession periods. This can be directly attributed to the gapclosing tendency of viscosity.To allow for an alternative view of the impact of di ff erent bi-nary and disc parameters on the dynamical structure of the discwe display in Fig. 19 the eccentricity (top) and the precessionperiod (bottom) of the gap plotted against the semi-major axisof the gap for all our simulations where di ff erent colours standfor di ff erent model series. P luto simulations are represented bydots whereas squares stand for R h d simulations. The refer-ence model is marked with a star. Interestingly the majority ofthe models lie on a main-sequence branch, which has the trendof increasing eccentricity and precession period with increasinggap width. The smaller branch that bifurcates near the referencemodel represents the low e bin sequence, which has higher e gap but smaller T prec . By coincidence, our reference model with theKepler-16 parameters lies near the bifurcation point.How the location of the bifurcation point depends of the discand binary parameter needs to be explored in subsequent numer-ical studies. A first simulation series with q bin = . e bin the models can have the tendency to switchbetween the two states during a single simulation. The physi-cal cause of the bifurcation needs to be analysed in subsequentsimulations.In addition, there are a more possibilities to further improveour model. First, the use of a non-isothermal equation should beconsidered to study the e ff ects of viscous heating and radiativecooling. Especially for simulations, which also evolve planetsinside the disc, Kley & Haghighipour (2014) showed that radia-tive e ff ects should been taken into account. Self-gravity e ff ects Article number, page 15 of 20 & A proofs: manuscript no. cb . . . . . . . . . e ga p Variation of e bin Variation of q bin Variation of h Variation of α a gap [ a bin ] T p r ec [ T b i n ] Fig. 19.
Eccentricity (top) and precession period (bottom) of the innergap plotted against the semi-major axis of the gap. Di ff erent simula-tion series are colour-coded. Both the measured gap eccentricity andgap semi-major axis are averaged over several thousand binary orbits.P luto simulations are represented by dots and squares stand for R h d simulations. The reference model with Kepler-16 parameters is markedwith the star. have been studied by Lines et al. (2015) and Mutter et al. (2017)and are important for galactic discs around a central binary blackhole. The backreaction of the disc onto the binary could also beconsidered in more detail. Angular momentum exchange withthe binary occurs through gravitational torques from the disc anddirect accretion of angular momentum from the material that en-ters the central cavity. The first part always leads a decrease inthe binary semi-major axis and an increase in eccentricity (Kley& Haghighipour 2015), while Miranda et al. (2017) pointed outthat angular momentum advection might even be dominant lead-ing to a binary separation.Another interesting question is the evolution of planets insidethe circumbinary disc and their influence on the disc structure.Since an in situ formation of the observed planets seems veryunlikely, the migration and final parking position of planets isan interesting topic that requires further investigations. In viewof our extensive parameter studies there is the indication that inseveral past studies (Pierens & Nelson 2013; Kley & Haghigh-ipour 2014, 2015) an inner radius was considered that was toolarge and that did not allow for full disc dynamics.Through fully three-dimensional simulations it will be pos-sible to study the dynamical evolution of inclined discs. Acknowledgements.
Daniel Thun was funded by grant KL 650 /
26 of the GermanResearch Foundation (DFG), and Giovanni Picogna acknowledges DFG supportgrant KL 650 /
21 within the collaborative research program “The first 10 Million Years of the Solar System”. Most of the numerical simulations were performedon the bwForCLuster BinAC, supported by the state of Baden-Württembergthrough bwHPC, and the German Research Foundation (DFG) through grantINST 39 / References
Andrews, S. M., Chandler, C. J., Isella, A., et al. 2014, ApJ, 787, 148Armitage, P. J. & Natarajan, P. 2002, ApJ, 567, L9Artymowicz, P., Clarke, C. J., Lubow, S. H., & Pringle, J. E. 1991, ApJ, 370,L35Artymowicz, P. & Lubow, S. H. 1994, ApJ, 421, 651Artymowicz, P. & Lubow, S. H. 1996, ApJ, 467, L77Baruteau, C. & Masset, F. 2008, ApJ, 672, 1054Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS, 336, 705Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307Cuadra, J., Armitage, P. J., Alexander, R. D., & Begelman, M. C. 2009, MNRAS,393, 1423de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529Di Folco, E., Dutrey, A., Le Bouquin, J.-B., et al. 2014, A&A, 565, L2D’Orazio, D. J., Haiman, Z., Du ff ell, P., MacFadyen, A., & Farris, B. 2016, MN-RAS, 459, 2379Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602Dunhill, A. C., Cuadra, J., & Dougados, C. 2015, MNRAS, 448, 3545Dutrey, A., Di Folco, E., Beck, T., & Guilloteau, S. 2016, A&A Rev., 24, 5Dutrey, A., di Folco, E., Guilloteau, S., et al. 2014, Nature, 514, 600Fleming, D. P. & Quinn, T. R. 2017, MNRAS, 464, 3343Georgakarakos, N. & Eggl, S. 2015, ApJ, 802, 94Guilloteau, S., Dutrey, A., & Simon, M. 1999, A&A, 348, 570Günther, R. & Kley, W. 2002, A&A, 387, 550Hawley, J. F., Smarr, L. L., & Wilson, J. R. 1984, ApJS, 55, 211Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90Kley, W. 1989, A&A, 208, 98Kley, W. 1999, MNRAS, 303, 696Kley, W. & Dirksen, G. 2006, A&A, 447, 369Kley, W. & Haghighipour, N. 2014, A&A, 564, A72Kley, W. & Haghighipour, N. 2015, A&A, 581, A20Kley, W. & Nelson, R. P. 2008, A&A, 486, 617Kostov, V. B., McCullough, P. R., Carter, J. A., et al. 2014, ApJ, 784, 14Kostov, V. B., Orosz, J. A., Welsh, W. F., et al. 2016, ApJ, 827, 86Leung, G. C. K. & Lee, M. H. 2013, ApJ, 763, 107Lines, S., Leinhardt, Z. M., Baruteau, C., Paardekooper, S.-J., & Carter, P. J.2015, A&A, 582, A5Liu, B., Muñoz, D. J., & Lai, D. 2015, MNRAS, 447, 747Lodato, G., Nayakshin, S., King, A. R., & Pringle, J. E. 2009, MNRAS, 398,1392Lubow, S. H. 1991, ApJ, 381, 259MacFadyen, A. I. & Milosavljevi´c, M. 2008, ApJ, 672, 83Masset, F. 2000, A&AS, 141, 165Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228Miranda, R., Muñoz, D. J., & Lai, D. 2017, MNRAS, 466, 1170Moriwaki, K. & Nakagawa, Y. 2004, ApJ, 609, 1065Müller, T. W. A. & Kley, W. 2012, A&A, 539, A18Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123Mutter, M. M., Pierens, A., & Nelson, R. P. 2017, MNRAS, 465, 4735Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012a, ApJ, 758, 87Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012b, Science, 337, 1511Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263Pierens, A. & Nelson, R. P. 2013, A&A, 556, A134Pringle, J. E. 1991, MNRAS, 248, 754Rozyczka, M. 1985, A&A, 143, 59Rozyczka, M. & Laughlin, G. 1997, in Astronomical Society of the Pacific Con-ference Series, Vol. 121, IAU Colloq. 163: Accretion Phenomena and RelatedOutflows, ed. D. T. Wickramasinghe, G. V. Bicknell, & L. Ferrario, 792Schwamb, M. E., Orosz, J. A., Carter, J. A., et al. 2013, ApJ, 768, 127Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337Welsh, W. F., Orosz, J. A., Carter, J. A., & Fabrycky, D. C. 2014, in IAU Sym-posium, Vol. 293, IAU Symposium, ed. N. Haghighipour, 125–132Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475Welsh, W. F., Orosz, J. A., Short, D. R., et al. 2015, ApJ, 809, 26Yang, Y., Hashimoto, J., Hayashi, S. S., et al. 2017, AJ, 153, 7 Article number, page 16 of 20aniel Thun et al.: Circumbinary discs: Numerical and physical behaviour
Table A.1.
Disc precession rate for di ff erent resolutions. Resolution T prec [ T bin ]448 ×
512 1277 . ×
580 1321 . × . 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 chal-lenge for numerical schemes of grid codes. We therefore study inthis appendix the dependence of the disc structure on numericalparameters, such as the numerical scheme or the grid resolution.First, we simulated our reference system, using an open innerboundary and an inner radius of R min = .
25 au, with the di ff er-ent codes described in Sec. 3.In Fig. A.1 the long-time evolution of the disc eccentric-ity (top) and of the longitude of periastron (bottom) are shown.All three codes produce comparable results, although we shouldnote that to get this agreement simulations performed with P luto need a slightly higher resolution. P luto simulations with a res-olution of 448 ×
512 grid cells produce a slightly smaller pre-cession period of the inner disc (red line in Fig. A.3, see alsoTable A.1). One possible explanation for this could be the sizeof the numerical stencil. F argo and R h d have a very compactstencil because they use a staggered grid. In contrast, P luto hasa wider stencil because a collocated grid is used, where all vari-ables are defined at the cell centres. In particular the viscosityroutine of P luto has a very wide stencil to achieve the correctcentring of the viscous stress tensor components. This could ex-plain why for P luto a slightly higher resolution is needed.The π -shift of the longitude of periastron calculated withR h d compared to the values calculated with P luto and F argo (bottom Panel of Fig. A.1) can be explained in the followingway. Simulations performed with P luto and F argo started thebinary at periastron with the pericentre of the secondary at ϕ = π (see red ellipses in Fig. 8), whereas R h d started the binary atapastron with the pericentre of the secondaries orbit at ϕ = h d simulations is shifted by π we can also see this π -shift in the disc’s longitude of periastron.Fig. A.2 shows the azimuthally averaged density profiles at t =
16 000 T bin . Even after such a long time span all three codesagree very closely. All codes produce roughly the same densitymaximum and the same density slopes.To check the influence of the numerical resolution on thedisc structure we run P luto simulations with di ff erent resolu-tions, summarised in Table A.1. The results of this resolutiontest are shown in Fig. A.3. The disc eccentricity converges to thesame value, whereas it seems that the disc precession rates in-creases with higher resolution. A closer look at the data showsthat after around 6000 orbits the precession rate for the higherresolution converged. Table A.1 summarises the measured pre-cession rate for di ff erent resolutions. Since there is no large vari-ation between 512 ×
580 and 600 × ×
580 grid cells for our Pluto simulations.As already discussed in Sect. 4.2 in some cases (for exam-ple high binary eccentricities) an outer radius of R max = . R max = . = . . . . . . . . e d i s c PlutoFargoRH2D t [ T bin ] . ππ π π $ d i s c t [yr] Fig. A.1.
Comparison of disc eccentricity (top) and disc longitudeof pericentre (bottom) for the standard model. The disc eccentricitywas calculated by summing over the whole disc ( R = R max in equa-tion (18)), whereas the disc longitude of pericentre was calculated onlyfor the inner disc ( R = . π -shift of the greencurve in the bottom panel is explained in the text. We used a resolutionof 512 ×
580 grid cells for P luto simulations and for F argo and R h d simulations a resolution of 448 × a bin (and to ensure the same resolution between 0 .
25 au and4 . × ff ect of di ff erent outer radii for thecase of varying binary eccentricities. For binary eccentricitiesgreater than 0 .
32, a larger outer radius makes a di ff erence. Also,around the bifurcation point, at e bin = .
18, a larger outer ra-dius changes the gap size and precession rate of the disc. Wealso added results from R h d simulations with an outer radius of R max = . a bin . These R h d results follow the same trend asthe P luto simulations, although for most binary eccentricities wecould not reproduce the very good agreement of the e bin = . e bin = .
18, which is the critical binary eccentricitywhere the two branches bifurcate, simulations tend to switch be-tween two states. This jump happens after roughly 3000 binaryorbits and does not depend on physical parameters. Changes innumerical parameters, like the Courant-Friedrichs-Lewy condi-tion or R max , can trigger the jump. This can be seen in Fig. A.4for e bin = .
16 where simulations with R max = . a bin (redcurve) and R max = . a bin (blue curve) produce di ff erent val-ues for the gap size and the precession period. An outer radiusof R max = . a bin is in the case of e bin = .
16 not too small,since a simulation with R max = a bin again produced the samevalues for the gap size and the precession period. For simulationswith e bin far away from e crit , we did not observe these jumps be- Article number, page 17 of 20 & A proofs: manuscript no. cb . . . . . . . . R [au] h Σ i ϕ (cid:2) g c m − (cid:3) PlutoFargoRH2D
Fig. A.2.
Azimuthally averaged density profiles calculated with di ff er-ent codes at t =
16 000 T bin . . . . . . . . . e d i s c × × × t [ T bin ] . ππ π π $ d i s c t [yr] Fig. A.3.
Comparison of di ff erent resolutions. Top:
Disc eccentricity.
Bottom:
Disc longitude of pericentre. All simulations were performedwith P luto . tween states. So far we have only observed these jumps in simu-lations carried out with P luto . . . . . . . . . . R . [ a b i n ] Pluto , R max = 18 . a bin Pluto , R max = 70 . a bin RH2D , R max = 18 . a bin . . . . . . . . e bin T p r ec [ T b i n ] Fig. A.4.
Influence of di ff erent outer radii on gap size and precessionrate while varying the binary eccentricity. Top:
Radius where the az-imuthally averaged surface density drops to 50 percent of its maximumvalue.
Bottom:
Precession period of the disc gap.
Appendix B: Comparison with massless testparticle
In our circumbinary disc evolution we found that the inner discbecomes eccentric with a constant precession rate. In this sectionwe compare this to test particle trajectories around the binary andinvestigate how a massless particle with orbital elements similarto the disc gap (Fig. 19) would behave under the influence of thebinary potential. In particular, we would like to know at whatdistance from the binary a test particle has to be positioned suchthat its precession period agrees approximately with that of thedisc. Using secular perturbation theory for the coplanar motionof a massless test particle around an eccentric binary star Mori-waki & Nakagawa (2004) derived the following formula for theparticle’s precession period for low binary eccentricities T prec =
43 ( q bin + q bin (cid:32) a p a bin (cid:33) / (cid:32) + e (cid:33) − T bin , (B.1)where a p is the semi-major axis of the particle. The same relationwas 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) whenvarying of binary and planet parameters, we performed seriesof three-body simulations with very low planet mass (10 − M (cid:12) ),where we varied individually q bin , e bin , and a p . In addition wevaried the planet eccentricity e p . For the simulations we used the Article number, page 18 of 20aniel Thun et al.: Circumbinary discs: Numerical and physical behaviour parameter of the Kepler-16 systems as a reference (see Table 1)and integrated the system for several 10 000 binary orbits. Ourresults of these three-body simulations are shown in Fig. B.1. Wefind that the precession rate scales with q bin , a p , and e bin exactlyas expected from relation (B.1). The fourth panel of Fig. B.1indicates that the precession period also depends on the particleeccentricity as T prec ∼ (1 − e ) , where we plotted the averageeccentricity of the orbit which is equivalent to the particle’s freeeccentricity. The agreement holds up to about e ∼ . a p . Clearly, for higher values of e p the particle’s orbit willbe unstable as this leads to close encounters with the binary. Thevalue of e p where this happens depends on the distance fromthe binary star. For lower a p the range will be more limited. Insummary, we find for the precession period of a particle arounda binary star T prec =
43 ( q bin + q bin (cid:32) a p a bin (cid:33) / (cid:16) − e (cid:17) (cid:16) + e (cid:17) T bin . (B.2)This relation is plotted in Fig. B.1 as the solid black line. Inaddition to the scaling behaviour we find that the numerical re-sults are about 3 −
5% lower than the theoretical estimate. Theagreement of the numerical results with the theoretical predic-tion becomes better for larger a p because eq. (B.2) is derivedfrom 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 sim-ulations. The best option for this comparison is to check thescaling of the precession rate for models where q bin has beenvaried because for binary mass ratio greater than q bin = . a gap = . a bin and e gap = .
27. As seen from eq. (B.2), a testparticle with these orbit elements would have a precession pe-riod that is far too short. However, if we assume that the testparticle has a semi-major axis of a p = . a bin (roughly 20 per-cent more than a gap ), the precession period of the gap and theparticle match very well for di ff erent q bin (Fig. B.2), at least forthe higher mass ratios. In general, however, due to the strongsensitivity of the precession period with a p , an exact agreementwill be di ffi cult to obtain. . . . . . . q bin T p r ec [ T b i n ] − body ∝ ( q bin +1) q bin a p [ a bin ] T p r ec [ T b i n ] − body ∝ (cid:16) a p a bin (cid:17) / . . . . . . . . . e bin T p r ec [ T b i n ] − body ∝ (cid:0) e (cid:1) − . . . . . . e p T p r ec [ T b i n ] − body ∝ (cid:0) − e (cid:1) Fig. B.1.
Scaling of the precession period of a particle orbiting a binarystars. As base binary parameter we chose the Kepler-16 parameters andfor the planet a p = .
50 au and e p = .
05. As an exception we chose forthe 1st panel a p = . & A proofs: manuscript no. cb . . . . . . . . . . q bin T p r ec [ T b i n ] gap3 − bodytheoretical Fig. B.2.
Comparison of the precession period of the inner disc withthe precession period of free particle. For the particle we used a p = . a bin , e p = .
25. To obtain the shown agreement the particle’s semi-major axis a p has to be roughly 20 percent longer than a gapgap