© The Authors 2026

1. Introduction

In recent years, the interaction between massive black hole binaries (MBHBs) and circumbinary discs has been extensively investigated with the aim of better understanding how the disc affects the evolution of the binary orbit. Despite significant progress, the diversity in numerical methods (e.g. 2D versus 3D codes or Eulerian versus Lagrangian codes) and adopted physical models (e.g. locally isothermal versus adiabatic equations of state, inclusion of the disc self-gravity) has resulted in a variety of results, making it difficult to piece together a coherent picture (Muñoz et al. 2019; Duffell et al. 2020; Tiede et al. 2020; Heath & Nixon 2020; Shi et al. 2012; Moody et al. 2019; Franchini et al. 2022). Within this context, Duffell et al. (2024) presented results obtained using a common, simple binary-disc setup and performing numerical simulations with a variety of hydrodynamic codes to measure, among other properties, the magnitude of the gravitational torque that the disc exerts on the binary. While the study shows general agreement on the sign of the gravitational torque between different codes, there are still significant differences in its absolute magnitude and in other aspects of the binary-disc evolution that essentially depend on the code properties (i.e. 2D versus 3D, Eulerian versus Lagrangian) and on the employed boundary conditions.

The choice of a physical model for the disc greatly influences the outcome of the interaction with the binary. For instance, the disc aspect ratio, H/R, was recently found to play a critical role in the binary orbital dynamics, determining whether the binary inspirals or outspirals as a result of the interaction with the gas (Tiede et al. 2020; Heath & Nixon 2020; Franchini et al. 2022). Additionally, self-gravity in circumbinary discs has been shown to regulate both the torque exerted onto the binary and disc temperature, leading the binary to shrink regardless of the initial temperature of the disc (i.e. initial aspect ratio; Cuadra et al. 2009; Roedig et al. 2012; Roedig & Sesana 2014; Franchini et al. 2021).

The majority of previous works used a locally isothermal equation of state (Farris et al. 2014; Zrake et al. 2021; D’Orazio & Duffell 2021; Franchini et al. 2022) and assumed a constant gas temperature profile through the disc over time, therefore limiting the ability to capture shock-induced heating and its effect on the disc and ultimately on the binary. For instance, some of the gas that leaks inside the cavity is flung back towards the inner edge of the cavity, producing shocks that can alter the disc aspect ratio and temperature profile (Artymowicz & Lubow 1996; Hayasaki et al. 2007).

In our previous work (Cocchiararo et al. 2024), we made a significant decision to include the gas self-gravity in the circumbinary disc together with a black-body-like cooling prescription for the gas thermodynamics evolution in 3D. We employed a live binary (Franchini et al. 2023) whose orbital parameters evolve over time under the influence of gravitational and accretion torques exerted by the disc. While we focused on the characteristics of the electromagnetic emission from the disc, it is worth mentioning here that we found the interaction between the binary and the circumbinary disc to cause the binary semi-major axis to decrease with time as the disc self-regulates and ultimately reaches a stable aspect ratio value in the inner parts (Franchini et al. 2021). We find that the initially circular binaries tend towards higher eccentricity values, in agreement with previous works (Roedig et al. 2011; D’Orazio & Duffell 2021; Siwek et al. 2022; Franchini et al. 2024b), whereas very eccentric binaries experience a negligible eccentricity evolution within the time frame of our simulations.

A critical aspect that has often been neglected in previous studies is the role of radiation pressure in the hydrodynamics evolution of the circumbinary disc. Previous works have already performed radiation-magnetohydrodynamics (RMHD) simulations (Sadowski 2016; Mishra et al. 2016) in both their Newtonian (Jiang et al. 2019; Huang et al. 2023) and general relativity (GR) versions (Zhang et al. 2025) around a single black hole as well as in the study of tidal disruption events (Shiokawa et al. 2015; Ryu et al. 2023). At sub-parsec scales, discs are likely dominated by radiation pressure; therefore, including this additional term in the hydrodynamics equations is very important to advancing the theoretical modelling of the MBHB-disc interaction. Radiation pressure plays an important role in the hotter inner regions of the disc, where it significantly affects the gas dynamics, potentially modifying its geometry (including its aspect ratio, cavity shape, and eccentricity) and therefore altering the inflow of gas towards the binary. Correctly modelling the inflow of gas inside the cavity is fundamental, as this affects the gravitational and accretion torques on the binary, the accretion rate, and therefore the evolution of the binary orbital parameters, i.e. semi-major axis and eccentricity.

In this work, we evaluate the impact of radiation pressure on the evolution of circumbinary discs around MBHBs for different values of the binary initial mass ratio and eccentricity. Building on our previous simulations (Cocchiararo et al. 2024), we implemented the contribution of the radiation pressure in our 3D numerical simulations with hyper-Lagrangian refinement and considered milli-parsec scale binaries. The simulations also account for the thermodynamics evolution of the gas using a radiative cooling prescription in the form of black-body radiation. We have also included the self-gravity of the disc and the Shakura-Sunyaev prescription for viscosity (Shakura & Sunyaev 1973).

This work is the first part of a larger project focusing on the role of radiation pressure in the evolution and emission of MBHBs at milli-parsec separation. Here, we investigate its impact on the binary and circumbinary disc dynamics. We have dedicated a follow-up paper (Cocchiararo et al. 2025) to the radiation pressure influence on the electromagnetic emission and the observational signatures of these systems.

We explored the time evolution of the binary and disc properties for three values of the binary eccentricity, e = 0, 0.45, 0.9, and mass ratio, q = 0.7, 1. We compare the evolution of the disc and binary orbital parameters with the simulations that included only gas pressure (Cocchiararo et al. 2024) in order to isolate the effect of radiation pressure.

The paper is organised as follows. In Sect. 2, we describe the numerical details of the simulations, how we implemented the radiation contribution to the total pressure, the physical parameters we used, and the circumbinary disc assumptions. We show and discuss the main results we obtained, including the time evolution of the main binary and disc properties, in Sect. 3. Finally, we report our conclusions in Sect. 4.

2. Numerical and physical model

