Issue 
A&A
Volume 670, February 2023



Article Number  A112  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202244828  
Published online  15 February 2023 
Threedimensional evolution of radiative circumbinary discs: The size and shape of the inner cavity
^{1}
Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux, CNRS,
B18N, allée Geoffroy SaintHilaire,
33615
Pessac,
France
email: arnaud.pierens@ubordeaux.fr
^{2}
Astronomy Unit, Queen Mary University of London,
Mile End Road,
London
E1 4NS,
UK
Received:
29
August
2022
Accepted:
5
November
2022
In this study we present the results of 3D hydrodynamical simulations of circumbinary discs that orbit around analogues of the Kepler16 and Kepler34 systems, including the effect of stellar heating and radiative cooling on the thermal disc structure. We find that, compared to their 2D counterparts, the structures of the cavities in 3D circumbinary disc models appear to reach a quasistationary state more rapidly, and in a subset of our runs the evidence for this is unambiguous. Furthermore, the sizes and eccentricities of the inner cavity are smaller in 3D compared to 2D. We attribute this difference to enhanced spiral wave dissipation in disc regions above the midplane, where the cooling time is of the order of the dynamical timescale, resulting in smaller inner cavity sizes in 3D disc models. Our results suggest that migrating planets should park closer to the central binary in 3D models of circumbinary discs, and point to the importance of including the 3D structure when simulating circumbinary discs and planets.
Key words: accretion, accretion disks / methods: numerical / hydrodynamics / planetdisk interactions / planets and satellites: formation
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
To date, 15 circumbinary planets have been discovered around eclipsing binary systems that contain main sequence stars. Most of these, such as Kepler16b and Kepler34b, were detected by the Kepler spacecraft (Doyle et al. 2011; Welsh et al. 2012), and recently TESS has also detected two circumbinary planets, TOI1338b (Kostov et al. 2020) and TIC172900988 (Kostov et al. 2021). The first detection of a circumbinary planet via the radial velocity method recently confirmed that Kepler16b has a mass similar to that of Saturn (Triaud et al. 2021), in agreement with the mass estimate obtained by fitting a photodynamical model to the transit data (Doyle et al. 2011).
The majority of these circumbinary planets (10 out of 15) orbit just outside the dynamical instability zone (Holman & Wiegert 1999), where perturbations due to the binary would cause planets to be ejected. Simulations show that the structure of the circumbinary disc would be highly perturbed in this region during the planet formation epoch, and numerous studies have shown that the formation of a significantly eccentric precessing cavity makes it highly unlikely that in situ planet formation via planetesimal accretion can occur there (Marzari et al. 2008; Lines et al. 2016). Furthermore, recent work has shown that the inner eccentric regions of circumbinary discs are prone to a parametric instability that generates hydrodynamical turbulence (Papaloizou 2005; Barker & Ogilvie 2014), whose impact on planet formation is twofold (Pierens et al. 2020, 2021): (i) it makes the collision velocities between small particles too high to allow grain growth, such that forming a massive planetesimal seed for pebble accretion becomes difficult; and (ii) it stirs up pebbles and renders pebble accretion very inefficient. Hence, the formation of circumbinary planets in situ near the dynamical instability zone via pebble accretion also appears to be very challenging. Instead, circumbinary planets could form at large distances from the binary, where dust growth and pebble accretion are more efficient, and then migrate in before stalling at their currently observed locations.
There is a substantial body of work that has examined the migration of planets in circumbinary discs, and a general picture that has emerged is that a planet embedded in a circumbinary disc migrates to the inner edge of the disc cavity. For moderate binary eccentricities, migration stalls at this location due to a positive corotation torque being exerted on the planet that counteracts the Lindblad torque (Masset et al. 2006; Pierens & Nelson 2007; Kley & Haghighipour 2014; Penzlin 2021). For binaries with higher eccentricities, the planet can acquire significant eccentricity due to its interaction with the highly eccentric circumbinary disc, which makes migration stall farther away from the central binary, at distances generally beyond the observed locations (Pierens & Nelson 2013; Kley & Haghighipour 2015). This indicates that models in which the disc eccentricity remains small are generally better for providing a good fit of the observed orbital parameters. For example, Mutter et al. (2017) and Penzlin (2021) recently found that a migrating planet with mass sufficient to open a partial gap tends to circularise the inner cavity and subsequently migrates closer to the central binary. In a recent paper, Coleman et al. (2022) arrived at the same conclusion in a scenario where the disc is circularised as a result of the effect of the dust feedback onto the gas.
In this paper, we examine whether or not 3D, rather than 2D, hydrodynamical simulations may also lead to a better agreement between simulated circumbinary systems and observations. Our motivation stems from the fact that, in our previous simulations of Kepler413 analogues (Pierens & Nelson 2018), we found that the density peak associated with the outer edge of the inner cavity was located at a similar location as that of the circumbinary planet Kepler413b when simulated in 3D, whereas the peak was located farther out in 2D simulations. In comparison with our previous works, in this paper we have significantly improved the modelling of the disc thermodynamics by incorporating irradiation heating from the central binary. This is important since it has been shown that accurate treatment of the thermodynamics inside the cavity has a significant impact on modelling circumbinary discs (Sudarshan et al. 2022). Although it remains uncertain to what degree some of our simulations have converged to a steady state, because of the computational expense of running 3D simulations, our results suggest that 3D models of circumbinary discs that include irradiation and cooling reach a quasisteady state much earlier than 2D models, and that both the size and eccentricity of the inner cavity are smaller in 3D compared to 2D.
This paper is organised as follows. In Sects. 2 and 3, we describe the physical and numerical models, respectively. In Sect. 4, we present the results of our 3D and 2D simulations of circumbinary discs around binaries with parameters equivalent to Kepler16(AB) and Kepler34(AB). We discuss our results in Sect. 5 and summarise our main findings in Sect. 6.
2 The physical models
2.1 The 3D model
We solved the hydrodynamical equations for the conservation of mass, momentum, and internal energy in spherical coordinates (r,θ,ϕ: radial, polar, azimuthal), with the origin of the frame located at the centre of mass of the binary. In particular, the energy equation that is incorporated in the code is given by (1)
where v is the velocity, and e = ρc_{υ}T is the internal energy per unit volume, with ρ the density, c_{υ} the specific heat capacity at constant volume, and T the temperature. In Eq. (1), γ is the adiabatic index, which is set to γ = 1.4.
We considered a cooling scheme where the temperature is relaxed towards a reference temperature, T_{irr}, on a cooling timescale, t_{cool}. In Eq. (1), the cooling rate, Q_{cool}, is therefore given by (2)
Here, we assume T_{irr} to correspond to the temperature set by stellar irradiation. During the simulations, this temperature was regularly updated using the RADMC3D Monte Carlo radiative transfer code (Dullemond et al. 2012), applying the procedure described in Sect. 3.2. For the 3D models, the temperature structure was updated every 100 binary orbits. The cooling timescale, t_{cool}, was calculated at every timestep by considering the timescale for radiative loss of energy from a Gaussian sphere with a characteristic length scale corresponding to the gas scale height, H (Lyra et al. 2016). It is given by (3)
where σ is the StefanBoltzmann constant and τ_{eff} is the effective optical depth given by (4)
Following Bae et al. (2016), the optical depth, τ, was calculated as 1/τ = 1/τ_{upper} + 1/τ_{lower} with (5)
where the opacity, κ, is computed using the Rosseland mean opacity of Zhu et al. (2009).
2.2 The 2D model
We also performed 2D simulations for direct comparison with the 3D models described above. These additional 2D simulations allow the circumbinary disc evolution to be calculated over tens of thousands of binary orbits for a reasonable computational cost, whereas evolving a 3D model beyond ~2000 binary orbits becomes prohibitive in terms of computational effort.
In the 2D models we considered the vertically integrated versions of the equations governing the disc evolution, and worked in 2D (R, φ) polar coordinates, with the origin of the frame located at the centre of mass of the binary. The energy equation that is solved is again given by Eq. (1), but now with the internal energy defined as e = Σc_{υ}T, with Σ the surface density, and with the optical depth given by (7)
Here, c_{1} is a correction factor introduced to account for the drop in opacity with vertical height Müller & Kley (2012). Again, a description of the procedure employed to compute the irradiation temperature that appears in Eq. (2) is provided in Sect. 3.2. Compared to 3D calculations in which T_{irr} are updated every 100 binary orbits, we note that T_{irr} is updated every 1000 binary orbits in the 2D models.
3 The numerical model
3.1 Setup
The simulations presented in this paper were performed using the multifluid version of FARGO3D (BenítezLlambay & Masset 2016). Computational units were chosen such that the total mass of the binary is M_{⋆} = 1, the gravitational constant is G = 1, and the radius R = 1 in the computational domain corresponds to the binary semimajor axis, a_{bin}. When presenting the simulation results, unless otherwise stated, we used the binary orbital period as the unit of time.
In the 3D models, the computational domain in the radial direction extended from R_{in} = 1.13 a_{bin} to R_{out} = 18 a_{bin} and we employed 995 logarithmically spaced grid cells. Ideally, the inner boundary should be set to R_{in} ≈ a_{bin} to allow for accurate modelling of mass flow and angular momentum transfer around the binary (Mutter et al. 2017). As shown by Thun et al. (2017), however, setting R_{in} = 1.13 a_{bin} appears to be a good compromise between an increased run time due to a smaller timestep and realistic modelling. For R_{in} ≲ 1.13 a_{bin}, the change in surface density appears to be very small upon further reduction of R_{in}. In the azimuthal direction, the simulation domain extended from 0 to 2π with 700 uniformly spaced grid cells. In the meridional direction, the simulation domain covered 3.5 disc pressure scale heights above and below the disc midplane, and we adopted 144 uniformly spaced grid cells.
In the 2D models, the computational domain in the radial direction extended from R_{in} = 1.13 a_{bin} to R_{out} = 40 a_{bin}, and we employed 1280 logarithmically spaced grid cells, such that the radial resolution was equivalent to that used in the 3D models. Compared to 3D models, the location of the outer boundary in the 2D simulations was set to a larger value due to a clear tendency for getting larger inner cavities when using a 2D setup. Numerical experiments have demonstrated that the location of the outer radius has no effect on the global evolution provided that the surface density at the outer boundary is not impacted by the formation and evolution of the inner cavity, namely provided it remains close to its initial value. In the azimuthal direction, the resolution is also set to be similar to that used in 3D.
Binary parameters employed in the simulations.
Fig. 1 Top: For Kepler34, 2D R–Z temperature distribution for the initial disc model. Bottom: 2D R–Z distribution of the dimensionless cooling time β = t_{cool}Ω_{k}. 
3.2 Initial conditions
The initial surface density profile in the 2D and 3D models is given by (8)
where Σ_{0} is the surface density at R_{0} = a_{bin} and is set to Σ_{0} = 1.3 × 10^{−4} in code units. The f_{gap} is a gap function used to initiate the disc with an inner cavity (assumed to be created by the binary) and is given by (9)
where R_{gap} = 2.5a_{bin} is the analytically estimated gap size (Artymowicz & Lubow 1994). The initial azimuthal velocity was set to the Keplerian velocity, whereas the initial radial velocity was set to zero. In the 3D runs, the meridional velocity was also set to zero.
The initial temperature, T_{0}, was determined by performing a Monte Carlo radiative transfer calculation of photon propagation using RADMC3D, in which multiple radiation sources can be incorporated. Here, we considered two radiation sources, whose stellar parameters correspond to the Kepler16 (AB) or Kepler34 (AB) systems. The stellar parameters that we employed for these two systems can be found in Table 1.
In RADMC3D, we assumed a grain size distribution that ranges from 0.1 µmto 10 µm with apowerlaw exponent of −3.5. The grains are 60% silicate and 40% amorphous carbon, with an internal density of 2.7 g cm^{−3}. The optical constants of silicate and amorphous carbon were obtained from Draine & Lee (1984) and Li & Greenberg (1997). For the 3D models, we first used the initial 3D gas density distribution to determine dust densities, assuming a net dusttogas ratio of 1%. RADMC3D then computed the dust temperature, and, since micrometer size grains tend to be strongly coupled to the gas, we can safely consider the gas and dust temperatures to be equal. The T_{0} can therefore be simply obtained by performing an azimuthal average of the calculated gas temperature. We illustrate this procedure by showing the resulting temperature structure for the Kepler34 analogue simulation calculated by RADMC3D in the upper panel of Fig. 1. The bottom panel of Fig. 1 displays the 2D (R–Z) distribution of the dimensionless cooling time β = t_{cool}Ω_{k}. At the location of the inner cavity, we see that β ≈ 100, meaning the disc behaves almost adiabatically, whereas it is close to isothermal near the disc surfaces. For the 2D models, we first assumed vertical hydrostatic equilibrium to deduce the 3D gas density distribution. We then applied a method similar to that described in the context of 3D models to get the initial gas temperature, but with the exception that in 2D it is defined as the vertically integrated temperature distribution (e.g. Bae et al. 2019): (10)
From this previous relation, the initial disc aspect ratio can be computed and is given by (11)
where ℛ is the perfect gas constant, µ is the mean molecular weight, which is set to µ = 2.4, and υ_{k} the Keplerian velocity. The initial surface density, temperature, and aspect ratio are displayed in Fig. 2. The factor of ~2 difference that is observed between the initial temperatures at the inner edge of the disc is in line with expectations. The received energy per unit time at this location is given by (Bitsch et al. 2013) (12)
where F_{⋆} is the stellar flux and can be approximated by (13)
where T_{1}, T_{2} (R_{1}, R_{2}) are the temperatures (radii) of the primary and secondary stars, respectively. On the other hand, the cooling at the surface of the inner disc edge is given by (14)
where T_{D} is the disc temperature. In thermodynamical equilibrium, we therefore have (15)
From this equation, and using the stellar parameters given in Table 1, it is straightforward to recover the factor of 2 shift between the disc temperatures of the Kepler16 and Kepler34 systems that are observed in Fig. 2.
For Kepler34, the temperature at the inner edge of the disc is likely to be higher than the threshold temperature for magnetorotational turbulence through thermal ionisation, which is estimated to be T_{MRI} ≈ 800 − 1200 K (e.g. Flock et al. 2016). For Kepler16, parametric instabilities may also give rise to turbulence, with an a viscous stress parameter (Shakura & Sunyaev 1973) α ≈ 5 × 10^{−3} (Pierens et al. 2020). Regardless of the origin of turbulence, the anomalous viscosity in the disc is modelled using the standard alpha prescription for the effective kinematic viscosity ν = αc_{s}H, where in this work, we used values of α = 10^{−3}, 10^{−2}, and 10^{−1}. Although a value of α = 10^{−1} is probably too large for a standard protoplanetary disc, it is comparable to the value expected in circumbinary discs around binary black holes (Moody et al. 2019; Muñoz et al. 2019; Tiede et al. 2020). Here, the calculation with α = 0.1 serves for us as a benchmark simulation used to validate the results from other runs, since it is expected to converge more rapidly due to a shorter viscous timescale. It is important to note that we do not incorporate in the energy equation the term corresponding to viscous heating so that the vertical thermal structure is similar for a given set of binary parameters. Our aim here is to examine the effect of viscosity on angular momentum transport rather than on the disc thermodynamics. Nevertheless, we discuss the effect of including viscous heating in Appendix A.
Fig. 2 Initial radial profiles of the surface density (top panel), temperature (middle panel), and aspect ratio (bottom panel). 
Fig. 3 Azimuthally averaged surface density (left panels), disc eccentricity (middle panels), and temperature (right panels) for the 2D and 3D Kepler16 analogues. 
Fig. 4 Temporal evolution of the disc mass for the Kepler16 simulations. Solid (respectively dashed) lines correspond to 3D (respectively 2D) models. 
3.3 Boundary conditions
We adopted an outflow boundary condition at the inner edge to allow mass to accrete onto the binary. For an outflow boundary condition, the velocity was set to 0 if directing towards the disc, such that gas material can leave the computational grid but no mass can flow into the domain. The values for the gas density and internal energy in the ghost zones have the same values as in the first active zone. In the context of 2D models, it has been shown that employing an open boundary together with a location of the inner disc edge R_{in} ≈ a_{bin} leads to a quasistationary disc structure after 𝒪(10^{4}) binary orbits. Compared to other boundary conditions, adopting an open boundary also appears to make the model less sensitive to numerical issues (Thun et al. 2017).
At the outer radial boundary, we employed a wavekilling zone for R > 0.88 R_{out} to avoid wave reflection or excitation at the disc’s outer edge (de ValBorro et al. 2006). The impact of this wavekilling zone on the disc shape and eccentricity is expected to be small since, as we mentioned earlier, the disc structure near the wavekilling zone remains close to the initial one.
At the meridional boundaries of the 3D domain, an outflow boundary condition was also employed, but with the exception that for the gas density we followed Bae et al. (2016) and maintained vertical stratification by solving the following condition for (16)
4 Results
4.1 Kepler16
In this section we describe the results for the Kepler16 analogue models, the parameters for which are given in Table 1.
4.1.1 Comparison between 3D and 2D models
In Fig. 3, we compare the surface density, disc eccentricity, and temperature profiles for the 2D and 3D models. The computational expense of running 3D simulations does not allow simulations to run for evolution times >2000 T_{bin}, so we restricted the comparison to times t ≤ 2000 T_{bin}. For the 3D runs, the disc eccentricity e_{d}(R) at radius R was defined as a densityweighted average, obtained by averaging the quantity e_{c} over θ and φ, where e_{c} is the eccentricity of a fluid element computed at the centre of each grid cell: (17)
where dV is the volume of one grid cell and where the integral over θ is performed over the total vertical extent of the disc. For the 2D models, e_{d}(R) is instead calculated following Pierens & Nelson (2013) and is given by (18)
where dS the surface area of one grid cell.
Consistent with previous work (Pierens & Nelson 2013; Kley et al. 2019; Penzlin 2021), the 2D models show larger and more eccentric inner cavities as α is increased. Interestingly, we find the opposite behaviour in the 3D models. Although the location of the maximum in the surface density seems similar for α = 10^{−3} and α = 10^{−2} in the 3D models, comparing the disc eccentricity panels reveals a higher disc eccentricity when α = 10^{−3}. Increasing α up to α = 0.1 leads to a smaller cavity size, while growth of the disc eccentricity is not observed at all in this case.
Due to higher cavity eccentricities, the 2D models tend to lose more mass through the inner boundary, since gas is more prone to approach the inner edge of the disc in that case. This is illustrated in Fig. 4, where we plot the evolution of the disc mass enclosed within the computational domain; both the 3D simulations and 2D models are shown. The increase in the disc mass observed in the 3D simulation with α = 0.1 is a consequence of (i) the wavekilling zone employed near the outer boundary, which acts as a mass reservoir because at each radius the radial velocity in the disc is towards the star, plus (ii) the fact there is no growth in the disc eccentricity for this model, meaning the inner disc material does not approach close to the inner boundary. A similar effect is not found in 2D due to the higher value for R_{out} which results in a longer viscous timescale at the outer boundary.
In Fig. 5, we present images of the surface density distributions for the different models, after an evolution time of 2000 T_{bin} when α = 10^{−3} and 10^{−2}, and 1500 T_{bin} when α = 0.1. The figure also shows the ellipse that best fits the structure of the inner cavity. The parameters of the ellipse are calculated following Thun et al. (2017), who assumed that the focus of the ellipse is located at the centre of mass of the binary and that the apocentre is in the direction where the surface density has its maximum value and at the location where the surface density is 10% of the maximum value. The values reported in the figure for the cavity semimajor axis, a_{gap}, and eccentricity, e_{gap}, show that the cavity structures in 2D and 3D are very similar when α = 10^{−3}; they also confirm that increasing α produces smaller cavities and disc eccentricities in 3D, whereas in 2D it produces larger cavity sizes and eccentricities (at least after a run time of 2000 T_{bin}). One implication of this trend is that for α > 0.01, not only the size but also the eccentricity of the inner cavity are smaller in 3D than in 2D, with the important consequence that migrating circumbinary planets should park closer to the binary when adopting a 3D setup and these higher values of α.
The time evolution of a_{gap} and e_{gap} over 2000 binary orbits is plotted for each model in Fig. 6. Both quantities show significant fluctuations because the procedure of finding the apocentre of the cavity may be misguided by local density maxima, as noticed by Penzlin (2021). Regarding the longterm trends, however, it is expected that the time to reach equilibrium should be shorter when α = 0.1 compared to the lower values, and it does indeed seem to be the case that the corresponding 3D run has converged after ~10^{3} T_{bin}. At that time, both a_{gap} and e_{gap} appear to have reached saturated values. For lower values of α, however, it is difficult to assess whether or not the 3D calculations have really converged after the total simulated time of 2000 orbits. For α = 10^{−3} it appears that a continuing evolution of a_{gap} and e_{gap} is occurring, whereas for α = 10^{−2} the values appear to be steadier on average towards the end of the run time. Despite the ambiguity about a steady state being achieved in the 3D runs, it is clear that there is a general tendency for the values of a_{gap} and e_{gap} obtained in the 3D runs to trend below the values obtained in 2D; additionally, for α ≥ 10^{−2} the cavities obtained in 3D are significantly smaller and more circular compared to those obtained in the 2D simulations. This can also be seen in Fig. 7, where we compare for each value of α the disc surface density and eccentricity profiles obtained at the end of the 3D runs with the corresponding 2D ones obtained at the same time.
Fig. 5 Surface density distributions for the 2D and 3D Kepler16 analogues for different α parameter. The snapshots correspond to simulation run times of 2000 T_{bin} when α = 10^{−3} and 10^{−2}, and 1500T_{bin} when α = 0.1. 
Fig. 6 Time evolution of the inner cavity semimajor axis, a_{gap}, and eccentricity, e_{gap}, for each Kepler16 analogue simulation. 
Fig. 7 For Kepler16 and for each value of α we considered, comparison of the surface density (top panel) and eccentricity (bottom panel) profiles obtained in the 3D (solid lines) and 2D (dashed lines) models. 
Fig. 8 Time evolution, over longer run times of 2 × 10^{4} binary orbits, of the innercavity semimajor axis, a_{gap}, and eccentricity, e_{gap}, for the 2D Kepler16 simulations. 
4.1.2 Long timescale evolution of the 2D runs
Previous work has shown that, in the context of a 2D setup, circumbinary discs typically settle to a stationary state after 𝒪(10^{4}) binary orbits. It can even require 𝒪(10^{5}) binary orbits to reach equilibrium in radiative discs, as shown recently by Kley et al. (2019) and Sudarshan et al. (2022). Our thermodynamics treatment of the disc in this work is different in comparison to these previous studies because the main source of heating is stellar irradiation rather than viscous heating. However, it is clear that a total evolution time of 2000 T_{bin} is not enough to allow the 2D discs to settle to a stationary state. For this set of runs, the time evolution of (a_{gap}, e_{gap}) over a longer simulated timescale of 2 × 10^{4} binary orbits is presented in Fig. 8. Not surprisingly, the model with α = 0.1 reaches a quasistationary state earlier compared to the two other runs, for which both a_{gap} and e_{gap} are still evolving at 2 × 10^{4} binary orbits. We note, however, that despite reaching a quasisteady state, the 2D run with α = 0.1 shows strong fluctuations in the computed values of a_{gap} and e_{gap}. Inspection of the results indicates that this arises because of the way that the algorithm used to calculate these quantities first identifies the location of the grid cell containing the maximum density. This is found to fluctuate strongly in this run, possibly because the cavity is relatively small and has a lower eccentricity. Hence, the fluctuations are a consequence of the algorithm somewhat misbehaving in this case. Since our primary intention here is to focus on 3D models, we did not explore the subsequent circumbinary disc evolution until a quasistationary state is reached for these 2D runs. Nevertheless, we can make use of the results of past studies on radiative circumbinary disc evolution to guess the evolution outcome in our 2D models with α = 10^{−3} and 10^{−2}. The case α = 10^{−3} has been examined by Sudarshan et al. (2022) using a βcooling prescription. For β ≈ 100, which, as mentioned in Sect. 3.2, corresponds to the value for β in the midplane of our disc models, these authors find a_{gap} = 4.35 a_{bin} and e_{gap} = 0.28 at equilibrium. We see that despite a slightly different treatment of the thermal physics, we obtain values for the cavity semimajor axis and eccentricity that are in rather good agreement with the ones reported by Sudarshan et al. (2022). The case with α = 10^{−2} has been considered by Kley et al. (2019), but only for a binary eccentricity e_{bin} = 0.1. For this value of e_{bin}, they And a_{gap} ≈ 4 a_{bin} and e_{gap} ≈ 0.38 once a quasisteady state is reached. Although we obtain a disc eccentricity similar to that found by Kley et al. (2019), we find a slightly higher value for the cavity semimajor axis, a_{gap} ≈ 5. This discrepancy probably arises because of the different values for e_{bin} that have been used, and/or possibly because we omit viscous heating in our simulations.
Fig. 9 Azimuthally averaged surface density (left panels), disc eccentricity (middle panels), and temperature (right panels) for the 2D and 3D models in the case of the Kepler34 system. 
Fig. 10 For Kepler34, 2D surface density distributions for the 2D and 3D models and for different α parameter values. These snapshots correspond to a simulated time of 2000 T_{bin} for calculations with α = 10^{−3}–10^{−2}, and of 1500 T_{bin} for the run with α = 0.1. 
Fig. 11 For Kepler34, time evolution of the inner cavity semimajor axis, a_{gap}, and eccentricity, e_{gap}, for each model. 
Fig. 12 For Kepler34 and for each value of α we considered, comparison of the surface density (top panel) and eccentricity (bottom panel) profiles obtained in the 3D (solid lines) and 2D (dashed lines) models. 
4.2 Kepler34
We now turn our attention to the circumbinary disc evolution for the Kepler34 analogue system (see Table 1 for the binary parameters). In Fig. 9, we show the surface density, disc eccentricity, and temperature profiles for the 2D and 3D models at times ≤2000 T_{bin}. Contrary to Kepler16 simulations, growth of disc eccentricity is now observed in the 3D run with α = 0.1, although the disc eccentricity is smaller by a factor of ~2 in comparison with its 2D counterpart. This tendency of the disc eccentricity being smaller in 3D is also verified in the simulations with lower values of α. Hence, as with the Kepler16 analogues, we find that the cavity sizes are systematically larger and more eccentric in 2D than in 3D. Looking at the 2D density distributions and the elliptical fits to the shape of the inner cavity in Fig. 10 confirms this tendency. Again, this would imply that migrating protoplanets should park closer to the central binary when adopting a 3D setup, but this needs to be checked in the future using dedicated simulations of embedded planets in 3D circumbinary disc models. Although we expect planets around highly eccentric binaries to stop their migration slightly beyond the location corresponding to the cavity semimajor axis (Penzlin 2021), it is interesting to notice that for α ≲ 0.01, the values obtained for a_{gap} in the 3D models are close to the observed semimajor axis of Kepler34b, which is a_{p} ≈ 4.7 a_{bin}.
With the exception of the run with α = 0.1, the values reported for the ellipse parameters in Fig. 10 suggest that a_{gap} and e_{gap} increase with α in both 3D and 2D discs, contrary to what has been obtained in the Kepler16 simulations. This similar behaviour between the 2D and 3D models can also be seen in Fig. 11, where we plot a_{gap} and e_{gap} as a function of time. Compared to the Kepler16 simulations, here it appears more evident that the 3D models reach a quasistationary state earlier than their 2D counterparts. While it is not possible to verify beyond all doubt that the 3D disc models have reached a quasisteady state for the runs with α < 10^{−1}, it is clear that for all values of α the cavity semimajor axes and eccentricities obtained in the 3D runs trend significantly below the values obtained in the corresponding 2D runs. Hence, the conclusion that the cavity sizes and eccentricities obtained in 3D are smaller than those obtained in 2D appears to be robust. For Kepler34, we compare in Fig. 12 the disc surface density and eccentricity profiles obtained at the end of the 3D runs with the 2D profiles at the same time. The fact that the cavity sizes and disc eccentricities tend to be smaller in 3D is clearly visible here.
For the 2D runs, the longterm evolution of the cavity semimajor axis and eccentricity is presented in Fig. 13. Similar to the discs around the Kepler16 analogues, an equilibrium state has not been reached after an evolution time of 2 × 10^{4} T_{bin}. At that time, the cavity size is estimated to be a_{gap} ~ 8a_{bin} for α = 10^{−2}, whereas we find a_{gap} ~ 7a_{bin} for α = 10^{−3}. This latter value is slightly higher than the one reported by Sudarshan et al. (2022), who find a_{gap} ~ 6a_{bin} for similar parameters in their simulations using a β prescription for cooling. This discrepancy possibly arises because we omit viscous heating, resulting in a smaller viscosity ν = αc_{s}H because of the lower temperature. As a consequence, the viscous torques that tend to counteract the gravitational torques from the binary are smaller in our simulations, leading to larger cavity sizes.
Overall, our results show that there is a trend for 3D models to reach equilibrium earlier than 2D models, and for the cavity sizes and eccentricities in 3D models to be smaller than in their 2D counterparts.
Fig. 13 For Kepler34, time evolution over 2 × 10^{4} binary orbits of the inner cavity semimajor axis, a_{gap}, and eccentricity, e_{gap}, for each 2D model. 
Fig. 14 Radial profile of the dimensionless cooling time for the initial disc models. 
Fig. 15 Top panel: surface density perturbations in the 3D and 2D discs using the perturbing potential given in Eq. (19) with A = 0.004. Bottom panel: same but for a perturbation amplitude A = 0.04. 
5 Discussion
In this section, we investigate the origin of the differences that are observed between the 2D and 3D models. We begin the discussion by emphasising that the spiral waves excited by the binary carry angular momentum that will be deposited in the disc to induce the formation of the inner cavity. Here, the larger cavity sizes obtained in the 2D simulations suggest that the angular momentum flux (AMF) carried by the spiral waves is larger in 2D than in 3D. This might occur, for example, if the 3D spiral waves propagating in the disc are subject to a more efficient radiative damping than in 2D discs. In the context of singlestar systems, previous work has demonstrated that the amplitude of the spirals depends significantly on the disc cooling timescale (Miranda & Rafikov 2020a,b; Zhang & Zhu 2020), and that the AMF may not be conserved under some circumstances. In particular, spiral waves tend to become weaker as the cooling time increases from a value corresponding to the isothermal limit to t_{cool} ~ Ω^{−1}. Spiral waves then become stronger again as the cooling time is increased from t_{cool} ~ Ω^{−1} to a value corresponding to the adiabatic limit. In adiabatic discs, the AMF is expected to be approximately constant (Miranda & Rafikov 2020a).
For our 2D models, the radial profile of the dimensionless cooling time, β, is shown in Fig. 14. We note that it also corresponds to the densityweighted vertically averaged t_{cool} in our 3D discs. At the edge of the inner cavity, β ≈ 10, meaning the disc behaves almost adiabatically. We therefore expect the AMF associated with waves excited there to be conserved in our 2D simulations. Going back to Fig. 1, however, we see that for the 3D models, β continuously decreases from β ≈ 100 in the equatorial plane to β ≈ 10^{−3} at the disc surfaces, so it is not guaranteed that the AMF is conserved in 3D.
To examine this issue in more detail, we studied the disc response to a perturbing potential of the form (19)
where A is the amplitude of the perturbing potential, for which we considered values of A = 0.004 and 0.04. We note that for this additional set of runs, we used the same thermodynamical treatment as for the simulations of circumbinary discs presented above. The results of these simulations are presented in Fig. 15. In the linear regime with A = 0.004, we see that the disc response is similar for the 2D and 3D models. In the nonlinear regime with A = 0.04, however, the spiral appears weaker in 3D, which indicates a higher (radiative) dissipation of the spiral wave. In the context of the simulations of circumbinary discs presented in Sect. 4, this would result in spiral waves with smaller amplitudes and consequently in smaller cavity sizes and eccentricities in 3D, in line with what has been observed. A more rapid dissipation of waves in the 3D calculations is also revealed by comparing contours of the surface density in Figs. 5 and 10, where we see that spiral waves indeed tend to be stronger in the 2D calculations.
In summary, it appears that the damping on nonlinear spiral waves occurs more rapidly in the 3D models compared to the 2D discs, and the origin of this is the different thermodynamical evolution in 3D versus 2D.
6 Summary
In this paper, we have presented the results of 3D hydrodynamical simulations of circumbinary discs around analogues of the Kepler16(AB) and Kepler34(AB) binary systems. We included stellar heating from the central binary, and radiative cooling by employing a cooling term in which the temperature is relaxed towards the irradiation temperature on a local cooling timescale, t_{cool}. The irradiation temperature was calculated by performing a Monte Carlo radiative transfer calculation of photon propagation using RADMC3D, whereas the cooling timescale, t_{cool}, was determined using realistic opacities.
We find that 3D circumbinary disc models show a tendency to reach a quasistationary state more rapidly than their 2D counterparts. Verification of this statement for all models that we have considered is difficult at the present time because of the computational expense of running 3D hydrodynamical simulations. However, for a subset of the parameters we have considered, the more rapid convergence of 3D simulations versus their 2D counterparts is beyond doubt. Furthermore, the sizes and eccentricities of the inner cavity are also found to be systematically smaller in 3D models compared to the equivalent 2D models. We interpret this result to be a consequence of the more rapid damping of the nonlinear spiral waves that are excited at the cavity edge in 3D versus 2D, or, equivalently, to be due to the enhanced nonconservation of the AMF associated with these waves as they travel into the disc. Although the disc behaves almost adiabatically in our 2D calculations, which would result in a constant AMF for a linear wave that does not undergo nonlinear damping, the 3D structure of the disc contains regions located at intermediate scale heights with t_{cool} ~ Ω^{−1}, and where the spiral wave can be subject to significant radiative damping. Compared to the 2D case, this appears to result in smaller inner cavities that are less eccentric in 3D calculations compared to 2D simulations.
A potentially important consequence of these results is that we would expect migrating protoplanets to park closer to the central binary when evolving in a 3D disc model compared to a 2D one, and to maintain a smaller eccentricity. In the case of the circumbinary planet Kepler34b, this might allow simulations to reproduce the observed orbital parameters of the planet, which has not been possible using 2D simulations (e.g. Pierens & Nelson 2013; Penzlin 2021). Determining whether or not this is the case will be the subject of future work.
Acknowledgements
Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and by HPC resources of Cines under the allocation A0110406957 made by GENCI (Grand Equipement National de Calcul Intensif). R.P.N. acknowledges support from Leverhulme Trust through grant number RPG2018418 and from STFC through grants ST/P000592/1 and ST/T000341/1.
Appendix A Effect of viscous heating
As mentioned in Sect. 3.2, we do not take the effect of viscous heating into account in this study. To estimate the importance of viscous heating relative to stellar heating in a simple way, we included viscous heating in RADMC3D and computed the temperature from the initial density field. For α = 0.01 and binary parameters corresponding to Kepler34, contours of the outputted temperature are presented in the top panel of Fig. A.1. Comparing this figure with Fig. 1, we see that including viscous heating leads to higher temperatures in the inner regions and in the disc midplane as well, which is not surprising as the strength of viscous heating scales with the value of the gas density. The bottom panel of Fig. A.1 shows the radial profile of the temperature. Including viscous heating typically leads to temperatures higher by a factor of ~ 2 compared to the case of stellar heating only. This would result in smaller cavity sizes in simulations where the effect of viscous heating is taken into account.
Fig. A.1 Top: For Kepler34, 2D R–Z temperature distribution for the initial disc model in the case where the effect of viscous heating is considered. Bottom: Corresponding radial temperature profile. Here, we assumed α = 0.01. 
Rerences
 Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [Google Scholar]
 Bae, J., Nelson, R. P., & Hartmann, L. 2016, ApJ, 833, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJ, 884, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Barker, A. J., & Ogilvie, G. I. 2014, MNRAS, 445, 2637 [NASA ADS] [CrossRef] [Google Scholar]
 BenítezLlambay, P., & Masset, F. S. 2016, ApJS, 223, 11 [Google Scholar]
 Bitsch, B., Crida, A., Morbidelli, A., et al. 2013, A&A, 549, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Coleman, G. A. L., Nelson, R. P., & Triaud, A. H. M. J. 2022, MNRAS, 513, 2563 [NASA ADS] [CrossRef] [Google Scholar]
 de ValBorro, M., Edgar, R.G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library. [ascl:1202.015] [Google Scholar]
 Flock, M., Fromang, S., Turner, N. J., et al. 2016, ApJ, 827, 144 [CrossRef] [Google Scholar]
 Georgakarakos, N., Eggl, S., & DobbsDixon, I. 2021, Front. Astron. Space Sci., 8, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [Google Scholar]
 Kley, W., & Haghighipour, N. 2014, A&A, 564, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., & Haghighipour, N. 2015, A&A, 581, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., Thun, D., & Penzlin, A. B. T. 2019, A&A, 627, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kostov, V. B., Orosz, J. A., Feinstein, A. D., et al. 2020, AJ, 159, 253 [Google Scholar]
 Kostov, V. B., Powell, B. P., Torres, G., et al. 2021, ApJ, 917, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Li, A., & Greenberg, J. M. 1997, A&A, 323, 566 [NASA ADS] [Google Scholar]
 Lines, S., Leinhardt, Z. M., Baruteau, C., et al. 2016, A&A, 590, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lyra, W., Richert, A. J. W., Boley, A., et al. 2016, ApJ, 817, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., Thébault, P., & Scholl, H. 2008, ApJ, 681, 1599 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. S., Morbidelli, A., Crida, A., et al. 2006, ApJ, 642, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Miranda, R., & Rafikov, R. R. 2020a, ApJ, 892, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Miranda, R., & Rafikov, R. R. 2020b, ApJ, 904, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Moody, M. S. L., Shi, J.M., & Stone, J. M. 2019, ApJ, 875, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, T. W. A., & Kley, W. 2012, A&A, 539, A18 [Google Scholar]
 Muñoz, D. J., Miranda, R., & Lai, D. 2019, ApJ, 871, 84 [CrossRef] [Google Scholar]
 Muñoz, D. J., Lai, D., Kratter, K., et al. 2020, ApJ, 889, 114 [Google Scholar]
 Mutter, M. M., Pierens, A., & Nelson, R. P. 2017, MNRAS, 465, 4735 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B. 2005, A&A, 432, 743 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Penzlin, A., Kley, W., & Nelson, R. P. 2021, A&A, 645, 68 [Google Scholar]
 Pierens, A., & Nelson, R. P. 2007, A&A, 472, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2013, A&A, 556, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2018, MNRAS, 477, 2547 [NASA ADS] [CrossRef] [Google Scholar]
 Pierens, A., McNally, C. P., & Nelson, R. P. 2020, MNRAS, 496, 2849 [CrossRef] [Google Scholar]
 Pierens, A., Nelson, R. P., & McNally, C. P. 2021, MNRAS, 508, 4806 [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
 Sudarshan, P., Penzlin, A. B. T., Ziampras, A., et al. 2022, A&A, 664, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tiede, C., Zrake, J., MacFadyen, A., et al. 2020, ApJ, 900, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Thun, D., Kley, W., & Picogna, G. 2017, A&A, 604, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Triaud, A. H. M. J. 2021, Plato Mission Conference 2021, 50 [Google Scholar]
 Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475 [Google Scholar]
 Zhang, S., & Zhu, Z. 2020, MNRAS, 493, 2287 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Hartmann, L., & Gammie, C. 2009, ApJ, 694, 1045 [Google Scholar]
All Tables
All Figures
Fig. 1 Top: For Kepler34, 2D R–Z temperature distribution for the initial disc model. Bottom: 2D R–Z distribution of the dimensionless cooling time β = t_{cool}Ω_{k}. 

In the text 
Fig. 2 Initial radial profiles of the surface density (top panel), temperature (middle panel), and aspect ratio (bottom panel). 

In the text 
Fig. 3 Azimuthally averaged surface density (left panels), disc eccentricity (middle panels), and temperature (right panels) for the 2D and 3D Kepler16 analogues. 

In the text 
Fig. 4 Temporal evolution of the disc mass for the Kepler16 simulations. Solid (respectively dashed) lines correspond to 3D (respectively 2D) models. 

In the text 
Fig. 5 Surface density distributions for the 2D and 3D Kepler16 analogues for different α parameter. The snapshots correspond to simulation run times of 2000 T_{bin} when α = 10^{−3} and 10^{−2}, and 1500T_{bin} when α = 0.1. 

In the text 
Fig. 6 Time evolution of the inner cavity semimajor axis, a_{gap}, and eccentricity, e_{gap}, for each Kepler16 analogue simulation. 

In the text 
Fig. 7 For Kepler16 and for each value of α we considered, comparison of the surface density (top panel) and eccentricity (bottom panel) profiles obtained in the 3D (solid lines) and 2D (dashed lines) models. 

In the text 
Fig. 8 Time evolution, over longer run times of 2 × 10^{4} binary orbits, of the innercavity semimajor axis, a_{gap}, and eccentricity, e_{gap}, for the 2D Kepler16 simulations. 

In the text 
Fig. 9 Azimuthally averaged surface density (left panels), disc eccentricity (middle panels), and temperature (right panels) for the 2D and 3D models in the case of the Kepler34 system. 

In the text 
Fig. 10 For Kepler34, 2D surface density distributions for the 2D and 3D models and for different α parameter values. These snapshots correspond to a simulated time of 2000 T_{bin} for calculations with α = 10^{−3}–10^{−2}, and of 1500 T_{bin} for the run with α = 0.1. 

In the text 
Fig. 11 For Kepler34, time evolution of the inner cavity semimajor axis, a_{gap}, and eccentricity, e_{gap}, for each model. 

In the text 
Fig. 12 For Kepler34 and for each value of α we considered, comparison of the surface density (top panel) and eccentricity (bottom panel) profiles obtained in the 3D (solid lines) and 2D (dashed lines) models. 

In the text 
Fig. 13 For Kepler34, time evolution over 2 × 10^{4} binary orbits of the inner cavity semimajor axis, a_{gap}, and eccentricity, e_{gap}, for each 2D model. 

In the text 
Fig. 14 Radial profile of the dimensionless cooling time for the initial disc models. 

In the text 
Fig. 15 Top panel: surface density perturbations in the 3D and 2D discs using the perturbing potential given in Eq. (19) with A = 0.004. Bottom panel: same but for a perturbation amplitude A = 0.04. 

In the text 
Fig. A.1 Top: For Kepler34, 2D R–Z temperature distribution for the initial disc model in the case where the effect of viscous heating is considered. Bottom: Corresponding radial temperature profile. Here, we assumed α = 0.01. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.