© P. Sudarshan et al. 2022

1. Introduction

About 50% of all stars are part of a binary configuration. Binary stars are believed to form together often, such that they have a common circumbinary disc in which planets can form and migrate. The Kepler mission has revealed about a dozen close-in circumbinary planets (Doyle et al. 2011; Orosz et al. 2012, 2019; Welsh et al. 2012, 2015; Schwamb et al. 2013; Kostov et al. 2014, 2020; Socia et al. 2020). The circumbinary protoplanetary disc (CBD/PPD) is key to understanding the evolution of the binary stars and circumbinary planets during their formation.

The structure of the disc and the interaction between discs and embedded objects is a subject of ongoing research. In the context of binary super-massive black holes (SMBHs), accretion from the disc onto the binary has been studied using various models. The effect of the disc on a binary object was recently studied by Muñoz et al. (2019, 2020) and Tiede et al. (2020), discovering binary expansion for thick, viscous discs (α ≥ 0.01, h = 0.1). Tiede et al. (2021) looked into the accretion dynamics with tracer particles embedded in a Cartesian grid simulation. The effects of the binary eccentricity on the momentum exchange between the binary components was investigated in models by D'Orazio & Duffell (2021). Dittmann & Ryan (2021, 2022) explored the effect of the disc height and binary mass ratios. A 3D study by Moody et al. (2019) tested how misalignment between the disc and binary changes the accretion rate. In this context of SMBHs, while disc heights are comparable to PPDs, the viscosity is much higher, leading to a faster evolution of the models, but also making the handling of viscous heating in the disc even more relevant.

Mösta et al. (2019) studied the flow inside the inner eccentric cavity of a thick disc (h = 0.2) and found that the stream pattern inside the inner cavity is likely the product of spiral shocks. For PPD-like systems, the recent work by Thun et al. (2017), Hirsh et al. (2020), and Ragusa et al. (2020) shows large eccentric inner cavities in circumbinary discs, the properties of which depend on the disc and binary parameters. The migration of planets in circumbinary discs was studied in Thun & Kley (2018) and Penzlin et al. (2021) without self-gravity and by Mutter et al. (2017) with self-gravity.

All of the mentioned studies so far have used vertically integrated, locally isothermal models, where the cooling timescale is negligible compared to the orbital period of the gas, ignoring changes to the thermal profile of the disc. Kley et al. (2019) introduced a radiatively cooled model for the circumbinary disc with migrating planets and Pierens et al. (2020, 2021) used a three-dimensional (3D) model with β- cooling including dust to study the impact of the turbulent disc on in situ planet growth close to the binary.

These studies confirm large eccentric inner cavities that are also derived analytically in Muñoz & Lithwick (2020). One caveat to simulating the disc, however, is the long time needed to fully evolve this inner cavity to a convergent state. As a result, simulations that include many advanced physical models can be unfeasible due to runtime constraints. The locally isothermal prescription can provide significant speed-up, at the cost of excluding the thermal evolution of the disc and cavity, which can be critical to their equilibrium state. As a result, the effects of a thermal model on the disc have not been compared in detail.

For single-star protoplanetary discs, recent work by Miranda & Rafikov (2020a,b), Zhang & Zhu (2020), and Ziampras et al. (2020) has shown that the cooling timescale can change the gap opening capabilities of embedded planets and the density slope at their gap edge, due to radiative cooling interfering with the angular momentum transport efficiency of spiral arms launched by such planets. Such effects could become relevant for the circumbinary cavity edge and therefore for the inner disc shape.

In this paper, we investigate how the locally isothermal model compares to the simple approach of a parametrised, constant cooling timescale and a more complex radiatively cooled model in the context of PPDs. We also construct simplistic models with adaptive parametrised cooling that aim to reach a compromise between computational speed and physical accuracy.

In Sect. 2, the model setup is explained and the radiative effects considered are described. The results are shown in Sect. 3, where all models are compared and contrasted throughout our parameter space. We discuss our findings as well as possible caveats and improvements in Sect. 4. Finally, our results are summarised in Sect. 5.

2. Model setup

We solved the vertically integrated hydrodynamics equations in cylindrical polar coordinates for a disc of ideal gas with surface density ∑, vertically integrated pressure P and velocity u = (u_r, u_ϕ). We used the same equations that are laid out in Kley et al. (2019):

$$\begin{array}{lll}\frac{\partial \Sigma }{\partial t}+\nabla \cdot \left(\Sigma u\right)\hfill & =\hfill & 0,\hfill \\ \frac{\partial \Sigma u}{\partial t}+\nabla \cdot \left(\Sigma u\right)\otimes u+PI-\overline{II}\hfill & =\hfill & \Sigma g,\hfill \\ \frac{\partial e}{\partial t}+\nabla \cdot \left(\left(e+P-\overline{II}\right)u\right)\hfill & =\hfill & \Sigma u\cdot g+Q.\hfill \end{array}$$(1)

Here, $\overline{\prod }$ denotes the viscous stress tensor, g refers to gravitational acceleration by the binary stars, and Q encompasses any additional source terms. For an ideal gas with adiabatic index γ, the isothermal and adiabatic sound speeds are related as $c_{\text{s},iso}=c_{\text{s}}/\sqrt{\gamma }=\sqrt{P/\sum }=H{\text{Ω}}_{\text{K}}$, where H is the gas pressure scale height and ${\text{Ω}}_{\text{K}}=\sqrt{\text{G}{M}_{\text{co}}/r^3}$ the Keplerian angular velocity at distance r from the central object(s) with mass M_co. The gravitational constant is denoted with G. The ideal gas relation also dictates that the specific internal energy of the gas is given by e = P/(γ − 1). The temperature T is then

$$T=\frac{\mu }{R}{c}_{\text{s,iso}}^2.$$(2)

Here, µ = 2.353 is the mean molecular weight and ℛ is the gas constant. These equations hold true for all models presented here.

2.1. Temperature descriptions