We used the 3D meshless finite mass (MFM) version of the code GIZMO (Hopkins 2015) to model the interaction between the MBHB and the circumbinary disc. We employed an adaptive hyper-Lagrangian refinement, i.e. particle splitting, to achieve a higher resolution inside the cavity carved by the binary (see Franchini et al. 2022, for details). We used the same refinement scheme as in our previous paper (Cocchiararo et al. 2024). In code units, the total mass of the binary is M_B = M₁ + M₂ = 1, where M₁ and M₂ are the masses of the primary and secondary black hole, respectively, and the initial semi-major axis is a₀ = 1. The two binary components were modelled as sink particles with an accretion radius of R_sink = 0.05a₀. As the binary orbit evolves with time (Franchini et al. 2023), the mass, linear, and angular momentum are conserved during each accretion event, following the method used in the PHANTOM code (Bate et al. 1995). We sampled the gas in the circumbinary disc with N = 10⁶ particles for a total disc mass of M_D = 0.01M_B. This mass is distributed with an initial surface density profile Σ ∝ R⁻¹ between R_in = 2a and R_out = 10a in the circular cases and between R_in = 3a and R_out = 10a in the eccentric cases. The disc is co-planar with the binary orbit and has an initial aspect ratio of H/R = 0.1. We generated the initial conditions of the disc using the SPH code PHANTOM (Price et al. 2018).

The thermodynamics evolution of the gas was governed by an adiabatic equation of state with index γ = 5/3. The gas was therefore allowed to heat and cool so that we could capture the effect of shocks. We modelled the angular momentum transport within the disc using the α viscosity parameter, which constitutes a simple parametrisation of the turbulence within the disc (Shakura & Sunyaev 1973). We included viscosity through the Shakura-Sunyaev parametrisation with kinematic viscosity ν = αc_sH, where c_s is the gas sound speed, α = 0.1, and the disc thickness (H) is H = c_s, i/Ω_K. We also included the gas self-gravity (Lodato 2007) and a black-body-like cooling. We started our discs in gravitational equilibrium by setting the initial Toomre parameter to Q > 1 (Toomre 1964).

We assumed radiative cooling in the form of black-body radiation. The cooling rate is given by

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{cool}}=\frac{8}{3}\frac{{\sigma }_{\mathrm{SB}}{T}^4}{\kappa \mathrm{\Sigma }},\end{array}$$(1)

where σ_SB is the Stephan-Boltzmann constant, κ is the opacity, Σ is the disc 3D surface density¹, and T is the temperature. We assumed that the opacity, κ, is a combination of Kramer opacity, κ_Kramer ∝ ρT^−7/2, and the electron scattering opacity, κ_es = 0.2(1 + X) cm² g⁻¹, with a hydrogen mass fraction of X = 0.59. We do not model optically thin regions; we assumed optical thickness larger than one and neglected any effect related to a possible atmosphere around the disc. Since we introduced a black-body cooling, we needed to choose meaningful physical units for our simulations. The total binary mass is M_B = 10⁶ M_⊙ and the initial separation is a₀ = 4.8 × 10⁻⁴ pc ≃ 1.2 × 10⁴ R_g, where R_g = GM_B/c² represents the gravitational radius of the binary.

We implemented the radiation pressure contribution in GIZMO using an approximate formulation that does not require solving the full radiation hydrodynamics equations. We first determined the initial temperature, T_0, i, by solving the following relation, which ensures that the initial pressure profile is consistent with hydrostatic equilibrium in the vertical direction:

$$\begin{array}{c}\hfill \frac{4}{3}\frac{{\sigma }_{\mathrm{SB}}}{c}{T}_{0,\mathrm{i}}^4+\frac{{\rho }_{\mathrm{i}}{k}_{\mathrm{B}}}{\mu m_p}{T}_{0,\mathrm{i}}-{\rho }_{\mathrm{i}}{c}_{\mathrm{s},\mathrm{i}}^2=0,\end{array}$$(2)

where c_s, i is the sound speed of the gas element i determined by the radial profile

$$\begin{array}{c}\hfill c_{\mathrm{s},\mathrm{i}}=\frac{H}{R}\sqrt{\frac{M_1+M_2}{a_0}}{r}^{-0.5}.\end{array}$$(3)

Here, H/R is the disc aspect ratio, M₁ and M₂ are respectively the masses of the primary and secondary black holes, a₀ is the initial semi-major axis, and r is the radius in units of a₀.

Using T_0, i, we calculated the initial internal energy, u_0, i, for each gas element, i, as the sum of the gas internal energy and the radiation internal energy. The values of u_0, i are included in the initial conditions file, ensuring that the contribution of the radiation pressure is accounted for at the beginning of each simulation. We then evolved the internal energy equation over time.

At each timestep of the simulation, the temperature, T_i, was re-calculated from the updated internal energy, u_i, by solving

$$\begin{array}{c}\hfill 4\frac{{\sigma }_{\mathrm{SB}}}{c{\rho }_{\mathrm{i}}}{T}_{\mathrm{i}}^4+\frac{3}{2}\frac{k_{\mathrm{B}}}{\mu m_{\mathrm{p}}}{T}_{\mathrm{i}}-u_{\mathrm{i}}=0,\end{array}$$(4)

and the ideal gas pressure was computed as P_gas = (ρ_ik_BT_i)/(μm_p). To account for the contribution of radiation pressure, we computed the parameter β, defined as the ratio of gas pressure to total pressure:

$$\begin{array}{c}\hfill \beta =\frac{1}{1+\frac{4}{3}\frac{{\sigma }_{\mathrm{SB}}}{c}\frac{T^4}{P_{\mathrm{gas}}}}.\end{array}$$(5)

We used this parameter to modify and update the effective adiabatic index, γ_β, in the equation of state, ensuring a self-consistent treatment of the coupled gas and radiation pressure during the evolution using the following definition:

$$\begin{array}{c}\hfill {\gamma }_{\beta }=\beta +{(4-3\beta )}^2\frac{\gamma -1}{\beta +12(\gamma -1)(1-\beta )}.\end{array}$$(6)

The adopted approach has some limitations, as it does not treat optically thin regions. However, it is computationally efficient and captures the main impact of radiation pressure in optically thick gas, which is the environment we are interested in regarding the context of this work. We ran six different simulations, three of which included the implementation of radiation pressure with the following combinations of binary eccentricity and mass ratio: e = 0 and q = 1, e = 0.45 and q = 0.7, e = 0.9 and q = 1. We ran the equal mass binary simulations for 1300 binary orbits and the unequal mass simulation for 1600 binary orbits to ensure that the disc reaches a quasi-steady state, i.e. the ratio between the change in binary angular momentum over accretion rate does not change with time.

3. Results

