© The Authors 2023

1 Introduction

To date, 15 circumbinary planets have been discovered around eclipsing binary systems that contain main sequence stars. Most of these, such as Kepler-16b and Kepler-34b, were detected by the Kepler spacecraft (Doyle et al. 2011; Welsh et al. 2012), and recently TESS has also detected two circumbinary planets, TOI-1338b (Kostov et al. 2020) and TIC-172900988 (Kostov et al. 2021). The first detection of a circumbinary planet via the radial velocity method recently confirmed that Kepler-16b has a mass similar to that of Saturn (Triaud et al. 2021), in agreement with the mass estimate obtained by fitting a photo-dynamical 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 Kepler-413 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 Kepler-413b 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 quasi-steady 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 Kepler-16(AB) and Kepler-34(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 $$\frac{\partial e}{\partial t}+\nabla \cdot \left(ev\right)=-\left(\gamma -1\right)e\nabla \cdot v+Q_{\text{cool}},$$(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 $$Q_{\text{cool}}=-\rho c_{\upsilon }\frac{T-T_{\text{irr}}}{t_{\text{cool}}}.$$(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 RADMC-3D 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 $$t_{\text{cool}}=\frac{\rho c_{\upsilon }H{\tau }_{\text{eff}}}{3\sigma T^3},$$(3)

where σ is the Stefan-Boltzmann constant and τ_eff is the effective optical depth given by $${\tau }_{\text{eff}}=\frac{3}{8}\tau +\frac{\sqrt{3}}{4}+\frac{1}{4\tau }.$$(4)

Following Bae et al. (2016), the optical depth, τ, was calculated as 1/τ = 1/τ_upper + 1/τ_lower with $${\tau }_{\text{upper}}={\displaystyle {\int }_z^{z_{\max }}\rho \left(z^{\prime }\right)\kappa \left(z^{\prime }\right)\text{d}{z}^{\prime }}$$(5)

and $${\tau }_{\text{lower}}={\displaystyle {\int }_{z_{\min }}^z\rho \left(z^{\prime }\right)\kappa \left(z^{\prime }\right)\text{d}{z}^{\prime }},$$(6)

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 $$\tau =\frac{c_1}{\sqrt{2\pi }}\kappa \text{Σ}.$$(7)

Here, c₁ 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 multi-fluid version of FARGO3D (Benítez-Llambay & 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 semi-major axis, a_bin. When presenting the simulation results, unless otherwise stated, we used the binary orbital period $T_{\text{bin}}=2\pi \sqrt{a_{\text{bin}}^3\text{/}G{M}_{\star }}$ 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.

Fig. 1 Top: For Kepler-34, 2D R–Z temperature distribution for the initial disc model. Bottom: 2D R–Z distribution of the dimensionless cooling time β = tcoolΩk.

3.2 Initial conditions

The initial surface density profile in the 2D and 3D models is given by $$\text{Σ}=f_{\text{gap}}{\text{Σ}}_0{\left(\frac{R}{R_0}\right)}^{-3/2},$$(8)

where Σ₀ is the surface density at R₀ = a_bin and is set to Σ₀ = 1.3 × 10⁻⁴ 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 $$f_{\text{gap}}={\left(1+\exp \left[-\frac{R-R_{\text{gap}}}{0.1R_{\text{gap}}}\right]\right)}^{-1},$$(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₀, was determined by performing a Monte Carlo radiative transfer calculation of photon propagation using RADMC-3D, in which multiple radiation sources can be incorporated. Here, we considered two radiation sources, whose stellar parameters correspond to the Kepler-16 (AB) or Kepler-34 (AB) systems. The stellar parameters that we employed for these two systems can be found in Table 1.

In RADMC-3D, we assumed a grain size distribution that ranges from 0.1 µmto 10 µm with apower-law exponent of −3.5. The grains are 60% silicate and 40% amorphous carbon, with an internal density of 2.7 g cm⁻³. 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 dust-to-gas ratio of 1%. RADMC-3D 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₀ 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 Kepler-34 analogue simulation calculated by RADMC-3D 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): $$T\left(R\right)=\frac{\int \rho \left(R,Z\right)T\left(R,Z\right)\text{d}Z}{\int \rho \left(R,Z\right)\text{d}Z}.$$(10)

From this previous relation, the initial disc aspect ratio can be computed and is given by $$h={\left(\frac{ℛT\left(R\right)}{\mu }\right)}^{1/2}{\upsilon }_k^{-1},$$(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) $$Q_{+}=2\pi R_{\text{in}}2H{F}_{*},$$(12)

where F_⋆ is the stellar flux and can be approximated by $$F_{*}=\sigma \frac{R_1^2{T}_1^4+R_2^2{T}_2^4}{R_{\text{in}}^2},$$(13)

where T₁, T₂ (R₁, R₂) 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 $$Q_{\text{cool}}=-4\pi R_{\text{in}}H\sigma T_{\text{D}}^4,$$(14)

where T_D is the disc temperature. In thermodynamical equilibrium, we therefore have $$T_{\text{D}}={\left(\frac{R_1^2{T}_1^4+R_2^2{T}_2^4}{R_{\text{in}}^2}\right)}^{1/4}.$$(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 Kepler-16 and Kepler-34 systems that are observed in Fig. 2.

For Kepler-34, 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 Kepler-16, parametric instabilities may also give rise to turbulence, with an a viscous stress parameter (Shakura & Sunyaev 1973) α ≈ 5 × 10⁻³ (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_sH, where in this work, we used values of α = 10⁻³, 10⁻², and 10⁻¹. Although a value of α = 10⁻¹ 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 Kepler-16 analogues.

Fig. 4 Temporal evolution of the disc mass for the Kepler-16 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 quasi-stationary disc structure after 𝒪(10⁴) 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 wave-killing zone for R > 0.88 R_out to avoid wave reflection or excitation at the disc’s outer edge (de Val-Borro et al. 2006). The impact of this wave-killing zone on the disc shape and eccentricity is expected to be small since, as we mentioned earlier, the disc structure near the wave-killing 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 $$\frac{1}{\rho }\frac{\partial }{\partial \theta }\left(c_{\text{s}}^2\rho \right)=\frac{{\upsilon }_{\varphi }^2}{\tan \theta }.$$(16)

4 Results

4.1 Kepler-16

In this section we describe the results for the Kepler-16 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 density-weighted 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: $$e_{\text{d}}\left(R\right)=\frac{{\int }_{{\theta }_{-}}^{{\theta }_{+}}{\int }_0^{2\pi }{\int }_R^{R+\text{d}R}\rho e_{\text{c}}\text{d}V}{{\int }_{{\theta }_{-}}^{{\theta }_{+}}{\int }_0^{2\pi }{\int }_R^{R+\text{d}R}\rho \text{d}V},$$(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 $$e_{\text{d}}\left(R\right)=\frac{{\int }_0^{2\pi }{\int }_R^{R+\text{d}R}\text{Σ}{e}_{\text{c}}\text{d}S}{{\int }_0^{2\pi }{\int }_R^{R+\text{d}R}\text{Σd}S},$$(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⁻³ and α = 10⁻² in the 3D models, comparing the disc eccentricity panels reveals a higher disc eccentricity when α = 10⁻³. 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 wave-killing 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⁻³ and 10⁻², 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 semi-major axis, a_gap, and eccentricity, e_gap, show that the cavity structures in 2D and 3D are very similar when α = 10⁻³; 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 long-term 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³ 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⁻³ it appears that a continuing evolution of a_gap and e_gap is occurring, whereas for α = 10⁻² 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⁻² 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 Kepler-16 analogues for different α parameter. The snapshots correspond to simulation run times of 2000 Tbin when α = 10−3 and 10−2, and 1500Tbin when α = 0.1.

Fig. 6 Time evolution of the inner cavity semi-major axis, agap, and eccentricity, egap, for each Kepler-16 analogue simulation.

Fig. 7 For Kepler-16 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 × 104 binary orbits, of the inner-cavity semi-major axis, agap, and eccentricity, egap, for the 2D Kepler-16 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⁴) binary orbits. It can even require 𝒪(10⁵) 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⁴ binary orbits is presented in Fig. 8. Not surprisingly, the model with α = 0.1 reaches a quasi-stationary state earlier compared to the two other runs, for which both a_gap and e_gap are still evolving at 2 × 10⁴ binary orbits. We note, however, that despite reaching a quasi-steady 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 quasi-stationary 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⁻³ and 10⁻². The case α = 10⁻³ 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 semi-major axis and eccentricity that are in rather good agreement with the ones reported by Sudarshan et al. (2022). The case with α = 10⁻² 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 quasi-steady 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 semi-major 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 Kepler-34 system.

Fig. 10 For Kepler-34, 2D surface density distributions for the 2D and 3D models and for different α parameter values. These snapshots correspond to a simulated time of 2000 Tbin for calculations with α = 10−3–10−2, and of 1500 Tbin for the run with α = 0.1.

Fig. 11 For Kepler-34, time evolution of the inner cavity semi-major axis, agap, and eccentricity, egap, for each model.

Fig. 12 For Kepler-34 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 Kepler-34

We now turn our attention to the circumbinary disc evolution for the Kepler-34 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 Kepler-16 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 Kepler-16 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 semi-major 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 semi-major axis of Kepler-34b, 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 Kepler-16 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 Kepler-16 simulations, here it appears more evident that the 3D models reach a quasi-stationary state earlier than their 2D counterparts. While it is not possible to verify beyond all doubt that the 3D disc models have reached a quasi-steady state for the runs with α < 10⁻¹, it is clear that for all values of α the cavity semi-major 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 Kepler-34, 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 long-term evolution of the cavity semi-major axis and eccentricity is presented in Fig. 13. Similar to the discs around the Kepler-16 analogues, an equilibrium state has not been reached after an evolution time of 2 × 10⁴ T_bin. At that time, the cavity size is estimated to be a_gap ~ 8a_bin for α = 10⁻², whereas we find a_gap ~ 7a_bin for α = 10⁻³. 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_sH 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 Kepler-34, time evolution over 2 × 104 binary orbits of the inner cavity semi-major axis, agap, and eccentricity, egap, 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 single-star 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 ~ Ω⁻¹. Spiral waves then become stronger again as the cooling time is increased from t_cool ~ Ω⁻¹ 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 density-weighted 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⁻³ 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 $$\text{Φ}\left(r,\phi \right)=\frac{A}{r^5}\cos \left(4\phi -\text{Ω}t\right),$$(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 non-linear 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 non-linear 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 Kepler-16(AB) and Kepler-34(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 RADMC-3D, whereas the cooling timescale, t_cool, was determined using realistic opacities.

We find that 3D circumbinary disc models show a tendency to reach a quasi-stationary 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 non-linear spiral waves that are excited at the cavity edge in 3D versus 2D, or, equivalently, to be due to the enhanced non-conservation 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 ~ Ω⁻¹, 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 Kepler-34b, 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 RPG-2018-418 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 RADMC-3D and computed the temperature from the initial density field. For α = 0.01 and binary parameters corresponding to Kepler-34, 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 Kepler-34, 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

All Tables

All Figures

	Fig. 1 Top: For Kepler-34, 2D R–Z temperature distribution for the initial disc model. Bottom: 2D R–Z distribution of the dimensionless cooling time β = tcoolΩ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 Kepler-16 analogues.
In the text

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

	Fig. 5 Surface density distributions for the 2D and 3D Kepler-16 analogues for different α parameter. The snapshots correspond to simulation run times of 2000 Tbin when α = 10−3 and 10−2, and 1500Tbin when α = 0.1.
In the text

	Fig. 6 Time evolution of the inner cavity semi-major axis, agap, and eccentricity, egap, for each Kepler-16 analogue simulation.
In the text

	Fig. 7 For Kepler-16 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 × 104 binary orbits, of the inner-cavity semi-major axis, agap, and eccentricity, egap, for the 2D Kepler-16 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 Kepler-34 system.
In the text

	Fig. 10 For Kepler-34, 2D surface density distributions for the 2D and 3D models and for different α parameter values. These snapshots correspond to a simulated time of 2000 Tbin for calculations with α = 10−3–10−2, and of 1500 Tbin for the run with α = 0.1.
In the text

	Fig. 11 For Kepler-34, time evolution of the inner cavity semi-major axis, agap, and eccentricity, egap, for each model.
In the text

	Fig. 12 For Kepler-34 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 Kepler-34, time evolution over 2 × 104 binary orbits of the inner cavity semi-major axis, agap, and eccentricity, egap, 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 Kepler-34, 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