In this work, we use three different prescriptions of radiative effects inside the disc: locally isothermal (no radiative effects), a parametrised cooling timescale, and self-consistent radiative cooling via thermal emission. In this section we describe each model in more detail.

2.1.1. Locally isothermal model

The locally isothermal prescription assumes that the temperature is a fixed radial profile, equivalent to a disc that responds instantly to heating or cooling and is always in thermal equilibrium. For a disc with an aspect ratio h = H/r, the temperature profile is given by

$$T=\frac{\mu G{M}_{\text{co}}}{R}\frac{h^2}{r}.$$(3)

We choose a constant h = 0.05, a value which is in good agreement with radiative simulations of binary systems by Kley et al. (2019).

2.1.2. Parametrised β-cooling

We use a β-cooling prescription similar to Gammie (2001), and follow the implementation detailed in Rometsch et al. (2021). The β relaxation model assumes that any temperature change is relaxed to a reference profile T₀ over a cooling timescale $t_{\text{cool}}=\beta {\text{Ω}}_{\text{K}}^{-1}$ via a source term Q_relax such that

$$Q_{\text{relax}}=-\Sigma c_v\frac{T-T_0}{t_{\text{cool}}}\Rightarrow \frac{\partial T}{\partial t}=-\frac{T-T_0}{\beta }{\Omega }_{\text{K}},$$(4)

where c_v = ^ℛ/µ(γ−1) is the heat capacity at constant volume. Assuming viscous heating with $Q_{\text{visc}}\approx \frac{9}{4}v\sum {\text{Ω}}_{\text{K}}^2$ (Tassoul 1978) and the α-viscosity model of Shakura & Sunyaev (1973), we have ν = αc_sH. Then, in a steady state with no compressional heating (∇ · u = 0) it follows that

$$\begin{array}{cc}{Q}_{\text{visc}}+Q_{\text{relax}}=0\Rightarrow T_{\text{eq}}=\frac{T_0}{1-k\alpha \beta },& k:=\frac{9}{4}\sqrt{\gamma }\left(\gamma -1\right)\end{array}$$(5)

with k ≈ 1.06 for γ = 7/5. In other words, the equilibrium temperature is higher for longer relaxation times or higher viscosity.

2.1.3. Radiative models

For radiative models we use the thermal cooling prescription of Kley et al. (2019). In a radiatively cooled system, the disc cools by emitting received heat as thermal radiation, the efficiency of which depends on the local environment. In these models, cooling is given by

$$\frac{\partial \Sigma c_{\text{v}}T}{\partial t}=-2{\sigma }_{\text{SB}}\frac{T^4}{{\tau }_{e\text{ff}}}=Q_{\text{cool}},$$(6)

where σ_SB is the Stefan-Boltzmann constant. The effective optical depth τ_eff can be calculated with the emission optical depth τ following Hubeny (1990) as

$$\begin{array}{cc}{\tau }_{e\text{ff}}=\frac{3}{8}\tau +\frac{\sqrt{3}}{4}+\frac{1}{4\tau +{\tau }_{\min }},& \tau ={\displaystyle \underset{0}{\overset{\infty }{\int }}{\kappa }_{\text{R}}\rho \text{ dz}\approx }\frac{c_1}{\sqrt{2\pi }}{\kappa }_{\text{R}}\Sigma .\end{array}$$(7)

The coefficient c₁ = 0.5 is obtained by comparing 2D disc models with fully 3D calculations following Kley et al. (2009). In this equation K_R is the Rosseland mean opacity as modelled by Lin & Papaloizou (1985), with the density and temperature power law relation detailed in Müller & Kley (2012). We use a minimum optical depth τ_min = 0.01 to account for the extremely optically thin regions of the disc.

We can estimate the cooling timescale via this cooling prescription following previous studies (e.g. Miranda & Rafikov 2020a; Ziampras et al. 2020) by writing

$$\frac{\partial e}{\partial t}~\frac{e}{t_{\text{cool}}}=\left|Q_{\text{cool}}\right|\Rightarrow \beta =\frac{\Sigma c_{\text{v}}T}{\left|Q_{\text{cool}}\right|}{\Omega }_{\text{K}}.$$(8)

Assuming thermal balance between viscous heating and thermal cooling, we can write that |Q_cool| = Q_visc and find that αβ = 1/k or β ≈ 1/α. As a result, we can estimate that in the absence of shock heating, the cooling timescale in our radiative model should beg β ≈ 10³. As we show later, this is indeed the case for the majority of the disc down to the cavity edge. Within the cavity region shock heating dominates due to proximity to the binary, while the very low gas density diminishes the efficiency of viscous heating which would modify the equilibrium state to |Qcool| = Q_visc − P∇ · u. For this reason, models where β is constant run the risk; of not treating cooling within the cavity region properly, especially tor larger values of β. We address this issue in Secs. 3.4;

Fig. 1 Time evolution of the 2D surface density normalised to its reference power-law profiles ∑re∝ r−0.5 for the locally isothermal test case with ebin = 0.15. The system eventually reaches an equilibrium by t = 50 000 Tbin. The dashed green curves represent the ellipse fitting. The binary spiral streams, which are present inside the cavity, are not visible within this colour range, but can be seen in Fig. 7.

2.2. Disc setup

We chose a binary setup that is comparable to the Kepler circumbinary planet-hosting systems discovered by the Kepler-telescope. The total mass of the binary is M_bin = M_co = 1 M_⊙, the binary semi-major axis is a_bin = 0.2 au, and She mass ratio is q = M₂/M₁ = 0.5. We vary the binary eccentricity e_bin ϵ (0.0, 0.15, 0.3, 0.5) to cover both e_bin-dependent branches and the bifurcation point in the cavity shape for circumbinary discs found by Thun & Kley (2018). The binary does not receive feedback from the disc, maintaining a fixed orbit around its common centre of mass. The total initial mass of the gas is 1.7% of M_bin.

We initialise the disc with a circular cavity in the central region, such that the initial gas density is set to be