In this section, we present the results that we obtained from our numerical simulations neglecting and including the radiation pressure contribution, as explained in Sect. 2. In particular, we compare the evolution of the disc properties, such as the time evolution of the disc surface density, temperature, and aspect ratio. We then compare the evolution of the binary orbital parameters and accretion rate. Finally, we focus on the interaction between disc morphology, the alternation of the accretion rate onto each black hole, and the preferential accretion in binary systems.

3.1. Disc evolution

For a more comprehensive physical understanding of circumbinary disc evolution, it is essential to account for radiation pressure. While it is often included in GRMHD simulations (Sadowski 2016; Mishra et al. 2016; Cattorini & Giacomazzo 2024; Tiwari et al. 2025), radiation pressure is usually neglected in studies of binaries at large separations. Its effect is particularly important in the hot inner regions of the disc, where it can significantly alter the gas dynamics and the disc geometry and ultimately influence the inflow of material towards the binary. Figure 1 shows the surface density maps in the x − y plane for the simulations with e = 0 and q = 1, e = 0.45 and q = 0.7, and e = 0.9 and q = 1 respectively at t = 1,  900,  1300 P_B for the unequal mass binary and t = 1,  900,  1600 P_B for the equal mass binaries. In order to highlight the radiation pressure effect, we show in each of the three panels the same simulation that includes only the gas pressure contribution in the upper row.

Fig. 1. Surface density map in code units (1 code unit = 200 AU) in the x − y plane for the simulations with e = 0, q = 1 (upper panel); e = 0.45, q = 0.7 (middle panel); and e = 0.9, q = 1 (bottom panel) at three different moments: t = 1,  900,  1300 PB for the equal mass cases and t = 1,  900,  1600 PB for the unequal mass case. In each panel, the first and the second row show results for Ptot = Pgas and Ptot = Pgas + Prad, respectively. In the eccentric simulations, the initial cavity radius is 3a0, while in the circular binary case, it is 2a0. When the radiation pressure is included, the binary begins accreting at the very early stages of its evolution; otherwise, the disc experiences a transient phase lasting approximately 900 PB, 700 PB, and 600 PB respectively for the cases e = 0, q = 1; e = 0.45, q = 0.7; and e = 0.9, q = 1. The spatial resolution at the edge of the cavity (r = 3a) is Δx = 0.025a.

The first noticeable effect of radiation pressure is to inhibit the cavity eccentricity growth, particularly in the circular equal mass binary case. The simulation with the circular equal mass binary and the implementation of radiation pressure shows a region of enhanced density in the inner part of the disc that extends up to R ∼ 5a. This feature also persists after t = 1300 P_B. In the eccentric binary cases, the presence of radiation pressure results in a denser over-density at the cavity inner edge. This feature, called a ‘lump’, has previously been found in the literature in circumbinary discs around close-to-equal-mass binaries, and it is the result of the combination of the flung-back streams impacting on cavity edge and of the disc eccentricity (Shi et al. 2012; Farris et al. 2014; Noble et al. 2021; Franchini et al. 2024a). In the radiation pressure case, the dense region that forms in the equal mass circular binary case possibly hinders the formation of the lump.

With our 3D simulations that include time-dependent thermodynamics, we were able to resolve shocks arising from gas streams that are flung back into the cavity inner edge. These shocks do indeed affect the inner disc thickness and temperature. The radiation pressure contribution can also play a role in this process, potentially influencing the disc structure.

Figure 2 shows the time evolution of the surface density, disc aspect ratio (H/R)², the midplane temperature, and the effective temperature (T_eff) as a function of the radius (R) for all of our simulations at four different times. In all the simulations, for each gas particle i, we obtained the temperature T_i by numerically solving the following implicit equation for T_i:

Fig. 2. Time evolution of the surface density (first row), disc aspect ratio (H/R; second row), the midplane temperature (third row), and the effective temperature (last row) as a function of radius for the binary case e = 0, q = 1 (first column); e = 0.45, q = 0.7 (central column); and e = 0.9, q = 1 (last column) at different times (distinguished by different colours). The solid and dashed lines refer to the simulation without and with the implementation of the radiation pressure, respectively. The initial condition of the disc is the same across all the cases except for the initial cavity radius, which is 3a0 in both of the eccentric simulations, while it is 2a0 in the circular binary cases. We calculated the midplane temperature and the effective temperature assuming that both the gas and the radiation pressure contribute to the hydrostatic equilibrium of the disc in all the cases. The radiation pressure contribution leads the disc to maintain a lower aspect ratio and temperature.

$$\begin{array}{c}\hfill P_{\mathrm{tot}}=P_{\mathrm{gas}}+P_{\mathrm{rad}}=\frac{\rho k_{\mathrm{B}}{T}_{\mathrm{i}}}{m_{\mathrm{p}}\mu }+\frac{4}{3}\frac{{\sigma }_{\mathrm{SB}}{T}_{\mathrm{i}}^4}{c}.\end{array}$$(7)

We divided the disc temperature domain into a 2D matrix in the x-y plane, corresponding to the binary orbital plane. For each pixel, the midplane temperature was calculated by averaging the temperatures of all the particles within the vertical range −0.05a < z < 0.05a, obtaining the midplane temperature matrix T. Finally, we calculated the effective temperature in the optically thick approximation (κΣ > 1) as

$$\begin{array}{c}\hfill T_{\mathrm{eff}}^4=\frac{4}{3}\frac{T^4}{\kappa \mathrm{\Sigma }},\end{array}$$(8)

where Σ is the surface density of each element of the matrix.

Since the Toomre parameter is Q > 1, the disc tends to cool down in order to reach a self-regulated state with Q ∼ 1. Regardless of the inclusion of the radiation pressure contribution, during the first ∼400 P_B, the cooling leads to a decrease in the disc aspect ratio to H/R ∼ M_B/M_D ∼ 0.01. When radiation pressure is neglected, the disc still experiences this cooling phase but the aspect ratio increases again after 900 orbits, settling at a much higher value, around 0.07. This evolution of the disc aspect ratio is driven by the change in the midplane temperature of the disc. As can be seen from the third row of Fig. 2, the midplane temperature initially decreases and then increases again, causing the vertical expansion of the disc. On the contrary, if radiation pressure is included, the midplane temperature decreases and does not increase again. The disc aspect ratio therefore remains at around H/R ∼ 0.02, except for the very inner parts of the disc that are influenced by the effect of the stream shocks. Since the disc is initially completely radiation pressure dominated, the cooling dominates heating, causing the disc’s vertical collapse. This behaviour is consistent with what was found in (Sadowski 2016; Mishra et al. 2016; Tiwari et al. 2025, 2025b). During this cooling phase, the surface density increases in both the radiation pressure and gas-pressure-only runs regardless of the binary mass ratio and eccentricity. In the simulations that include the radiation pressure effect, the surface density reaches a stable profile on a shorter timescale. This possibly implies that the additional vertical support provided by the radiation pressure rapidly regulates both the thermal balance and geometrical structure of the disc. The effective temperature decreases from initial values between T_eff ∼ 6 × 10³ − 10⁴ K to T_eff ∼ 1.5 − 7 × 10³ K, with the hottest region lying, as expected, within R ≲ 2.5 − 3 a₀. The final effective temperature of the disc inner edge in the simulations with radiation pressure is ∼2 − 3 × 10³ K, which is a factor by at most three lower than the simulations that neglect the radiation pressure contribution. Consistently with the evolution of the disc aspect ratio, the effective temperature remains lower in the simulations that include the effect of radiation pressure.

3.2. Binary evolution

In the previous section, we showed that the radiation pressure has a non-negligible effect on the thermal properties of the disc inner edge, which ultimately regulates how much material is able to enter the cavity and accrete onto the binary. The radiation pressure contribution to the hydrostatic equilibrium of the disc also affects the migration of material from the outer to the inner regions of the circumbinary disc. This could result in a different contribution of the accretion torque, L_acc, to the total angular momentum equation, potentially changing the binary fate. Indeed, the conservation of angular momentum in our simulations implies that

$$\begin{array}{c}\hfill \frac{d\mathbf{L}}{\mathit{dt}}={\mathbf{T}}_{\mathrm{grav}}+\frac{d{\mathbf{L}}_{\mathrm{acc}}}{\mathit{dt}},\end{array}$$(9)

where L is the total binary angular momentum and T_grav is the gravitational torque (i.e. dL_grav/dt). The contribution from the accretion of gas particles onto the two binary components (L_acc) and the contribution of the gravitational torque exerted by the disc elements onto each individual massive black hole (T_grav) determine whether the binary orbit shrinks or expands as a result of the interaction with the circumbinary disc (Roedig et al. 2012).

We computed the gravitational torque exerted by N gas particles on each binary component as

$$\begin{array}{c}\hfill {\mathbf{T}}_{\mathrm{grav}}={\displaystyle \sum _{k=1}^2\sum _{i=1}^N}{\mathbf{r}}_{\mathrm{k}}\times \frac{G{M}_{\mathrm{k}}{m}_{\mathrm{i}}({\mathbf{r}}_{\mathrm{i}}-{\mathbf{r}}_{\mathrm{k}})}{|{\mathbf{r}}_{\mathrm{i}}-{\mathbf{r}}_{\mathrm{k}}{|}^3},\end{array}$$(10)

where m_i and r_i are the mass and the position of the gas particle i, M_k and r_k are the mass and the position of the sink particle k. The gas particles and the massive black hole position vectors were computed with respect to the centre of mass of the system. The accretion torque can be computed as

$$\begin{array}{c}\hfill d{\mathbf{L}}_{\mathrm{acc}}={\mathbf{r}}_i\times m_{\mathrm{i}}{\mathbf{v}}_{\mathrm{i}}-\frac{m_{\mathrm{i}}{M}_{\mathrm{k}}}{m_{\mathrm{i}}+M_{\mathrm{k}}}[({\mathbf{r}}_{\mathrm{i}}-{\mathbf{r}}_{\mathrm{k}})\times ({\mathbf{v}}_{\mathrm{i}}-{\mathbf{v}}_{\mathrm{k}})],\end{array}$$(11)

where r_i, v_i, and m_i are the position, velocity, and mass of the accreted particle i, while r_k, v_k, and M_k are the position, velocity, and mass of the sink k. We calculated the cumulative variation of the binary angular momentum over the entire duration of the simulations as

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathbf{L}={\displaystyle \sum _{\mathit{dt}}}({\mathbf{T}}_{\mathrm{grav}}dt+d{\mathbf{L}}_{\mathrm{acc}}),\end{array}$$(12)

where dt is the timestep between two consecutive snapshots. Figure 3 shows the cumulative change in angular momentum in our six simulations. We note that we can self-consistently prove that we conserve the angular momentum as we evolve a live binary and are therefore able to track its orbital parameters evolution. Indeed, since our binaries are live, we could compute the evolution of the semi-major axis, a, and eccentricity, e, directly from the simulation by using the black hole orbital energy and angular momentum at each timestep. We computed the (instantaneous) evolution of the semi-major axis and eccentricity from the binary energy and angular momentum. In particular, we measure the binary energy directly from the simulations using

Fig. 3. Angular momentum conservation for all six of our simulations: three with gas pressure only and three with the addition of radiation pressure. The first column shows the circular equal mass binary simulations, the middle column shows binaries with e = 0.45 and q = 0.7, and the right column shows equal mass highly eccentric binaries with e = 0.9. The brown and yellow lines show the gravitational and accretion torque contributions, and their sum is shown with the dotted red line for the gas-pressure-only simulations. The green line represents the binary angular momentum as measured directly from the simulation. For simulations that include radiation pressure, the same quantities are plotted with dashed cyan lines (gravitational torque), dashed purple lines (accretion torque), dotted blue lines (sum of gravitational and accretion contribution), and orange lines (binary angular momentum calculated directly from the live evolution).

$$\begin{array}{c}\hfill E(t)=\frac{1}{2}\mu (t){|\mathbf{v}|}^2-\frac{GM(t)\mu (t)}{|\mathbf{r}(t)|},\end{array}$$(13)

where μ is the binary reduced mass, and we then computed its eccentricity,

$$\begin{array}{c}\hfill e(t)=\sqrt{1+\frac{2E(t)L^2(t)}{G^2{M}^2(t){\mu }^3(t)}},\end{array}$$(14)

and semi-major axis,

$$\begin{array}{c}\hfill a(t)=\frac{L^2(t)}{GM(t){\mu }^2(t)(1-e^2(t))}.\end{array}$$(15)