$$\begin{array}{cc}{\Sigma }_{\text{ini}}={\Sigma }_0{\left(\frac{r}{a_{\text{bin}}}\right)}^{-0.5}\cdot f_{\text{cav}},& f_{\text{cav}}\end{array}={\left[1+\exp \left(-\frac{R-2.5a_{\text{bin}}}{0.1a_{\text{bin}}}\right)\right]}^{-1},$$(9)

with ∑₀ = 357433 g cm⁻². The reference temperature profile corresponds to a constant h = 0.05 as shown in Eq. (3). For β models, the initial temperature follows Eq. (5)) directly such that the pressure is given by

$$P_{\text{ini}}=\frac{R{T}_{\text{eq}}{\Sigma }_{\text{ini}}}{\mu }.$$(10)

In all models, we use α = 10⁻³.

We model the disc with the GPU-capable version of Pluto (Thun et al. 2017; Mignone et al. 2007) in a 2D cylindrical grid with N_r × N_ϕ = 684 × 1168 cells, which corresponds to a resolution of ≈9 cells per scale height in both directions. The radial domain extends between 1 –10 a_bin with logarithmic spacing, thereby excluding the binary from the simulation domain for low eccentricities. However, this does not affect the structure of the disc in a significant way. We investigated this by lowering the inner radius in locally isothermal test runs and found that 1 a_bin offers accurate results while maintaining a reasonable timestep, in agreement with Tiede et al. (2021)) who showed that material within inside 1 a_bin is accreted and does not change the dynamics of the outer disc.

At the inner boundary we implement a strict outflow boundary condition for density and radial velocity (∂,∑ = 0, v_r,_in = −|v_r(r_min)| The azimuthal velocity there is set to ∂_rΩ = 0 and the pressure is reflected at the inner boundary. At the outermost 10% of the radial domain we damp Σ and P to their initial profiles and u_r to zero. We use a density floor of 10⁻⁷ Σ_ini. In the radiative models we include a pressure limit equivalent to a maximum temperature of 3000 K and a minimum of 3 K.

3. Results

In this section, we lay out the influence of cooling on the structure of the inner circumbinary disc with results from our numerical simulations. We first showcase the evolution of the circumbinary disc and the properties of its cavity with a locally isothermal model in Sect 3.1. We then look at models with different cooling prescriptions in Sect. 3.2 and later analyse the effect of the parametrised cooling timescale β on the disc shape by comparing to a model where β drops over time in Sect. 3.3. Finally, we demonstrate a simple model that can reproduce a temperature profile that is comparable to the radiative runs using β as a function of Σ in Sect. 3.4.

3.1. Example of disc evolution in the locally isothermal limit

Most of the observed circumbinary planets orbit a binary with a moderate binary eccentricity e_bin ≤ 0.21. Therefore) we choose simulations with a locally isothermal prescription, previously studied for Kepler binaries (see Thun et al. 2017; Thun & Kley 2018) with e_bin = 0;15 as a standart of comparison tor our models in the next section.

In Fig. 1, we Display the disc evolution in 2D snapshots at t = {0, 1000, 5000, 50 000} T_bin. The binary clears out an eccentric inner cavity in the disc that eventually stabilises to an ellipse with a semi-major axis of a_cav = 4.25 a_bin and an eccentricity e_cav = 0.26, as can also be seen in Fig. 2.

In all 2D plots, a dashed green bounded curve traces an ellipse that is fit to the edge of the cavity, which is defined as the region where Σ drops by a factor of 10% compared to the current azimuthal density maximum. This is done similarly to Thun et al. (2017), using the peak density location to compute the argument of the ellipse’s apocenter and constructing the ellipse by finding the cavity edges along the line connecting its apsides, using the centre of mass as the focal point.

In Fig. 2, the orbital parameters of the cavity oscillate around a mean value. This oscillation is caused by the precession of the disc around the binary. Especially for large e_bin, the cavity shape experiences large oscillations caused by the precessing motion of the CBD around the binary and the changing alignment between the argument of periastron of the disc and the binary orbits. At times of minimal cavity size and eccentricity, the disc and the binary orbits are in alignment. Thun & Kley (2018) have shown that the precession is directly linked to the cavity shape in the locally isothermal case by relating it to the theoretical precession of a massless third body. We discuss this in the context of our models in the next section.

Fig. 2 Time evolution of the semi-major axis and cavity eccentricity for the locally isothermal test case with ebin = 0.15. The dashed line represents the mean value around which the orbital parameters oscillate.

3.2. Comparison of cooling models

In the previous section, we analysed a sample circumbinary disc model as it evolves towards a quasi-steady state and laid out our method of measuring the cavity parameters in equilibrium with our locally isothermal example model. Here, we investigate a range of models with different e_bin end cooling prescriptions as listed in Table 1 using this method. Each model lias run for at least 50 000 T_bin, until the cavity profile has converged.

Figure 3 shows 2D snapshots of the disc surface density at a date at 50 000 T_bin. The cavity becomes small and circular for β = 1 independent of e_bin, and its size grows with binary eccentricity. The latter is in agreement with the instability region depending en binary eccentricity as seen in the early dynamical analysis by Artymowicz & Lubow f1994) or tire numerically evaluated instability region by Holman & Wiegert (1999), giving a minimum stable region for the disc or planets. This is a hint that the conditions for the β = 1 case lead to reduced hydrodynamical interaction inside the disc. Radiative and locally isothermal models produce larger and more eccentric cavities in all cases, as they correspond to an effective β ≈ 10³ and β ➞ 0 respectively. The high- and low-β models are accordingly comparable to the radiative or locally isothermal case with the exception of {e_bin = 0, β= 100}, where the disc remains circular, while the locally isothermal and radiative models show very eccentric cavities. To a lesser degree this also affects the highly eccentric disc caused by the eccentric binary in our model with {e_bin = 0.5, β = 100). The possible causes and solutions for problems with the β = 100 cases will be addressed in Sect. 3.4. The cavity size and eccentricity in radiative and isothermal models are smallest for e_bin = 0.15, in agreement; with Thun & Kley (2018). For very eccentric binaries (e_bin = 0.5) the cavity becomes very wide and eccentric, which also leads to wider planetary parking orbits around such systems as is the case for Kepler-34b (Welsh et al. 2012. Penzlin et al. 2021).