Regardless of the binary properties or the implementation of radiation pressure, the total angular momentum is conserved, with only minor fluctuations in the e = 0.9, q = 1 binary case. In all the simulations, the negative gravitational torque dominates, and the total angular momentum variation is ΔL < 0, which means that the binary undergoes orbital contraction. In particular, in the circular cases, the accretion torque contribution to the total angular momentum evolution is nearly negligible. Indeed, even if ΔL_acc > 0, this positive contribution is not sufficient to counterbalance the stronger negative gravitational torque, resulting in a net angular momentum loss. We found a similar evolution in the medium eccentric case (e = 0.45). In the more eccentric case (e = 0.9), even if in the gas-pressure-only simulation the accretion contribution is higher than in the other cases, the negative gravitational torque still dominates, and the binary loses angular momentum. The inclusion of radiation pressure does not change the binary fate, and the net negative torques lead the binary to inspiral. These results are in agreement with those in Dittmann & Ryan (2021), where the negative gravitational torque contribution dominates the positive accretion torque whenever the circumbinary disc aspect ratio is small enough.

Figure 4 shows the time evolution of the semi-major axis, eccentricity, and accretion rate over time for all of our simulations. We observed that in all the cases, the binary either undergoes orbital contraction or stalls. In the circular case, the presence of radiation pressure slows down the binary inspiral. This is consistent with the fact that the accretion torque component, ΔL_acc, is slightly stronger in the simulations with radiation pressure. In the moderate eccentric case (e = 0.45), after a transient phase of ∼300 orbits, the semi-major axis is still approximately equal to one. In the highly eccentric case (e = 0.9), we found orbital contraction regardless of the presence of radiation pressure. However, in the radiation pressure case, the binary shrinks between 900 − 950 orbital periods and then stabilises, maintaining a nearly constant semi-major axis around ∼0.997. Without the radiation pressure, the semi-major axis decreases and still reaches a value of ∼0.997, but with some oscillations.

Fig. 4. Time evolution of the semi-major axis (top panel), eccentricity (centre panel), and accretion rate (bottom panel) over the last 400 orbital periods for the binary case e = 0, q = 1 (first column) and e = 0.9, q = 1 (last column) and over the last 700 orbital periods for the binary case e = 0.45, q = 0.7 (central column). The orange and the blue lines refer to the simulation with and without the radiation pressure contribution, respectively. The semi-major axis has been normalised to set a(t = 900 PB) = 1 since this is the point where our binaries start to accrete material.

In the initially circular binary, the eccentricity (e) experiences a slight increase, reaching values of ∼0.003 independent of radiation pressure presence. In the e = 0.45 case, the eccentricity increases marginally to e  ∼  0.46 − 0.47, while with e = 0.9, the eccentricity remains essentially constant. We therefore did not find a unique value at which the binary eccentricity saturates.

The accretion onto the binary is generally lower in the presence of radiation pressure, except during a transient phase of ∼400 P_B in the middle eccentric case. This result is consistent with the fact that the disc is geometrically thinner and colder, as shown in previous works (Ragusa et al. 2016; Tiede et al. 2025a) if we include radiation pressure in our simulations; see Fig. 2. The radiation pressure contribution causes the accretion rate to decrease from values ≳ Ṁ_edd to average values between Ṁ ∼0.1−0.01 M_⊙/yr. This result is consistent with what has been found in Tiwari et al. (2025) and in Tiwari et al. (2025b), where radiation pressure contribution decreases the accretion rate from 1 Ṁ/Ṁ_edd to 0.1 Ṁ/Ṁ_edd for circular equal and unequal mass binaries with a total mass of ∼10⁷ M_⊙.

3.3. Interplay between disc shape and binary preferential accretion

Here, we investigate the radiation pressure effect on the disc morphology and how this can be reflected on the accretion of gas onto the binary components. As can be seen from Fig. 1, during the binary evolution, the cavity becomes more eccentric when radiation pressure is excluded from simulations. Additionally, the precession of the cavity also appears to be different between simulations with and without radiation pressure.

To compare the evolution of the cavity eccentricity and its precession under different conditions and binary configurations, we computed the Laplace-Runge-Lenz eccentricity vector (LRL) for each gas element i:

$$\begin{array}{c}\hfill {\mathbf{e}}_{\mathrm{i}}=\frac{|{\mathbf{v}}_{\mathrm{i}}{|}^2{\mathbf{r}}_{\mathrm{i}}-({\mathbf{v}}_{\mathrm{i}}·{\mathbf{r}}_{\mathrm{i}}){\mathbf{v}}_{\mathrm{i}}}{G{M}_{\mathrm{B}}}-\widehat{{\mathbf{r}}_{\mathrm{i}}}.\end{array}$$(16)

Following this definition, we computed the orbital eccentricity as

$$\begin{array}{c}\hfill |\mathbf{e}|=\sqrt{e_{\mathrm{x}}^2+e_{\mathrm{y}}^2},\end{array}$$(17)

where e_x and e_y are the Cartesian components of the vector. We determined the longitude of the pericentre of the disc, ϖ_d, as

$$\begin{array}{c}\hfill {\varpi }_{\mathrm{d}}=arctan\frac{e_{\mathrm{y}}}{e_{\mathrm{x}}}.\end{array}$$(18)

We computed the mass weighted eccentricity vector between r = a and r = 4a of all the gas elements with a bound orbit (i.e. E_i < 0) over 1300 and 1600 orbital periods for the equal and unequal mass simulations, respectively.

Figure 5 shows the evolution of the disc eccentricity magnitude (top panels) and the longitude of pericentre (bottom panels) for the six simulations. The orange and blue lines refer to the simulations with and without the implementation of radiation pressure, respectively. In the circular binary case, when the radiation pressure contributes to the hydrostatic equilibrium of the disc, the cavity remains circular and the precession is not significant, unlike in the case without the radiation pressure. In the eccentric binary cases, the radiation pressure effect is less dramatic, and the disc cavity becomes eccentric and undergoes precession. The evolution of the cavity eccentricity is similar in the simulations with and without the radiation pressure. This suggests that radiation pressure has an important effect on counteracting cavity eccentricity growth and precession in circular equal mass systems. While the effect of radiation pressure is to significantly reduce the disc temperature and locally isothermal simulations find more eccentric cavities for colder discs, we did not find any disc eccentricity growth in our circular binary simulation in the presence of radiation pressure. The main difference is indeed the inclusion of radiation pressure in our work, but we caution that to pinpoint the exact origin of this eccentricity suppression, we would need to carry out more of these expensive simulations in order to investigate the dependence of this result from the choice of initial conditions of the circumbinary disc. We note that a similar result was also obtained with radiation MHD simulations at smaller binary separations (Tiwari et al. 2025, their Fig. 9).

Fig. 5. Eccentricity vector evolution over time for the simulations with e = 0, q = 1 binaries (left column); e = 0.45, q = 0.7 binaries (centre column); and e = 0.9, q = 1 binaries (right column). The top panels show the magnitude, and the bottom panels show the phase in radians. The orange and blue lines refer to the simulation with and without the radiation pressure contribution, respectively.

We investigated how the precession of the binary with respect to the precession of the cavity affects the accretion dynamics. In fact, the angle between the binary eccentricity vector and the disc cavity eccentricity vector plays a crucial role in determining the amount of preferential accretion experienced by each binary component (Dunhill et al. 2015; D’Orazio & Duffell 2021; Siwek et al. 2022). When this angle settles to a nearly constant value over time, due to synchronised precession between the binary and the disc, we refer to it as a locking angle. In particular, Siwek et al. (2022) showed that non-zero locking angles that depend on both the binary mass ratio and eccentricity can enhance or dampen preferential accretion. For instance, apsidal misalignment between the binary and the cavity eccentricity vectors, especially at higher eccentricities (e_B > 0.6), may reduce accretion efficiency, slowing mass ratio growth. Conversely, apsidal alignment at lower eccentricities leads to more efficient accretion, promoting mass ratio equalisation.

We computed the alternation of the accretion rate onto each black hole and the preferential accretion ratio λ = Ṁ₂/Ṁ₁ across all the binary configurations with and without the implementation of radiation pressure. The results are shown in the first two rows of Fig. 6. For the circular equal mass system, the accretion onto the binary components tends to be symmetric. We found λ = 1 both with and without the implementation of radiation pressure, which is consistent with the expected value from the empirical relation λ_fit = q^−0.9 (Duffell et al. 2020) and in agreement with previous gas-only works (e.g. Farris et al. 2014; Muñoz et al. 2020; Duffell et al. 2020; Dittmann & Ryan 2021). For the highly eccentric equal mass binaries e = 0.9, q = 1, we obtained λ = 1.2 and λ = 0.98 with only gas pressure and with radiation pressure, respectively. This result is broadly consistent with Siwek et al. (2022).

Fig. 6. Time evolution of the individual accretion rate (first and second row) and the relative precession phase between the disc and the binary (ϖ = ϖd − ϖb; third and last row) for the equal mass circular case e = 0, q = 1 (left column); e = 0.45, q = 0.7 (middle column); and the equal mass highly eccentric one e = 0.9, q = 1 (right column). The first two rows show the accretion onto the primary (green line) and the secondary (pink line) black hole over the last 400 and 200 orbits for the equal and unequal mass cases, respectively. The upper row shows the results of simulations without the radiation pressure, while the second row includes the radiation pressure. The third and last row display the evolution of the relative precession phase over the radius. The third row shows results without radiation pressure, and the last row includes it.

The most interesting regime to investigate was the unequal mass binary case, as prolonged preferential accretion is expected to occur in unequal mass systems. In the intermediate eccentricity unequal mass case, we found λ = 2.28 and λ = 1.18 without and with the radiation pressure contribution in the simulations. The stronger preferential accretion we observed in the gas-only pressure case results from the favourable locking angles between the binary and cavity eccentricity vectors, which significantly enhances the accretion onto the secondary binary component. We note that Siwek et al. (2022), for a similar case with e = 0.4 and q = 0.7, reported λ ≈ 2. The inclusion of the radiation pressure again leads to a more balanced accretion rate between the binary components. Our results are somewhat in contrast with the RMHD simulations of a circular unequal-mass binary presented in Tiwari et al. (2025b), where they found that radiation pressure reduces the total accretion rate but does not suppress preferential accretion onto the secondary (λ ∼ 1.8).

To better understand the effect of the behaviour of the alternation of the accretion rate onto each black hole and the preferential accretion, we analysed the relative precession phase between the disc and the binary, defined as ϖ = ϖ_d − ϖ_b. We computed the binary eccentricity vector e_b and its phase ϖ_b following Eqs. (16)–(18). We binned and averaged the LRL vectors over the mass in the radial direction before computing the disc eccentricity vector phase using Eq. (18). The third and last row of Fig. 6 show the time evolution of this relative precession for all of our binary configurations. Note that the bottom-left panels of the third and last row have little meaning, as the disc is quasi-circular. In the circular equal mass case (left panel), the phase difference between the binary eccentricity vector and the disc eccentricity vector exhibits a periodic behaviour, repeating over ∼200 orbital periods. When radiation pressure is included, the disc eccentricity growth is strongly suppressed (see Fig. 5), and the disc experiences a periodic behaviour on a shorter timescale. The disc eccentricity remains small, and therefore, even with an alternating pattern in the relative angle, this does not affect accretion. In the highly eccentric equal mass scenario (i.e. e = 0.9, q = 1) shown in the rightmost panel of Fig. 6, the angular difference exhibits a similar behaviour in the gas-only and radiation pressure simulations. In both cases, ϖ precesses over approximately 260 orbital periods. In the moderate eccentric unequal mass case (i.e. e = 0.45 q = 0.7) shown in the middle panel of Fig. 6, the relative precession phase evolves more slowly, with a full cycle spanning ∼400 orbital periods. The full cycle duration was determined by analysing the relative precession phase over a longer time window, while Fig. 6 shows results in a narrow time window (t = 1400 − 1600 P_B) that is consistent with the time window used in other figures. This results in a temporary locking phase associated with a strong preferential accretion onto the secondary component. In this simulation, we found λ ∼ 2. In this case, the inclusion of radiation pressure does not lead to a significant change in the difference between the relative precession phase between the disc and the binary.

To further investigate the dynamics of the streams that provide material to the binary component during the alternation of the accretion rate onto each black hole and the preferential accretion phases, which are essentially determined by the angle between the disc and binary eccentricity vectors, we show the disc surface density maps for the two eccentric binary simulations in Fig. 7. From top to bottom, the panels correspond to angles of ϖ ∼ −π, −π/2, 0, except for the e = 0.45 radiation pressure case, where we report 0, π/4, π/2. These configurations respectively represent accretion preferentially onto the primary black hole, symmetric accretion onto both components, and accretion predominantly onto the secondary.

Fig. 7. Surface density maps in the x − y plane of the circumbinary disc around the equal mass eccentric binary e = 0.9 and q = 1 (left panel), and the unequal mass eccentric binary e = 0.45 and q = 0.7 (right panel) at different relative precession angles. The left and right columns respectively present the simulations without and with the radiation pressure. The solid blue and orange circles indicate the primary and the secondary black hole, respectively, while the semi-transparent circles marks each orbit. The green arrows indicate the disc eccentricity vector at different radii.