The time evolution of the cavity parameters a_cav and e_cav is shown in Fig. 4. The average values of the cavity shape over the last 10 000 T_bin are listed in Table 1 and can differ slightly from the values at; the time of the snapshot in Fig. 3. The curves in Fig. 4 show oscillations around these average values, an effect caused by the precession of the disc. As visible in the oscillations in Fig. 4, typical precession periods are in the order of ~10³ T_bin (see also Thun et al. 2017). The oscillation amplitude increases with higher e_bin, as the disc pericenter approaches the binary apocenter more closely once per precession.

The cavity in the radiative model with e_bin = 0.15 reaches its maximum size after 20 000 T_bin but continues evolving, shrinking to slightly lower final values of a_cav and e_cav over 50 000 T_bin as the disc surface density and temperature profiles readjust to an equilibrium state that corresponds to a smaller but self-consistently calculated aspect ratio h ≈ 0.04−0.05 in the inner disc (also shown in Kley et al. 2019). Given this slightly lower h, a smaller and less eccentric cavity is to be expected (Penzlin et al., in pres.).

Contrary to other models, where the cavity reaches an equilibrium state within ≈ 50k T_bin, the cavity of the radiative model with e_bin = 0 is still increasing in size even after > 100 000 T_bin. While it is not yet dully understood why this is the case, see speculate that it might be the result of a combination od a very long effective cooling timescale β ≈ 10³ (see Sect. 2.1.3) with the generally longer evolution timescale of models with e_bin = 0 due to the accumulation of long-lived eccentricity modes created by the symmetrically rotating binary that are trapped by the steeply truncated disc (Muñoz & Lithwick 2020).

The precession rate of the discs is of the order of a few 10007 T_bin and the exact number for the models with e_bin = 0.3 derived from the argument of periastron of the inner disc is displayed in Fig. 5 (coloured dots). When the density peaks are narrow as Fig. 3 shows, this precession can be approximated with a three-body orbit at the cavity edge tor locally isothermal models like in Thun & Kley et al 2018)

$$\begin{array}{cc}{T}_{\text{prec}}=\frac{3}{4}\frac{{\left(q+1\right)}^2}{q}{\left(\frac{a_{\text{cav}}}{a_{\text{bin}}}\right)}^{3.5}\frac{{\left(1-e_{\text{cav}}\right)}^2}{1+1.5e_{\text{bin}}^2}{T}_{\text{bin}},& q=M^2/M_1.\end{array}$$(11)

However, in discs with a wider density maximum like in the high β or radiative cases (Kley et al. 2019), the period of the precession increases as the disc precesses rigidly and the associated third-body orbit shafts outwards. Therefore, even for equal cavity sizes between radiative and locally isothermal models, the precession rate linked to the width and position of the peak density can still differ, and is usually higher for longer cooling timescales. This can be understood by comparing the precession raies of the emulations with those obtained from Eq. (11) for ellipses that are defined at locations in the cavity where the density is 10% and 50% of the peak density, as shown in Fig. 5. For locally isothermal models and those with β ≤ 1 the density peak is so narrow that both values match the precession well. However, for β = 00 the density peak is much wider so that the 3-body precession for the 10% case is significantly lower and does not match the precession rate anymore. This difference in peak position and precession rate is a result of the higher pressure in the β = 100 case near the inner edge, which influences the propagation of excitations in the disc similar to the effect studied in Teyssandier & Ogilvie (2016). The radiative model behaves differently, showing a minor increase in precession period. This is likely caused by the reduced pressure in the inner disc (Kley et al. 2019), leading to less rigid disc rotation that affects a smaller section of the disc.

Fig. 3 Comparison of the 2D surface density at steady state of all models in Table 1 after 50 000 Tbin. The surface density is normalised to the reference power-law profile (Σref ∝ r−05). The annotated values for acav (in abin) and ecav refer to the ellipse fit in dashed green lines.

Fig. 4 Time evolution of acav and ecav for all models presented in Table 1. The dashed vertical line marks the timestamp at t = 50 000 Tbin, which is used for comparison in Fig. 3. The data shown has been smoothed with a rolling mean over 300 Tbin for a clearer view.

Fig. 5 Precession period of the disc (coloured dots) for all models with ebin = 0.3. For comparison, the 3-body precession period for the 10% and 50% peak density ellipse are plotted with “+” and “x”. For the 3-body precession, the semi-major axis and eccentricity have been averaged over 4000 Tbin.

3.3. The influence of β on the cavity shape

As we have seen in Sect. 3.2, the cavity for β = 1 becomes small and nearly circular, while generally circumbinary disc cavities are eccentric even for large variations in aspect ratio and viscosity for locally isothermal (see e.g. Thun et al. 2017; Ragusa et al. 2020; Muñoz et al. 2020) as well as radiative discs (Kley et al. 2019). We ran models for e_bin = 0.15 with β = {10⁻¹, 10¹} in addition to the values in Table 1 and calculated the converged cavity shape for all five β models. The resulting V-shape, displayed by the points in Fig. 6, is centred around β = 1 and verifies that this value most likely corresponds to a total minimum.

Based on our results using different models with a constant β, we constructed a model of a disc where β varies in time from 10² to 10⁻² over 40 000 T_bin with

$$\log \beta =2-\frac{t}{{10}^4{T}_{\text{bin}}},$$(12)

using the equilibrium state of the e_bin = 0.15, β = 10² model at 50 000 T_bin as our initial state. The evolution of the cavity shape for this model is displayed in Fig. 6. The averaged curves show a clear minimum at β = 1 for both a_cav and e_cav, with the cavity shrinking and circularising as β approaches this value. Due to the initially long cooling timescale, the disc’s response to the cavity-shrinking effect as β approaches unity is significantly slowed, so th_et the orange curve does not exactly match our data points during its descent towards β = 1. As a result, e_cav is higher than the value predicted for our fixed-β model when β ≈ 10. For β ≲ 1, however, our variable-β model agrees very well with our models with a constant β.