In the eccentric equal mass case (e = 0.9, q = 1), regardless of the presence of radiation pressure, when the disc and the binary eccentricity vectors are anti-parallel, in agreement with Fig. 6, the accretion is preferentially onto the black hole located near the cavity edge. When the vectors are perpendicular, the two binary components accrete the same amount of gas. Finally, when the disc and the binary eccentricity vectors are parallel, accretion is favoured onto the black hole further away from the cavity edge.

In the middle eccentric case (e = 0.45, q = 0.7), the behaviour is more complex. In the gas-pressure-only simulation, we find a similar result to the equal mass case: when ϖ_d and ϖ_b are anti-parallel, accretion is preferentially onto the primary black hole, and when they are aligned, it favours the secondary. However, the configuration at π/2 results in a preferential inflow onto the secondary, which is located closer to the cavity edge. When the radiation pressure is included, the trend is reversed: Accretion occurs preferentially onto the primary when the vectors ϖ_d and ϖ_b are aligned. As in the only-gas simulation, for ϖ ∼ π/2, the black hole closer to the cavity accretes more gas from the circumbinary disc.

4. Conclusions

In this work, we have studied the impact of radiation pressure in the evolution of accreting MBHBs at milli-parsec-scales embedded in thin circumbinary gaseous discs using 3D numerical simulations with hyper-Lagrangian refinement. We described the discs using an adiabatic equation of state and allowed the gas to change its temperature through shocks, PdV work, and black-body radiation. We explored binary eccentricities e = 0, 0.45, 0.9 and mass ratios q = 1, 0.7. We investigated the role that radiation pressure plays in the evolution of the disc’s aspect ratio, H/R; surface density, Σ; and effective temperature, T_eff. Additionally, we measured the torques exerted by the disc onto the binary and the effect they have on the orbital evolution of the binary by directly measuring the change in eccentricity, e; semi-major axis, a; and accretion rate. Additionally, we studied the alternation of the accretion rate onto each black hole and the preferential accretion, focusing on its dependence on the angle between the disc and the binary eccentricity vectors.

We find that radiation pressure suppresses the growth of cavity eccentricity and its precession in circular equal-mass binaries, while it has a negligible impact in eccentric systems. In the circular equal mass case, the effect of radiation pressure is to effectively suppress disc eccentricity growth and therefore prohibit the formation of the ‘lump’, i.e. over-density at the cavity inner edge. This result is in agreement with a recent RMHD simulation of circular equal-mass binaries presented in Tiwari et al. (2025), where they found that radiation contribution leads to a less eccentric cavity and a weaker over-density at the inner edge of the disc due to the reduced impact of streams and overall disc cooling.

Conversely, we find that if the binary is eccentric, the effect of radiation pressure is to lead to a more pronounced lump (see Fig. 1). This is consistent with the disc aspect ratio being significantly lower in the simulations with radiation pressure, which, consistent with other studies (Franchini et al. 2024a), contributes to enhancing the over-density contrast. In these cases, the high binary eccentricity can overcome the effect of radiation pressure suppressing the cavity eccentricity growth, while in the circular binary case, radiation pressure maintains the circularity of the cavity.

Radiation pressure plays a crucial role in regulating both the vertical structure and thermal evolution of the circumbinary disc. In all simulations, the disc initially cools over the first ∼400 orbits, reducing its aspect ratio from H/R = 0.1 to H/R ∼ 0.01, as it attempts to reach a self-regulated state with a Toomre parameter of Q ∼ 1. When radiation pressure is included, the disc aspect ratio remains at a value of ∼0.01, except for the innermost regions, where H/R ∼ 0.05 due to shocks driven by gas streams that are flung back from the binary and impact the cavity edge. Since the disc is initially completely radiation pressure dominated, the cooling dominates the heating, and the disc collapses in the vertical direction, consistent with previous more sophisticated numerical simulations in the literature (Sadowski 2016; Mishra et al. 2016). The presence of radiation pressure leads to effective temperatures up to a factor of three lower than in simulations with only gas pressure. The accretion rate onto the binary is lower when radiation pressure is included (Ṁ ∼0.01−0.1 M_⊙/yr), which is consistent with the colder thinner disc structure (Ragusa et al. 2016; Tiede et al. 2025a; Tiwari et al. 2025, 2025b).

Since we evolved the binary orbit with time, we could track the evolution of the binary orbital elements and assess the dynamical impact of radiation pressure. We find that radiation pressure does not significantly affect the conservation of total angular momentum. The gravitational torque exerted by the disc is always negative, and this is consistent with the fact that our circumbinary discs self-regulate towards an aspect ratio lower than 0.1, i.e. the value considered to be the threshold for binary outspiral (Duffell et al. 2024). The torque due to the accretion of particles onto the binary components is always positive and slightly higher in the only-gas-pressure simulation very eccentric e = 0.9 case. However, its magnitude is never sufficiently large to counteract the negative gravitational torque, and we therefore find our binaries to always decrease their semi-major axis.

The binary eccentricity remains mostly constant, with only minor variations (e ∼ 0.003 in the circular case and e ∼ 0.46 in the e = 0.45 and q = 0.7 case). We therefore did not find convergence towards a unique value regardless of the presence of radiation pressure in the simulations. While it is possible that the systems have not yet reached an equilibrium eccentricity due to the limited duration of the runs, we note that our results are consistent with those of Siwek et al. (2023), who similarly found no evolution of the binary eccentricity in the circular equal mass case and an eccentricity growth in the e = 0.4, 0.5 and q = 0.7 cases over a longer timescale.

We investigated the interplay between the disc morphology and the alternation of accretion onto each black hole. In the unequal mass binary case we explored, the radiation pressure contribution reduces the preferential accretion onto the secondary, resulting in a more balanced mass growth between the binary components. In particular, in the gas-pressure-only eccentric e = 0.45 simulation, we find the ratio between the secondary and the primary black hole to be λ = Ṁ₂/Ṁ₁ ∼ 2.3, in agreement with Siwek et al. (2022). However, when radiation pressure is included, the accretion is more balanced and λ = 1.18, and the mass ratio growth is effectively suppressed. This is in contrast with the RMHD simulations of a circular unequal-mass binary presented in Tiwari et al. (2025b), where they found that radiation pressure does not suppress preferential accretion onto the secondary. In order to better understand the behaviour of the alternation of the accretion rate onto each black hole and the preferential accretion, we analysed the relative precession phase between the disc and the binary defined as ϖ = ϖ_d − ϖ_b. In equal mass binary cases, we find that ϖ cycles over ∼200 and ∼260 orbital periods for eccentricities of e = 0 and e = 0.9, respectively, and with shorter cycles in cases with the radiation pressure. For the eccentric unequal mass binary (e = 0.45 and q = 0.7), ϖ evolves more slowly, allowing temporary locking phases between the disc and the binary. These phases correspond to episodes of enhanced accretion onto one binary component over time. In the radiation pressure simulations, accretion occurs onto the secondary black hole when ϖ ∼ 0 (i.e. disc and binary eccentricity vectors are aligned). For intermediate angles such as π/2, we find that preferential accretion can still occur in the unequal mass case, where the black hole closer to the cavity tends to accrete more.

Finally, we note that since radiation pressure in hydrodynamics simulations of MBHBs regulates the thermal structure of the circumbinary disc, it may also impact the electromagnetic emission of MBHBs, possibly affecting the periodic signatures in emitted light used to identify these objects in observations. This is the subject of a follow up work (Cocchiararo et al. 2025).

Acknowledgments

We thank Julian Krolik, Tamara Bogdanović and Maria Charisi for their feedback that improved the clarity of this manuscript. We thank Daniel Price for providing the PHANTOM code for numerical simulations and acknowledge the use of SPLASH (Price 2007) for the rendering of the figures. We thank Phil Hopkins for providing the GIZMO code for numerical simulations. AS acknowledges financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691). AS also acknowledges the financial support provided under the European Union’s H2020 ERC Advanced Grant “PINGU” (Grant Agreement: 101142079). AL acknowledges support by the PRIN MUR “2022935STW”. AF acknowledges financial support from the Unione europea – Next Generation EU, Missione 4 Componente 1 CUP G43C24002290001.

References

All Figures

Fig. 1. Surface density map in code units (1 code unit = 200 AU) in the x − y plane for the simulations with e = 0, q = 1 (upper panel); e = 0.45, q = 0.7 (middle panel); and e = 0.9, q = 1 (bottom panel) at three different moments: t = 1,  900,  1300 PB for the equal mass cases and t = 1,  900,  1600 PB for the unequal mass case. In each panel, the first and the second row show results for Ptot = Pgas and Ptot = Pgas + Prad, respectively. In the eccentric simulations, the initial cavity radius is 3a0, while in the circular binary case, it is 2a0. When the radiation pressure is included, the binary begins accreting at the very early stages of its evolution; otherwise, the disc experiences a transient phase lasting approximately 900 PB, 700 PB, and 600 PB respectively for the cases e = 0, q = 1; e = 0.45, q = 0.7; and e = 0.9, q = 1. The spatial resolution at the edge of the cavity (r = 3a) is Δx = 0.025a.

In the text

Fig. 2. Time evolution of the surface density (first row), disc aspect ratio (H/R; second row), the midplane temperature (third row), and the effective temperature (last row) as a function of radius for the binary case e = 0, q = 1 (first column); e = 0.45, q = 0.7 (central column); and e = 0.9, q = 1 (last column) at different times (distinguished by different colours). The solid and dashed lines refer to the simulation without and with the implementation of the radiation pressure, respectively. The initial condition of the disc is the same across all the cases except for the initial cavity radius, which is 3a0 in both of the eccentric simulations, while it is 2a0 in the circular binary cases. We calculated the midplane temperature and the effective temperature assuming that both the gas and the radiation pressure contribute to the hydrostatic equilibrium of the disc in all the cases. The radiation pressure contribution leads the disc to maintain a lower aspect ratio and temperature.

In the text

Fig. 3. Angular momentum conservation for all six of our simulations: three with gas pressure only and three with the addition of radiation pressure. The first column shows the circular equal mass binary simulations, the middle column shows binaries with e = 0.45 and q = 0.7, and the right column shows equal mass highly eccentric binaries with e = 0.9. The brown and yellow lines show the gravitational and accretion torque contributions, and their sum is shown with the dotted red line for the gas-pressure-only simulations. The green line represents the binary angular momentum as measured directly from the simulation. For simulations that include radiation pressure, the same quantities are plotted with dashed cyan lines (gravitational torque), dashed purple lines (accretion torque), dotted blue lines (sum of gravitational and accretion contribution), and orange lines (binary angular momentum calculated directly from the live evolution).

In the text

Fig. 4. Time evolution of the semi-major axis (top panel), eccentricity (centre panel), and accretion rate (bottom panel) over the last 400 orbital periods for the binary case e = 0, q = 1 (first column) and e = 0.9, q = 1 (last column) and over the last 700 orbital periods for the binary case e = 0.45, q = 0.7 (central column). The orange and the blue lines refer to the simulation with and without the radiation pressure contribution, respectively. The semi-major axis has been normalised to set a(t = 900 PB) = 1 since this is the point where our binaries start to accrete material.

In the text

	Fig. 5. Eccentricity vector evolution over time for the simulations with e = 0, q = 1 binaries (left column); e = 0.45, q = 0.7 binaries (centre column); and e = 0.9, q = 1 binaries (right column). The top panels show the magnitude, and the bottom panels show the phase in radians. The orange and blue lines refer to the simulation with and without the radiation pressure contribution, respectively.
In the text

Fig. 6. Time evolution of the individual accretion rate (first and second row) and the relative precession phase between the disc and the binary (ϖ = ϖd − ϖb; third and last row) for the equal mass circular case e = 0, q = 1 (left column); e = 0.45, q = 0.7 (middle column); and the equal mass highly eccentric one e = 0.9, q = 1 (right column). The first two rows show the accretion onto the primary (green line) and the secondary (pink line) black hole over the last 400 and 200 orbits for the equal and unequal mass cases, respectively. The upper row shows the results of simulations without the radiation pressure, while the second row includes the radiation pressure. The third and last row display the evolution of the relative precession phase over the radius. The third row shows results without radiation pressure, and the last row includes it.

In the text

Fig. 7. Surface density maps in the x − y plane of the circumbinary disc around the equal mass eccentric binary e = 0.9 and q = 1 (left panel), and the unequal mass eccentric binary e = 0.45 and q = 0.7 (right panel) at different relative precession angles. The left and right columns respectively present the simulations without and with the radiation pressure. The solid blue and orange circles indicate the primary and the secondary black hole, respectively, while the semi-transparent circles marks each orbit. The green arrows indicate the disc eccentricity vector at different radii.

In the text