The periodic variations in the blue curve in that figure are an artefact of changes during one precession period of the disc around the binary, as analysed in Thun et al. (2017), Thun & Kley (2018) and Kley et al. (2019) (see also Fig. 5). Therefore, even though the cavity sizes at both ends of the curve are comparable, the precession period varies due to the difference in the inner disc structure as described in Sect. 3.2.

Fig. 6 Cavity eccentricity and semi-major axis of the models with ebin = 0.15 and constant β values (black dots) and the model with time-varying β (blue curve) following Eq. (12). The orange curve is the rolling average of the blue model over 5000 Tbin.

Fig. 7 Heatmaps of β in radiative models. The cavity is marked by the green dashed line. We highlight the uniformly shaded disc, indicating thermal equilibrium between viscous heating and β-cooling, and the shock-dominated cavity where β deviates significantly from the expected value of 103. Regions with a very low β inside the cavity act as an indicator of spiral streamers, where heating is strongly dominated by shocks.

3.4. Accounting for the cavity in models with slow cooling

In the models presented, the actual pressure scale height can vary from the initially prescribed profile as the disc equilibrium temperature is changed by spiral shock heating or radiative cooling. This also affects models with β ≫ 1 significantly, as the slow cooling cannot counteract spiral shock heating efficiently near the binary. The effective β of radiatively cooled models, computed using Eq. (8), is displayed in Fig. 7. As described in Sect. 2.1.3, it corresponds to β ≈ 10³ inside the disc, but drops significantly within the cavity with the exception of spikes due to a combination of viscous and spiral shock heating around very low density streamers. “This indicates that viscous heating is not efficient enough within the cavity, consistent with previous works (Kley et al. 2019), and suggests that the cooling timescale should be much lower in this low-density region leading to unnaturally increased temperatures for the high-β models near the cavity.

As already seen in Sect. 3.2, some models with β = 100 display unexpectedly circular cavity structures. To understand the actual disc properties of the models which influence the dynamics, we plotted the aspect ratio h for the β = 100 model and the radiative model with e_bin = (tin Fig. 8. The aspect ratio here is derived from the disc with h = c_sr/Ω_K as shown in Sect. 2. Additional models that are described further below aim to reconcile the differences between the two main runs. Coloured dots mark the radial distance where the density is lower than t0% of the peak density, which is the criterion we use to define the cavity edge. High-β regions inside the cavity as seen in Fig. 7 contain so tittle material that they show little relevance on the density-weighted h. The radiative model has a nearly constant, low h within the disc reciting into the cavity. In radiative models, the aspect ratio spikes only close to the binary due to shocks created by the binary motion as also seen in Fig. 7. For that reason, density weighting was used to remove spikes tue to shock heating inside the low density region in Fig. 8. The azimuthal dependence of h is better highlighted in Fig. 9.

However, for the β = 100 model the aspect ratio increases substantially within the disc near the cavity edge. The very high aspect ratio (h > 0.1) of the β = 100 model at the inner edge leads to a stronger local viscosity $\left(v\propto h^2\sqrt{r}\right)$, which can suppress the excitation of the eccentric cavity by smearing out the wakes in the cavity caused by the binary motion while also refilling the forming cavity more efficiently. This hints at the temperature profile inside the cavity playing an important role in the final shape of said cavity.

In light of this deviation of the aspect ratio profiles, we implemented two additional models. First, a density-dependent prescription β(Σ) where β is constant within the disc but drops in low density regions, when the density drops below the threshold value Σ_th = 5% Σ₀. Then, for a given input β_ref, the density-dependent β(Σ) is given as follows:

$$\begin{array}{cc}\log \beta \left(\Sigma \right)=\log {\beta }_{\text{ref}}-\left(2+\log {\beta }_{\text{ref}}\right)\cdot \frac{{\Sigma }_{\text{th}}-\Sigma }{{\Sigma }_{\text{th}}},& \Sigma \le {\Sigma }_{\text{th}}\end{array},$$(13)

such that it drops smoothly towards 0.01 within the cavity. The resulting aspect ratio profile of this model is given by the orange curve in Fig. 8 and avoids the unrealistic long cooling time near that cavity edge that leads to the unphysical increase in temperature in the inner disc.

The second additional model uses instead an adaptive β(Σ, T) that; depends on the local thermodynamical properties of the gas following Eq. (8), while relaxing the temperature profile to one where h = 0.044 in order to mimic the radiative model more closely. To stabilise this method in the cavity, where β varies strongly with position, as well as to prevent runaway heating when β > 1/α (see Sect. 2.1.3), we lim it the value of β between 10⁻³ (locally isothermal limit)) and 10³ (adiabatic limit). This prescription, which attempts to emulate the results of the radiative model, is shown with red (turves in Figs. 8 and 10 and shows very good agreement with the radiative model in terms of cavity properties. This model also closely follows the envelopment of e_cav(t) in the radiative model as can be seen on Fig. 10, while the simpler density-dependent β(Σ) models both lag behind the radiative evolution. Nevertheless, e_cav in all adaptive β models follows an increasing trend and is closer to the radiative model than the locally isothermal model discussed in Sect. 3.2.

The heatmap in Fig. 9 also shows that; the structure and size of the inner cavity of the improved models are comparable, while the dimple β = 100 model deviate, significantly. The plotted maps of h(r, 0) illustrate this behaviour, with the β = 100 model having a substantially higher aspect ratio in the cavity. The uniformly high pressure in the cavity for β = 100 suppresses the wakes launched by the binaries that are otherwise visible ins all other models by stream s of high aspect ratio caused by shocks inside the cavity.

However, the adjusted β(Σ,T) model requires a priori knowledge of the thermal profile of the disc, while at the same time results in a strictly hotter disc in legions where shock heating is negligible due to the very long cooling timescale in the bulk of the disc (see Fig. 7 and Eq. (5t). This excess heating becomes evident lay looking at the aspect ratio profile of this model in Figs. 8 and 9, which show that this adaptive model behaves similarly to the radiative model outside of the cavity but its h is consistently higher in value.

Fig. 8 Comparison of azimuthally averaged, density-weighted aspect ratio profiles of the fiducial model with β = 100 and h = 0.05 (blue) against, density-dependent. β(Σ) for h = 0.05 (orange) and h = 0.044 (green), adaptive β(Σ, T) according to Eq. (8) (red), and radiative models (purple) after 50 000 Tbin for = ebin.= 0 Horizontal dashed lines denote the equilibrium state for β models where shock f eating would be absent Coloured dote mark. the respective acav.

Fig. 9 Comparisons of surface density (top) and aspect ratio (bottom) heatmaps between (from left to right) β = 100, density-dependent ß(Σ) (h = 0.05 and h = 0.044), adaptive β(Σ, T) and radiative models for ebn = 0, at time 50 000 Tbin. The top panels have been normalised to their reference power-law profiles Σref- ∝ r−0.5.

Fig. 10 Time evolution of ecav for the models shown in Fig. 8. The data shown has been smoothed with a rolling mean over 300 Tbin for clarity.

4. Discussion

In this section, we consider the influence of different physical effects as well as binary parameters on the topic of the applicability of our results on observed systems. We also highlight some caveats of our models and discuss possible improvements in future modelling, with particular focus on larger circumbinary discs.

4.1. Effective range of thermal models

In this study we focus on configurations that represent the circumbinary planet systems observed by the Kepler mission, and therefore we chose close binary configurations with a_bin = 0.2 au. In this regime, viscous heating is the dominant heat source as also shown in Kley et al. (2019), while stellar irradiation has a negligibly small contribution. Assuming thermal equilibrium, the cooling time inside the disc in this regime is then β ≈ 1/α, as shown in Sect. 3.4. Therefore, the disc of these close-in systems can be reasonably well-represented by density-dependent β(Σ) models with β_ref ≲ 10² if one takes into account that the cooling timescale is much shorter within the cavity and that the aspect ratio in radiative models can deviate towards slightly lower values in thermal equilibrium.

For very large discs such as GG Tau with a cavity size of a_cav ~ 180 au, irradiation heating (e.g. Chiang & Goldreich 1997; Menou & Goodman 2004) would become more relevant compared to viscous heating, which would result in lower temperatures and conditions closer to locally isothermal as the disc becomes optically thinner. While β models are easily scalable, viscously heated models only represent the circumbinary disc around close binaries. To investigate realistic thermal conditions around observed large CBDs, a model should account for stellar irradiation in addition to viscous heating. To adapt a β model to such a more realistic case, β would need to depend on radius or even azimuth, as well as drop in optically thin regions such as inside the cavity or at large radii, far from the binary. In that sense, a lower β corresponds to the case of an even wider disc. As a result, the different cases introduced in this work bear relevance to wide circumbinary discs, especially ones where locally isothermal conditions are likely. Models including irradiative terms would be the topic of future work.

However, as we have demonstrated, the main difference between our radiative and locally isothermal models even for short scales is the precession rate, while conditions close to a β = 1 case cause the disc to radically change its cavity structure. This effect may be exacerbated in this work as our models do not include the effects of in-plane cooling via flux-limited diffusion (Miranda & Rafikov 2020b), which could possibly mediate the interaction between cavity and disc in this case. This could be investigated further in future work.

4.2. Planet formation scenario

The models presented here are created with the close-in circumbinary planets in mind. We cover the full range of binary eccentricities relevant to the planet case from circular binaries like Kepler-47 (Orosz et al. 2019) and Kepler-413 (Kostov et al. 2014) to the most eccentric case of Kepler-34 (Welsh et al. 2012). Such discs are a representation of environments of planet formation around binaries. However, given the masses of circumbinary planets of ~0.1 M_Jup, an inserted planet will change the dynamics of the disc (Penzlin et al. 2021). Especially when it reaches an orbit inside the cavity, the planet will shield the disc from further excitation by the binary, leading to a circularisation that may appear similar to the β = 1 case, albeit caused by different interactions.

4.3. On the binary mass ratio

Our simulations all use a mass ratio of M₁/M₂ = 0.5. Using locally isothermal models, Thun & Kley (2018) have shown that variations of the mass ratio only have a minor effect on the cavity shape. D'Orazio et al. (2016) found a behaviour transition in discs around binaries with mass ratios above M₁/M₂ > 0.04 that leads to asymmetric, time-variable discs with a large cavity. As a result, for stars of comparable size the large, eccentric cavities explored in this work are expected. We chose an intermediate mass ratio of 1:2 such that our general results can be expected to be comparable to binaries within a range of mass ratios. The exact shape and properties of the cavity as a function of the binary mass ratio is beyond the scope of this study.

4.4. Three-dimensional effects

We use a 2D cylindrical grid. As our simulations require up to 10⁵ T_bin for the cavity to converge to a steady state, running 3D models is simply not feasible. Shi et al. (2012) performed 3D MHD simulations of a circumbinary disc for up to 1000 T_bin, and found an increasing trend in disc eccentricity. Pierens et al. (2020) showed the same trend in their 3D simulations of CBDs, where comparably eccentric cavities develop. Using 3D smooth particle hydrodynamics, Ragusa et al. (2020) also reproduced the growth of eccentric cavities for locally isothermal models. Hence, neglecting the vertical direction in our simulations does not appear to affect the general picture of an eccentric circumbinary disc with a large central cavity. One drawback of 3D simulations is that the computational cost does not allow for the long simulation times necessary to reach a converged state with the precessing disc, as can be seen in our 2D models. As the cavity precesses, its properties vary as shown in Fig. 2. As precession times can be in the order of 10³ T_bin it is difficult to assess the convergence between 2D models that run for 50 000 T_bin and 3D models that run for only 1000–3000 T_bin The eccentricity reached at the end of the aforementioned 3D simulations is compatible with that of our 2D simulations at the same time. Nevertheless, the exact structure and properties of the cavity could be different in a fully 3D model in equilibrium. Such a model is beyond the scope of this study.

4.5. Computational domain

In this study, we chose an inner domain boundary at 1 a_bin such that the full interaction of the binary is not included for runtime reasons. Tiede et al. (2021) shows that inside this radius nearly all mass is accreted and does not alter the dynamics of the outer disc. As Fig. 7 also reveals, the mass of the streamers inside the cavity is 3 to 4 orders of magnitude smaller than the disc density at the converged state for the chosen viscosity. In very viscous or thick systems for example, advection can become the more relevant term for angular momentum transport and can play an important role in the dynamics of such discs (e.g. Shi et al. 2012; Muñoz et al. 2019). In these protoplanetary disc cases the gravitational torque drives the eccentricity of the disc as described in Muñoz & Lithwick (2020). At the same time, the effect of the disc on the binary motion is not taken into account in our model. For the chosen viscosity the orbital change of the binary is in the order of 10⁻⁷ a_bin/T_bin (Penzlin et al. 2022) and thereby much longer than the relevant timescales of the simulation.

4.6. Different cooling prescriptions

As mentioned in Sect. 3.4, we implemented additional models with β-cooling in an effort to construct a prescription that can reproduce the results of our radiative runs while also gaining the speed advantage of β methods, which do not involve using a stand-alone cooling module in comparison. In particular, we compared the radiative model against setups with a density-dependent β(Σ) designed to quickly drop inside the cavity using a simple exponential cut-off, as well as a setup where β is computed self-consistently through Eq. (8) by taking into account the local density, temperature, and opacity of each cell.

We found that both adaptive models are approximately 30% faster compared to the radiative model, maintaining the same timestep but requiring significantly less computational time per step. The β(Σ) model was marginally faster than the β(Σ,T) model by 3–5%, as it did not require the use of an opacity function, which case can be quite computationally expensive.

We note that the β(Σ,T) model resulted in temperatures noticeably higher than in the radiative case. Nevertheless, β(Σ) model with β_ref ➞ 10³ instead of the used 10² would result in a hotter disc as well (see Eq. (5)). Hence, a β(Σ) model depends on an appropriate choice for β_ref which introduces additional uncertainty. For a larger-scale system where irradiation dominates over viscous heating such that β would be of the order of 10⁻¹−10², the β(Σ,T) model would approach the correct radiative profile more accurately while maintaining its significant speed-up factor over it. Nevertheless, β(Σ) model would most likely also provide comparable results if a radius-dependent β_ref is used. Of course, in both cases, a radiative model would provide the best results at the cost of a moderate computational overhead.

Based on our results, we stress that the choice of the target aspect ratio h0 is important as it will dictate the aspect ratio in equilibrium heq through Eq. (5). It might be possible to reach a better agreement between β and radiative models by intentionally choosing a lower h₀ such that h_eq ≈ 0.044. Nevertheless, our adaptive models match the radiative models in terms of cavity shape reasonably well even with our prescription of h₀.

5. Conclusions

In this study, we have investigated three thermodynamical models in the context of circumbinary discs using 2D numerical simulations of locally isothermal, β-cooled, and radiative models. The circumbinary environment naturally introduces a large inner cavity, the size and shape of which depends on the specific prescription of disc and binary parameters.

We find that the cavity itself develops to comparable sizes and shapes between the locally isothermal and radiative models. This happens because in the very low-density, optically thin cavity the cooling timescale of the radiative model is short enough to be comparable to a locally isothermal cooling prescription. However, the structure inside the inner disc and the width of the high-density region at the cavity edge differs between locally isothermal and radiatively cooled models. This divergence in the inner disc is also reflected in the different precession rates between the locally isothermal and radiative models.

By using β models that simulate conditions of a nearly locally isothermal disc (β ≪ 1), we reproduce results that match the locally isothermal model well. However, when naively using a high value of β = 100 to match the radiatively cooled prescription, models can deviate significantly. The reason is an unphysical treatment of the cooling timescale inside the cavity, leading to overheating of the inner disc. Such excess heat leads to an overly strong pressure increase in the inner disc, pushing the disc farther towards the binary. This leads to circularisation of the cavity especially for e_bin = 0 as the wakes caused by the binary motion are no longer able to overcome the pressure-supported inner disc structure.

We introduce two possible solutions of adaptive β prescriptions to achieve parametrised models that can potentially be more physically accurate than constant β models: an empirical, opacity-model-agnostic prescription in which densities below a threshold value lead to a drop in the cooling timescale such that the effective β(Σ) profile emulates both the optically thin cavity and optically thick disc as would be expected in a radiative model, and a more physically motivated model where β(Σ,T) is evaluated locally in a self-consistent manner by computing the thermal cooling timescale based on the current density, temperature, and opacity of each cell. We found that the first model produces reasonably good results in terms of both temperature profiles and cavity properties, while the second model matches the cavity shape of our radiative runs more closely at the cost of a slightly hotter disc far from the binary. The latter happens due to the fundamentally different approach of a relaxation prescription using β – where the background temperature profile is assumed to be known in advance – as opposed to a truly radiative setup where the gas could, in principle, cool down to very low temperatures. The differences between our prescriptions of adaptive β and radiative models underline the importance of correctly modelling the binary cavity when comparing to observations.

After investigating models with different, constant cooling timescales β ϵ [10⁻²,10²] we find that the size and eccentricity of the cavity becomes minimal for values of β close to 1. This is likely due to the same effect observed in the work by Miranda & Rafikov (2020b), where radiative cooling of the spiral arms launched by the binary interferes with the angular momentum transport capabilities of said spirals for β ~ 1, thereby resulting in smaller and more circular cavities. We note that our models did not include the effects of in-plane cooling, which might affect our results according to their study.

Acknowledgements.

We dedicate this paper to our late mentor and friend Willy Kley. Without his support and guidance, this work would not have been possible. We will always be grateful for all the time he invested in teaching us and helping us grow. AP was funded by grant 285676328 of the German Research Foundation (DFG). The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC, and the German Research Foundation (DFG) through grant no. INST 37/935-1 FUGG. AZ and RPN are supported by STFC grant ST/P000592/1, and RPN is supported by the Leverhulme Trust through grant RPG-2018-418. All plots in this paper were made using the Python library matplotlib (Hunter 2007).

Appendix A Highly viscous environments

Accretion discs around black holes are believed to be much more viscous than typical PPDs. For a first idea on the effect of cooling in such systems, we tested out a small series of models using the parameters given in Sec. 2.2 with an e_bin = 0 15 and a high Shakura–Sunyaev parameter of α = 0.1. For this setup, we ran models with locally isothermal, radiative, and β = [1,9] prescriptions. Due to the shorter viscous timescale for this choice of α, simulations already show a fully converged state at t = 700 T_bin. Figure A.1 shows the density and thermal structure of the inner disc.

For highly viscous disc models, we see that the β prescription is ill-suited to simulate the inner cavity structure. There are minor differences, between the locally isothermal and radiative model with the locally isothermal model opening the widest eccentric cavity. In contrast, the radiative model shows a slightly smaller and diffused cavity due to the heating induced by the high viscosity. However, the β = 1 model shows a more shallow and almost completely symmetric cavity that does not compare well with the structure of the locally isothermal and radiative models. Moreover, the aspect ratio of the material and the streams inside the cavity increases to unrealistic levels. The β = 9 case does not reproduce the locally isothermal or radiative model appropriately. Cooling in the inner region is slow by construction such that the cavity growth is almost completely inhibited. As a result, the overall temperatures in the disc rise to extreme levels (see Eq. (5)), puffing up the disc to aspect ratios of nearly 0.2, which cannot describe the bulk of the disc realistically.

For such extremely viscous environments, treating radiative effects can produce different results depending on the exact physical parameters compared to the often-used locally isothermal prescription. However, a β-cooling model should not be used for such systems.

Fig. A.1. Comparisons of surface density (top) and aspect ratio (bottom) heatmaps for a viscous α = 0.1 between (from left to right) isothermal, β = 1, β = 9 and radiative models for ebin = 0.15, at time 700 Tbin. The top panels have been normalised to their reference power-law profiles Σref ∝ r−0.5.

References

All Tables

All Figures

Fig. 1 Time evolution of the 2D surface density normalised to its reference power-law profiles ∑re∝ r−0.5 for the locally isothermal test case with ebin = 0.15. The system eventually reaches an equilibrium by t = 50 000 Tbin. The dashed green curves represent the ellipse fitting. The binary spiral streams, which are present inside the cavity, are not visible within this colour range, but can be seen in Fig. 7.

In the text

	Fig. 2 Time evolution of the semi-major axis and cavity eccentricity for the locally isothermal test case with ebin = 0.15. The dashed line represents the mean value around which the orbital parameters oscillate.
In the text

	Fig. 3 Comparison of the 2D surface density at steady state of all models in Table 1 after 50 000 Tbin. The surface density is normalised to the reference power-law profile (Σref ∝ r−05). The annotated values for acav (in abin) and ecav refer to the ellipse fit in dashed green lines.
In the text

	Fig. 4 Time evolution of acav and ecav for all models presented in Table 1. The dashed vertical line marks the timestamp at t = 50 000 Tbin, which is used for comparison in Fig. 3. The data shown has been smoothed with a rolling mean over 300 Tbin for a clearer view.
In the text

	Fig. 5 Precession period of the disc (coloured dots) for all models with ebin = 0.3. For comparison, the 3-body precession period for the 10% and 50% peak density ellipse are plotted with “+” and “x”. For the 3-body precession, the semi-major axis and eccentricity have been averaged over 4000 Tbin.
In the text

	Fig. 6 Cavity eccentricity and semi-major axis of the models with ebin = 0.15 and constant β values (black dots) and the model with time-varying β (blue curve) following Eq. (12). The orange curve is the rolling average of the blue model over 5000 Tbin.
In the text

Fig. 7 Heatmaps of β in radiative models. The cavity is marked by the green dashed line. We highlight the uniformly shaded disc, indicating thermal equilibrium between viscous heating and β-cooling, and the shock-dominated cavity where β deviates significantly from the expected value of 103. Regions with a very low β inside the cavity act as an indicator of spiral streamers, where heating is strongly dominated by shocks.

In the text

Fig. 8 Comparison of azimuthally averaged, density-weighted aspect ratio profiles of the fiducial model with β = 100 and h = 0.05 (blue) against, density-dependent. β(Σ) for h = 0.05 (orange) and h = 0.044 (green), adaptive β(Σ, T) according to Eq. (8) (red), and radiative models (purple) after 50 000 Tbin for = ebin.= 0 Horizontal dashed lines denote the equilibrium state for β models where shock f eating would be absent Coloured dote mark. the respective acav.

In the text

	Fig. 9 Comparisons of surface density (top) and aspect ratio (bottom) heatmaps between (from left to right) β = 100, density-dependent ß(Σ) (h = 0.05 and h = 0.044), adaptive β(Σ, T) and radiative models for ebn = 0, at time 50 000 Tbin. The top panels have been normalised to their reference power-law profiles Σref- ∝ r−0.5.
In the text

	Fig. 10 Time evolution of ecav for the models shown in Fig. 8. The data shown has been smoothed with a rolling mean over 300 Tbin for clarity.
In the text

	Fig. A.1. Comparisons of surface density (top) and aspect ratio (bottom) heatmaps for a viscous α = 0.1 between (from left to right) isothermal, β = 1, β = 9 and radiative models for ebin = 0.15, at time 700 Tbin. The top panels have been normalised to their reference power-law profiles Σref ∝ r−0.5.
In the text