EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 707 (March 2026) A&A, 707 (2026) A369 Full HTML
Open Access
Issue
A&A
Volume 707, March 2026
Article Number A369
Number of page(s) 8
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202558109
Published online 19 March 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

In the next decade, the Laser Interferometer Space Antenna (LISA) mission is expected to detect gravitational wave (GW) signals in the millihertz regime emitted by massive black hole binaries (MBHBs) approaching a merger (Kocsis et al. 2006; Dotti et al. 2006; Mayer et al. 2007; Amaro-Seoane et al. 2012; Capelo et al. 2015; Amaro-Seoane et al. 2017). These sources reside in the centre of galaxies and are expected to evolve and eventually merge in gas-rich environments. The presence of gas and the highly dynamic space-time suggest the possibility of accretion onto the binary components, resulting in the production of strong electromagnetic (EM) signals. Therefore, MBHBs close to merging represent the perfect scenario for multi-messenger astronomy with massive compact objects. The coincident detection of a GW signal and an EM counterpart could provide unique information about the system’s environment and evolution.

Furthermore, the pulsar timing array (PTA) community recently announced the detection of a stochastic GW background signal in the nanohertz frequency band (EPTA Collaboration 2023; Agazie et al. 2023; Reardon et al. 2023; Xu et al. 2023), likely produced by a population of 108 − 109 M⊙ black hole (BH) binaries. These PTA results showcase our deep understanding of the GW signal produced by massive compact objects while also providing a first proof of the existence of close-separation MBHBs.

Still, there is considerable uncertainty about the signature features of an EM signal produced by a MBHB that make the source distinguishable from a single active galactic nucleus (Bogdanović et al. 2011; Dotti et al. 2012; Bogdanović et al. 2022) and that cannot be explained by other mechanisms, for example asymmetric emissions or a periodic modulation of an active galactic nucleus disc (Rigamonti et al. 2025; Dotti et al. 2023). Major differences may exist in the accretion flows close to the central binary. The non-axisymmetric potential of a binary could carve an empty cavity around the BHs due to strong torques while thin streams of gas fall into the centre and are accreted by the compact objects (MacFadyen & Milosavljević 2008; Noble et al. 2012), yielding characteristic EM features. The dynamics of these processes may depend strongly on the spins of the BHs and on the properties of the gaseous environment (Bowen et al. 2018; Bode et al. 2010; Cattorini et al. 2021; Ruiz et al. 2023). Variability in the emitted light curves is also observed in many simulations (Krauth et al. 2023; Gutiérrez et al. 2022; Franchini et al. 2024a; Cocchiararo et al. 2024, 2025) and could be a smoking gun for detecting an evolving MBHB with optical facilities such as the Vera Rubin Observatory’s Large Synoptic Survey Telescope (LSST) (Chiesa et al. 2025).

Numerical simulations are powerful tools for studying and understanding such extreme environments and can help predict the EM emission from a MBHB approaching a merger. Many studies aim to reach a stationary regime to study the properties of a circumbinary disc. Therefore, starting from idealised initial conditions, it is necessary to simulate the system for hundreds of orbits, as discussed in the binary disc code comparison presented in Duffell et al. (2024). Hence, many studies have focused on long-term Newtonian simulations, both in 2D and in 3D (Muñoz et al. 2020; Krauth et al. 2023), also considering post-Newtonian (PN) corrections to the binary evolution (Franchini et al. 2024a).

On the other hand, many interesting accretion features exclusive to a MBHB system may reside close to the event horizons (d’Ascoli et al. 2018; Gold 2019; Cattorini & Giacomazzo 2024). To accurately explore those regions, 3D numerical relativity simulations are necessary. However, these kinds of simulations are extremely computationally expensive and are therefore limited to evolving BH binaries only for a few orbits before the merger (Bode et al. 2010; Giacomazzo et al. 2012; Gold et al. 2014; Paschalidis et al. 2021; Fedrigo et al. 2024). Only recently have hybrid and approximated methods been employed to solve this issue and evolve the system for a longer time (Mignon-Risse et al. 2023; Avara et al. 2024; Ennoggi et al. 2025a,b).

A simpler solution is presented by Combi et al. (2021) and Combi & Ressler (2026) in the form of an approximated dynamic metric created by the superposition of two boosted BHs in Kerr-Schild (KS) coordinates, called the superposed Kerr-Schild (SKS) metric. It represents an accurate yet computationally cheap way to evaluate the dynamic metric of a BH binary, and it can be used to perform general relativistic magnetohydrodynamic (GRMHD) or ray-tracing simulations. Thanks to the high versatility of the SKS metric, it is now possible to explore in depth the evolution and dynamics of the gaseous environment in which MBHBs reside using different state-of-the-art astrophysical codes.

In this work we present the implementation of the SKS metric in the relativistic meshless version of the GIZMO code (Hopkins 2015; Lupi 2023; Fedrigo & Lupi 2025). Our goal is to build a code that can evolve relativistic circumbinary discs around BH binaries for thousands of orbits, to study their properties and search for unique signatures.

In Sect. 2.1 we briefly review the SKS formulation, and in Sect. 2.2 we present some details regarding its implementation in GIZMO. Then, in Sect. 3, we report results from two tests we performed to verify the metric was correctly implemented and verify its stability within the code. The first test compared our results with a numerical relativity simulation of a uniform gas distribution accreting onto a BH binary (Sect. 3.1). With the second test, we explored the differences between a Newtonian and a relativistic circumbinary disc (Sect. 3.2). Finally, we draw our conclusions in Sect. 4.

2. Black hole binary metric

In this section we describe the implementation of the SKS metric in GIZMO. The relativistic scheme implemented in GIZMO is presented in Lupi (2023, hydrodynamics) and Fedrigo & Lupi (2025, magnetohydrodynamics). Note that we do not consider magnetic fields in the tests discussed in this work, deferring the analysis of magnetised plasma in a MBHB environment to a future study. The superposed binary metric was first introduced in Lopez Armengol et al. (2021), later refined in Combi et al. (2021), and finally extended to include the merger phase in Combi & Ressler (2026, hereafter Co24). Here, we briefly review the main concepts of the Co24 metric formulation, and we then focus on some implementation details. Throughout this work, we employe natural units, i.e. G = c = 1.

2.1. SKS formulation

The metric of a single spinning BH of mass M and spin a can be expressed in the ‘Cartesian’ KS coordinate system {Xμ}:

g μ ν = η μ ν + 2 H l μ l ν , Mathematical equation: $$ \begin{aligned} g_{\mu \nu } = \eta _{\mu \nu } + 2\mathcal{H} l_\mu l_\nu , \end{aligned} $$(1)

where the null vector l is defined by

l μ d X μ = d t + 1 r 2 + a 2 [ r X i − ϵ jk i a j X k + a i X j δ ij r a i ] d X i , Mathematical equation: $$ \begin{aligned} l_\mu \mathrm{d}X^\mu = \mathrm{d}t + \frac{1}{r^2+a^2}\left[rX^i - \epsilon ^i_{jk}a^jX^k + \frac{a^iX^j\delta _{ij}}{r}a^i\right]\mathrm{d}X_i, \end{aligned} $$(2)

the function ℋ is

H ( X i ) = M r 3 r 4 + ( a i X j δ ij ) 2 , Mathematical equation: $$ \begin{aligned} \mathcal{H} (X^i) = \frac{Mr^3}{r^4+(a^iX^j\delta _{ij})^2}, \end{aligned} $$(3)

the Boyer-Lindquist radius is

r = R 2 − a 2 + ( R 2 − a 2 ) 2 + 4 ( a i X j δ ij ) 2 2 , Mathematical equation: $$ \begin{aligned} r=\sqrt{\frac{R^2-a^2+\sqrt{(R^2-a^2)^2+4(a^iX^j\delta _{ij})^2}}{2}}, \end{aligned} $$(4)

and R2 = X2 + Y2 + Z2.

Following Co24, a generalised Lorentz boost transformation is then applied to the KS metric (Eq. 1), going from the non-inertial frame at rest with the BH {Xμ} to an inertial frame where the BH is moving on an inspiralling orbit, {xμ}. Defining the BH trajectory with s(t), the coordinate transformation with velocity v = βn(t) is given by

T = t − W β ( n x x ¯ + n y y ¯ + n z z ¯ ) , X = x ¯ [ 1 + n x 2 ( W − 1 ) ] + y ¯ ( W − 1 ) n x n y + z ¯ ( W − 1 ) n x n z , Y = y ¯ [ 1 + n y 2 ( W − 1 ) ] + x ¯ ( W − 1 ) n x n y + z ¯ ( W − 1 ) n y n z , Z = z ¯ [ 1 + n z 2 ( W − 1 ) ] + y ¯ ( W − 1 ) n y n z + x ¯ ( W − 1 ) n x n z , Mathematical equation: $$ \begin{aligned}&T = t-W\beta (n_x\bar{x} + n_y\bar{y} + n_z\bar{z}),\nonumber \\&X = \bar{x}[1+n_x^2(W-1)] + \bar{y}(W-1)n_xn_y + \bar{z}(W-1)n_xn_z,\nonumber \\&Y = \bar{y}[1+n_y^2(W-1)] + \bar{x}(W-1)n_xn_y + \bar{z}(W-1)n_yn_z,\nonumber \\&Z = \bar{z}[1+n_z^2(W-1)] + \bar{y}(W-1)n_yn_z + \bar{x}(W-1)n_xn_z, \end{aligned} $$(5)

where x ¯ i = x i − s i ( t ) Mathematical equation: $ \bar{x}^i = x^i - s^i(t) $ and W is the Lorentz factor. To conclude, the superposed metric of two inspiralling spinning BHs in KS coordinates can be written as

g μ ν = η μ ν + [ 2 ∂ X ρ ∂ x μ ∂ X σ ∂ x ν H l ρ l σ ] 1 + [ 2 ∂ X ρ ∂ x μ ∂ X σ ∂ x ν H l ρ l σ ] 2 . Mathematical equation: $$ \begin{aligned} g_{\mu \nu } = \eta _{\mu \nu } + \left[2\frac{\partial X^\rho }{\partial x^\mu }\frac{\partial X^\sigma }{\partial x^\nu }\mathcal{H} l_\rho l_\sigma \right]_1 + \left[2\frac{\partial X^\rho }{\partial x^\mu }\frac{\partial X^\sigma }{\partial x^\nu }\mathcal{H} l_\rho l_\sigma \right]_2. \end{aligned} $$(6)

This approximate metric formulation exhibits good accuracy and low constraint violation for intermediate binary separations, but it begins to fail at BH separations ≲8M. The SKS metric performs well and is smooth at the innermost stable circular orbit regions, which are important features to take into consideration in order to correctly evolve the accretion disc dynamics (see the space-time analysis in Combi et al. 2021). While the Co24 metric formulation also encompasses the merger phase, we did not take this feature into consideration in this work, with the exception of the Bondi-like test in which we evolved the binary until the merger moment (Sect. 3.1).

2.2. SKS implementation in GIZMO

The SKS formulation, as presented in Co24, has been coded in a public, ready-to-use C implementation (Combi & Ressler 2024) that takes the BH parameters as input data and returns the covariant metric from Eq. (6) at a given position. At each timestep, the needed inputs are the positions, velocities, spins, and masses of each BH. In our implementation, we took advantage of the BH particle structure already present in the Newtonian version of the code (see Franchini et al. 2022, 2023, 2024a), from which we exported the BH parameters. Hence, the binary is described by two BH particles that are evolved live, and their properties can be influenced by the interaction with the surrounding gas. Simulations presented in Franchini et al. (2023) show the importance of evolving live binaries, highlighting crucial differences in the conservation of the total angular momentum of the system, in the binary separation evolution, and in the effects of gravitational torques. In this work, however, we neglected the change in BH masses and spin driven by the gaseous disc for simplicity. This choice is reasonable as we assumed our discs to have much lower masses compared with the binary and therefore to have a negligible effect on its dynamical evolution. The focus here is instead to investigate the changes that the general relativity treatment brings to the gas dynamics.

We integrate the equation of motion considering PN corrections up to the 2.5 PN order, i.e. including the contribution of the GW emission to the orbital evolution (Blanchet 2014), to evaluate the BH trajectories. We refer to Franchini et al. (2024a) for a complete presentation of the PN implementation in GIZMO1.

During the system evolution, gas can be accreted by the BHs. In our implementation, at each timestep, the code evaluates the distance of each gas particle to both BHs and identifies the closest one. When such a distance becomes smaller than a user-defined excision radius, the gas particle is removed from the simulation. Since the superposed metric is written in the horizon-penetrating KS coordinate system, we usually set the excision radius to reside inside the event horizon. As mentioned above, in this work we did not change the BH properties when gas particles are accreted onto the horizons. This represents a sufficiently good approximation if MBH ≫ Mgas, as we assumed in our work. To keep track of the accretion onto the BH during the simulation, at each timestep we also store the properties of the particles crossing the closest BH’s event horizon radius.

We used the on-the-fly adaptive particle splitting approach, as presented in Franchini et al. (2022). With this procedure, particles are split according to a scheme that depends on the particle distance to the domain centre. In practice, the mass resolution is doubled each time the gas crosses a set of defined radii, to better resolve the disc edge and the accretion streams inside the inner disc cavity. We refer to Franchini et al. (2022) for an in-depth presentation of this refinement approach. The only difference with the relativistic version is that, in the latter, while the two child particles inherit the hydrodynamic properties of the parent one, the metric coefficients are recomputed at the location of the two new particles.

Finally, when the resolution is poor, gas particles close to the BH’s event horizon might cause the neighbour search to be extended to particles located on the other side of the BH. Such a case would either produce spurious results (as the gas cells should be causally disconnected from each other) or, in the worst-case scenario, completely break the code. To avoid such events in our algorithm, the neighbour search in our implementation is limited to a maximum distance equal to the current distance of the particle from the closest BH. Note that we performed accretion tests while varying the excision radius and with different neighbour search prescriptions, for example letting the kernel extend over the BH but preventing particle pairs whose line of sight is directly obscured by a BH from interacting. All the tests yielded comparable results, without any statistically significant differences.

3. Metric validation tests

In this section we present the two general relativistic (GR) hydrodynamic tests we performed to confirm the correct implementation of the SKS metric in our code2. We studied a set of inspiralling binaries at close and intermediate separations, surrounded by two different gaseous environments: a uniform plasma distribution and a thin circumbinary disc. We compared our results with those obtained using different evolution schemes. We focused on a comparison of the gas distribution and the mass accretion rate onto the horizons.

3.1. Bondi accretion

As a first test, we evolved an equal-mass binary initially placed on a circular orbit with a separation of 12M. For simplicity, in this test we did not include any live BH particle, but we defined the BH location through pre-computed trajectories using CBWAVES, which included PN corrections up to order 3.5. We embedded the binary in a homogeneous gas distribution initially at rest extending up to 1500M and let the system evolve towards a Bondi-like solution around the binary system, similarly to the un-magnetised case studied by Cattorini et al. (2021). Having such a large domain is important as we want the simulation to last for several orbits without including any particle injection at the outer boundary. At the same time, sampling the gas distribution with enough resolution to properly resolve the event horizon with several gas cells is crucial, as preliminary tests at low resolution showed spurious density fluctuations at the centre of the domain when the number of cells was too low.

Fulfilling both requirements at the same time would require an extremely high number of particles, which would make the simulation unfeasible. To overcome this limitation, we decided to include adaptive particle splitting during the simulation, which operates between 50M and 1000M, increasing the resolution by a factor of 210 following a parabolic splitting function (see Franchini et al. 2022 for details). Moreover, we also mapped the initial conditions using the same splitting function employed during the simulation, which resulted in 20 million particles with variable masses, instead of the tens of billions of particles required with a constant mass resolution everywhere. As we neglected the contribution of the gas self-gravity (see Franchini et al. 2021, 2024b for discs with self-gravity), the particle splitting does not introduce any spurious effect in the calculations while guaranteeing a very high resolution close to the binary from the beginning. To prevent particles at the outer edge from leaving the domain, we further constrained the location and the hydrodynamic properties of the outermost five particle spherical shells to remain constant over time.

For this specific test, we evolved the system for the last 5 hours before the merger, which occurs at t = 1772 for a binary with M = 106 M⊙, assuming an ideal gas law with Γ = 4/3. In Fig. 1 we show the accretion rate onto each BH measured during the simulation, compared with the numerical relativity simulation results from Cattorini et al. (2021). Given the symmetry of the system, the accretion rates onto the BHs are exactly the same, as expected. After an initial rise, the accretion rate reaches a plateau that corresponds to the Bondi-like solution, in agreement with the results from Cattorini et al. (2021). At later times, both simulations exhibit an increasing accretion rate as the two BHs approach the merger. Compared with the numerical relativity results, the approximated metric slightly overestimates the gravitational pull onto the BHs, which results in a roughly 10 percent higher accretion rate (Combi & Ressler 2026). Considering (i) the different space-time evolution in the two cases, (ii) differences in the numerical techniques employed and the discretisation scheme, and (iii) the maximum resolution close to the horizon, which in our case is about 0.4M, much lower than the resolution reached by Cattorini et al. (2021), i.e. 0.02M, these differences do not appear particularly significant, especially when we consider that our simulation required only 40 hours on 48 CPUs. Note also that in numerical relativity simulations, the space-time perturbations due to the two BHs close to merging are expected to differ from the analytic metric from Co24, potentially causing inaccuracies in the accretion rate calculation as well. Moreover, the accretion rates in grid-based simulations are typically estimated by integrating the mass flux at the horizon surface, whereas in GIZMO we can directly tag the particles crossing it. We can therefore conclude that our treatment of the binary dynamics is sufficiently robust, and we can proceed by running simulations of binaries embedded in circumbinary discs.

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Mass accretion rate as a function of time in the Bondi-like test. The solution from Cattorini et al. (2021) is shown with the orange line.

3.2. Circumbinary disc evolution

We directly compared a circumbinary disc evolved with the Newtonian version of GIZMO with one evolved with our GR implementation in order to study the effects of the SKS metric implementation on the gas evolution. For both runs, we considered the same initial conditions, which correspond to those used in Franchini et al. (2022, 2023) and produced with the code PHANTOM (Price et al. 2017). We initialised 106 particles, sampling a Newtonian accretion disc extending from Rin = 150M to Rout = 1000M. The initial aspect ratio, H/R, was 0.1. The BH binary was initialised on a circular orbit with a separation d = 100M, and the BHs were non-spinning. We note that the initial configuration was fully Newtonian. This assumption is well justified because the initial binary separation was sufficiently large and the disc was far enough from the BHs to ensure that relativistic corrections played only a minor role. As an equation of state, we used the ideal gas relation, with Γ = 5/3, both in the initial conditions and during the evolution.

We performed the first run using the Newtonian version of the code, where the gas is evolved according to the classical hydrodynamic equation and the binary potential is given by the superposition of the two Keplerian potentials generated by the BH particles. The second run was evolved with the GR scheme, where the gas follows the equations of GR hydrodynamics as presented in Lupi (2023) and with the improvements described in Fedrigo & Lupi (2025), evaluated on the SKS metric. In both cases, we employed the meshless finite-volume (MFM) method and adaptive particle splitting (Franchini et al. 2022), and the BH dynamics was evaluated considering PN correction terms up to the 2.5 order. For simplicity, we decided to evolve both discs as inviscid, thus without either viscosity or magnetic field prescriptions (the latter will be the focus of an upcoming work). Although idealised, this test allowed us to focus solely on the relativistic effects of the binary metric on the gas, without any additional processes that might affect the comparison. The angular momentum transport therefore occurs only via asymmetries in the flow and possible numerical dissipation.

A more comprehensive description of the physics within accretion discs would require solving the radiative transport equations during the hydrodynamic evolution of the system, which would possibly affect the thermodynamics and the structure of the disc (Kaaz et al. 2025). However, such a prescription would be computationally very demanding. A much simpler approach would be to artificially cool the gas to a target entropy, as done in Noble et al. (2012) and Combi et al. (2022), where the evolved discs appeared thinner and with less vertical pressure support (d’Ascoli et al. 2018). In our simulations, however, we decided to neglect any gas cooling prescription to keep the problem as simple as possible. As a consequence, we expected our discs to become warmer and thicker over time. We let the two systems evolve until t = 200T0, where T0 is the orbital period of the binary at t = 0.

In Fig. 2 we show density maps of the two evolved circumbinary discs at t = 200T0. In the Newtonian snapshot (right panel), an empty, slightly elliptical cavity is formed around the binary, carved by the torques produced by the orbiting BHs. The edge of the cavity is found at approximately r ≈ 200 − 250M, i.e. ∼2d, consistent with previous simulations. A faint gas stream can be seen accreting from the disc onto the horizon of one of the BHs.

Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Density map of the central part of the circumbinary disc for the GR run (left) and the Newtonian run (right), at t = 200T0. Density is volume-averaged along the z direction. The white circles mark the BHs’ event horizons, scaled up by a factor of 10 for visualisation purposes.

On the other hand, in the GR evolution (left panel), there is a slightly smaller and more circular internal cavity and a lower density in the inner region of the disc. The inner cavity extends up to r ≈ 150 − 200M, while a second under-dense region can be seen extending up to r ≈ 500M. As expected, both discs become thicker over time, reaching H/R ratios of ∼0.15.

At late times, after ≈100 orbital periods, both discs appear asymmetric, with a m = 1 weak overdensity, which is usually referred to as ‘lump’ (Shi et al. 2012; Noble et al. 2012; Mignon-Risse et al. 2023). This feature is believed to be caused by streams of gas falling into the central cavity and redirected outwards towards the inner edge of the circumbinary disc, and it is observed in both magnetised and non-magnetised circumbinary discs (Noble et al. 2021) around either equal-mass or unequal-mass binaries, provided that the mass ratio is high enough. The gas accumulates in a region that is usually a few times denser than the rest of the inner disc. The origin of the lump is often associated with the Rossby-wave instability (RWI), which can occur in the innermost regions of the circumbinary disc (Mignon-Risse et al. 2023). The RWI is expected to develop in the presence of local maxima in the fluid vortensity (Lovelace et al. 1999). A simpler criterion for the development of the RWI requires the minima of epicyclic frequency (κ) to be min(κ2)≲0.5 − 0.6Ω2 (Chang & Youdin 2024). This condition can be satisfied in the presence of strong density gradients, such as the one between the cavity and the inner edge of the circumbinary disc.

We computed the ratio κ2/Ω2 within our circumbinary disc, finding values of ∼0.5 in the innermost regions, with smaller values limited to a very narrow radial range and only in our GR run. Our density gradient at the cavity edge does not seem strong enough to develop RWI, i.e. we find a ‘lump’ with a density contrast well within a factor of 2. The reason for such a weak ‘lump’ likely lies in the choice of the equation of state, i.e. adiabatic without cooling, and the dimensionality of the problem. For instance, prominent lumps have been found in circumbinary disc simulations that employ an isothermal fluid (Shi et al. 2012) or cooling prescriptions that efficiently dissipate extra heating from shocks (Noble et al. 2012). In the adiabatic case, a prominent lump is found in 2D simulations, which prevents the gas from redistributing along the z direction, boosting density gradients and, as a consequence, more easily triggering the RWI (Mignon-Risse et al. 2023).

For a more quantitative analysis, in Fig. 3 we plot the surface density (Σ) radial profile for both runs at t = 200T0. In the innermost region, the Newtonian evolution yields a surface density that is 2 − 3 times higher than in the GR case, indicating that less gas is entering the central cavity. Between r ≈ 250M and r ≈ 500M, a small drop in the Σ profile highlights the under-dense region discussed above. This feature is more pronounced in the relativistic run, with a drop by a factor of 2. We elaborate on this discrepancy further below in this section.

Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Surface density (Σ) radial profile for the GR (blue), Newtonian (red), and PW (lime) runs at t = 200T0. The initial Σ profile is shown with a dashed black line.

Finally, we evaluated the total mass accretion rate onto the BHs’ event horizons (Ṁ) and show it in Fig. 4. As expected, after an initial transient (grey shaded region), the Newtonian binary (red line) exhibits low gas accretion, with a slow decay over time. We recall that there is no viscosity or magnetic field prescription in the simulations, which are usually the principal processes invoked to transport angular momentum and promote accretion. We note that the adaptive particle splitting procedure can introduce a small diffusivity in the evolution scheme at the radii where the resolution is refined3. However, its adoption in our work enabled us to better resolve the gas morphology within the cavity (Franchini et al. 2022; Duffell et al. 2024). We opted to use the adaptive resolution to correctly resolve the internal structure at the cost of introducing a small diffusivity.

The GR binary (blue line in Fig. 4) shows a non-zero accretion rate characterised by a constant average trend and a quasi-periodic variability related to the inner disc orbital motion. These kinds of oscillations have been seen and studied in depth in many other studies (Paschalidis et al. 2021; Combi et al. 2022; Cattorini et al. 2022; Bright & Paschalidis 2023; Fedrigo et al. 2024; Manikantan et al. 2025). We note that the presence of a non-zero accretion rate in the GR run is consistent with the density drop found in the Σ profile between r ≈ 200 and 500M (Fig. 3). The accreting gas comes from this inner part of the disc, which becomes less dense over time4.

Given the differences between the Newtonian and GR numerical schemes, such a difference might arise due to multiple factors, both numerical or physical. For this reason, we performed tests using more diffusive slope limiters in the Newtonian binary evolution, which resulted in a slightly higher mass accretion rate compared with the reference case. However, even with a higher numerical diffusion, the Newtonian accretion rate remains lower than in the GR run, suggesting that the diffusivity of the scheme cannot explain the different inflow onto the binary in the two cases5.

3.2.1. Post-Newtonian validation test

A physical explanation would be that the higher gas accretion rate is instead due to the non-Keplerian nature of the potential. To test this hypothesis, we performed a run using the Paczynski-Wiita (PW) potential implemented in the Newtonian scheme. The PW potential for a binary can be written as

Φ ( x ) = − G M 1 ( | x − x 1 | − r 1 EH ) − G M 2 ( | x − x 2 | − r 2 EH ) , Mathematical equation: $$ \begin{aligned} \Phi (\mathbf x ) = -\frac{GM_1}{\left(|\mathbf x -\mathbf{x }_1|-\mathbf{r }^\mathrm{EH}_{1}\right)} -\frac{GM_2}{\left(|\mathbf x -\mathbf{x }_2|-\mathbf{r }^\mathrm{EH}_{2}\right)}, \end{aligned} $$(7)

where |x − xi| is the distance between the point where the potential is evaluated and the i-th BH, and riEH is the radius of the i-th BH’s event horizon. This PN potential (Paczyńsky & Wiita 1980) only produces two-thirds of the apsidal precession produced by a relativistic Schwarzschild potential (in the single central object scenario), but it is extremely easy to implement and fast to compute. For this test, we started from the same initial conditions employed in the other two runs and let the system evolve for ∼200 orbits. The resulting Σ profile is shown in Fig. 3, and the Ṁ in Fig. 4. We clearly see that the surface density profile in the PN run is very similar to the GR case. After the usual initial transient, the mass accretion rate displays a trend comparable with the GR case, with stronger oscillation but similar averaged values. In fact, the average Ṁ computed between 50 and 200 orbits is ⟨Ṁ⟩ = 5.4 · 10−15 for the PW run, while it amounts to 6.4 ⋅ 10−15 and 1.2 ⋅ 10−15 in the GR and Newtonian runs, respectively.

Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Mass accretion rate as a function of time for the GR (blue), Newtonian (red), and PW (lime) binaries. The transient phase, caused by the initial data relaxation, is marked with grey shading.

We suggest that relativistic apsidal precession of the gas orbiting a BH binary may play a significant role in the circumbinary disc scenario (Hénon 1965; Paczynski 1977; Rudak & Paczynski 1981). When the trajectories of the orbiting gas are forced by apsidal precession to cross, shocks are generated. During shocks, part of the kinetic energy is converted into thermal energy and angular momentum is transferred outwards, resulting in a net gas inflow towards the centre. To substantiate this idea, we show the density-weighted gas temperature distribution at t = 200T0 in Fig. 5, obtained by measuring the pressure and density from the simulation and then using P = ρkbT/mpμ. As expected, the GR disc tends to be about a factor of 2 − 3 hotter than its Newtonian equivalent in the region affected by shocks, in agreement with our hypothesis. We also note that due to relativistic apsidal precession, the hottest area found at the inner edge of the relativistic disc is shifted forwards in the direction of the disc rotation relative to the Newtonian case. A gas cooling prescription could have mitigated this heating process in both our runs.

Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Density-weighted gas temperature distribution at t = 200T0 for the GR (top) and Newtonian (bottom) runs. The inner cavity (r ≲ 225M) is masked with a white patch. The temperature is normalised for a binary of Mb = 106 M⊙.

The presence of apsidal precession can therefore explain the higher mass accretion rate found in our GR and PW disc runs when compared with the Newtonian simulation. Even though this effect only slightly perturbs the gas trajectory, it can be enough to trigger a faint but non-negligible accretion process. Note that when a more realistic scenario is considered, for example including the presence of an effective viscosity or magnetic fields triggering magneto-rotational instability, the inflow rate we observed will likely be subdominant, thus going unnoticed in the presence of more efficient accretion mechanisms.

3.2.2. Angular momentum transport and effective viscosity

To measure the angular momentum transport within the disc, we computed the effective viscosity (αeff) in our discs in different ways. First, we assumed the prescription by Shakura & Sunyaev (1973), who modelled the viscous stress tensor as simply proportional to the pressure in the disc, α being the constant of proportionality. From this follows the expression for kinematic viscosity, ν ≡ αeffcsH. For a disc in a steady state, the kinematic viscosity is related to the mass accretion rate onto the central object and the disc surface density as ν = Ṁ/3πΣ. We computed the accretion rate radial profile within the disc by measuring the inward flux of material through cylindrical surfaces at different radii. We averaged the results over the last ten binary orbits, and we obtained values of αeff ranging between 10−3 and 10−2. The values of the GR run are 3 − 5 times higher than the Newtonian one. We note that this quantity estimates the global angular momentum transport, taking the contributions from both hydrodynamic stresses and numerical dissipation into account.

Secondly, we computed the effective viscosity using its definition ⟨R ϕr⟩=αeff⟨P⟩, where ⟨R ϕr⟩=⟨ρh δur δuϕ⟩ is the Reynolds stress tensor, δuμ = uμ − ⟨ρuμ⟩/⟨ρ⟩ is the velocity fluctuation relative to the mean motion, and ⟨ ⋅ ⟩ is a time and azimuthal average operator. We used the relativistic formalism presented in Noble et al. (2012), and we averaged again over the last ten binary orbits. Using this prescription, we obtained effective viscosity values of up to 5 × 10−3 and 2 × 10−3 in the outermost part of both discs for the GR and Newtonian run, respectively. A somewhat irregular behaviour is instead seen in the inner part of the disc, between r ≈ 2.5d and r ≈ 5d, where the development of spiral waves and asymmetries in the gas distribution results in a noisier local effective viscosity, with values reaching 10−2 (Mignon-Risse et al. 2025). At the cavity edge, αeff reaches a maximum value of ∼2.5 × 10−2. As suggested in Noble et al. (2012), we also computed the Reynolds stress tensor as ⟨R ϕr⟩=⟨T ϕr⟩−⟨A ϕr⟩, where T is the stress-energy tensor, ⟨A ϕr⟩=⟨ρl⟩⟨ρhur⟩/⟨ρ⟩, and l = −uϕ/ut is specific angular momentum. As expected, this equivalent formulation yielded results compatible with the αeff estimation via the fluid velocity fluctuations.

Comparing the different αeff estimations, we can conclude that the angular momentum transport in our discs is dominated by the hydrodynamic stresses from fluid asymmetries. However, this effect yields an effective viscosity 1 − 2 orders of magnitude lower than values usually adopted in viscous disc simulations, for example α ∼ 0.1 (Franchini et al. 2024a).

3.2.3. Minidisc morphology and dynamics

As a final analysis, we considered the morphology of the accreting structures that form close to the BHs’ horizons, usually called ‘minidiscs’, which are expected to be affected by relativistic corrections. To define these structures, we considered the gas present within each BH’s Roche lobe (Eggleton 1983). In the case of an equal-mass binary, this region can be approximated as a sphere centred on the BH, with an equivalent-volume radius of Req ≃ 0.38d (Leahy & Leahy 2015). We note that the very low angular momentum transport and the consequently low density in the central cavity lead to transient and highly variable ‘minidisc’ structures (Fig. 6) that are hardly comparable in morphology to those found in simulations with stronger angular momentum transport (e.g. GRMHD simulations in Noble et al. 2012; Combi et al. 2022 or viscous simulations in Franchini et al. 2022). However, they serve our purpose to study GR effects caused by the metric.

Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Density maps of the minidiscs during a filling phase (after ≃60 BH orbits from the start of the run) for the GR (top) and Newtonian (bottom) runs. The black dots mark the event horizons, the solid lines the innermost stable circular orbit radii, and the dashed lines the Roche lobe radii. Snapshots are aligned with the binary axis.

We compared the time evolution of minidiscs in the two runs, finding a more uniform and homogeneous gas distribution in the relativistic case, and a more asymmetric distribution in the Newtonian one. Furthermore, the latter is more prone to strong density fluctuations, indicating moments where the Roche lobes are almost empty. The more stable evolution and the more homogeneous gas distribution of the GR minidiscs can be understood in terms of the higher accretion rate in the cavity. Relativistic apsidal precession plays an important role here too. In our simulation, the gas streams falling towards the relativistic BHs circularise, producing more homogeneous gas structures around the BHs and effectively mitigating density asymmetries in the minidiscs.

We evaluated the minidisc eccentricity in our simulations from the gas energy within the Roche lobes. We performed this analysis between t = 50T0 and t = 150T0, i.e. excluding the initial transient and the final part of the simulations in which the Newtonian cavity is the emptiest, and we plot the results in Fig. 7. While the Newtonian minidiscs show an eccentricity that oscillates around e ≈ 0.6, the GR and PW ones average approximately to e ≈ 0.5. Hence, the relativistic minidiscs tend to be slightly more circular than the Newtonian ones. We note that the circularisation effect is small because we are considering non-spinning BHs, whose apsidal precession effect is less strong when compared with spinning BHs.

Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Average gas eccentricity within the minidiscs as a function of time for the GR (blue), Newtonian (red), and PW (lime) runs.

Finally, we note that, while the PW prescription can produce results similar to our GR formulation in the scenario we considered here, it cannot emulate the spin of the BH and its effects on the gas dynamics. Also, the PW potential is divergent at the event horizon, which can lead to numerical instabilities and artefacts.

To conclude, our results highlight the importance of considering GR effects in the evolution of circumbinary discs and BHs’ minidiscs, as they play an important role even at distances as large as hundreds of gravitational radii and thousands of orbits before the binary merger.

4. Conclusions

In this paper we have presented our implementation of the SKS metric from Co24 in the relativistic scheme of the GIZMO code. The code can now simulate BH binaries at intermediate separations (101M ≲ d ≲ 103M) with great computational efficiency and accuracy, taking relativistic effects in the gas and magnetic field evolution into account.

We tested our implementation by comparing our results with a numerical relativity simulation of a BH binary embedded in a uniform gas distribution, obtaining good agreement in the gas evolution behaviour and in the estimated mass accretion rate onto the BH horizons. Then, we compared the evolution of two circumbinary discs evolved using the Newtonian and relativistic schemes in GIZMO. Both simulations were stable and produced reasonable results. The GR run, however, yielded a higher rate of gas accreted onto the horizons and a lower-density region in the innermost part of the disc, between r ≈ 2.5d and r ≈ 5d. We find that the presence of a non-Keplerian potential is the cause of this phenomenon, as demonstrated by a PN simulation performed with the standard version of GIZMO.

With this work, we highlight the importance of including relativistic effects in simulations of circumbinary discs around BH binaries approaching a merger, even at scales not yet accessible with numerical relativity simulations. In future work, we plan to explore the parameter space of long-evolved magnetised circumbinary discs, in order to search for characteristic features that could suggest EM signatures.

Acknowledgments

We thank the anonymous referee for very useful comments that improved the quality of the manuscript. We thank Daniel Price and Enrico Ragusa for useful discussions. AL acknowledges support from PRIN MUR ‘2022935STW’. GF thanks Sormani Mattia Carlo and Zoltán Haiman for their fruitful suggestions on the angular momentum transport interpretation. MB acknowledges support from the Italian Ministry for Universities and Research (MUR) programme ‘Dipartimenti di Eccellenza 2023–2027’, within the framework of the activities of the Centro Bicocca di Cosmologia Quantitativa (BiCoQ). AF acknowledges financial support from the Unione europea – Next Generation EU, Missione 4 Componente 1 CUP G43C24002290001.

References

  1. Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023, ApJ, 951, L8 [NASA ADS] [CrossRef] [Google Scholar]
  2. Amaro-Seoane, P., Aoudia, S., Babak, S., et al. 2012, ArXiv e-prints [arXiv:1201.3621] [Google Scholar]
  3. Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, ArXiv e-prints [arXiv:1702.00786] [Google Scholar]
  4. Avara, M. J., Krolik, J. H., Campanelli, M., et al. 2024, ApJ, 974, 242 [Google Scholar]
  5. Blanchet, L. 2014, Liv. Rev. Relat., 17, 2 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bode, T., Haas, R., Bogdanović, T., Laguna, P., & Shoemaker, D. 2010, ApJ, 715, 1117 [Google Scholar]
  7. Bogdanović, T., Bode, T., Haas, R., Laguna, P., & Shoemaker, D. 2011, Class. Quant. Grav., 28, 094020 [CrossRef] [Google Scholar]
  8. Bogdanović, T., Miller, M. C., & Blecha, L. 2022, Liv. Rev. Relat., 25, 3 [Google Scholar]
  9. Bowen, D. B., Mewes, V., Campanelli, M., et al. 2018, ApJ, 853, L17 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bright, J. C., & Paschalidis, V. 2023, MNRAS, 520, 392 [Google Scholar]
  11. Capelo, P. R., Volonteri, M., Dotti, M., et al. 2015, MNRAS, 447, 2123 [NASA ADS] [CrossRef] [Google Scholar]
  12. Cattorini, F., & Giacomazzo, B. 2024, Astropart. Phys., 154, 102892 [Google Scholar]
  13. Cattorini, F., Giacomazzo, B., Haardt, F., & Colpi, M. 2021, Phys. Rev. D, 103, 103022 [NASA ADS] [CrossRef] [Google Scholar]
  14. Cattorini, F., Maggioni, S., Giacomazzo, B., et al. 2022, ApJ, 930, L1 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chang, E., & Youdin, A. N. 2024, ApJ, 976, 100 [Google Scholar]
  16. Chiesa, A., Izquierdo-Villalba, D., Sesana, A., et al. 2025, A&A, 707, A89 [Google Scholar]
  17. Cocchiararo, F., Franchini, A., Lupi, A., & Sesana, A. 2024, A&A, 691, A250 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Cocchiararo, F., Franchini, A., Lupi, A., & Sesana, A. 2025, ArXiv e-prints [arXiv:2510.21919] [Google Scholar]
  19. Combi, L., & Ressler, S. 2024, https://doi.org/10.5281/zenodo.10841021 [Google Scholar]
  20. Combi, L., & Ressler, S. M. 2026, Phys. Rev. D, 113, 044023 [Google Scholar]
  21. Combi, L., Lopez Armengol, F. G., Campanelli, M., et al. 2021, Phys. Rev. D, 104, 044041 [NASA ADS] [CrossRef] [Google Scholar]
  22. Combi, L., Lopez Armengol, F. G., Campanelli, M., et al. 2022, ApJ, 928, 187 [NASA ADS] [CrossRef] [Google Scholar]
  23. Csizmadia, P., Debreczeni, G., Rácz, I., & Vasúth, M. 2012, Class. Quant. Grav., 29, 245002 [Google Scholar]
  24. d’Ascoli, S., Noble, S. C., Bowen, D. B., et al. 2018, ApJ, 865, 140 [CrossRef] [Google Scholar]
  25. Dotti, M., Colpi, M., & Haardt, F. 2006, MNRAS, 367, 103 [Google Scholar]
  26. Dotti, M., Sesana, A., & Decarli, R. 2012, Adv. Astron., 2012, 1 [Google Scholar]
  27. Dotti, M., Bonetti, M., Rigamonti, F., et al. 2023, MNRAS, 518, 4172 [Google Scholar]
  28. Duffell, P. C., Dittmann, A. J., D’Orazio, D. J., et al. 2024, ApJ, 970, 156 [NASA ADS] [CrossRef] [Google Scholar]
  29. Eggleton, P. P. 1983, ApJ, 268, 368 [Google Scholar]
  30. Ennoggi, L., Campanelli, M., Zlochower, Y., et al. 2025a, Phys. Rev. D, 112, 063009 [Google Scholar]
  31. Ennoggi, L., Campanelli, M., Krolik, J., et al. 2025b, ArXiv e-prints [arXiv:2509.10319] [Google Scholar]
  32. EPTA Collaboration, InPTA Collaboration, Antoniadis, J., et al. 2023, A&A, 678, A50 [CrossRef] [EDP Sciences] [Google Scholar]
  33. Fedrigo, G., & Lupi, A. 2025, A&A, 700, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Fedrigo, G., Cattorini, F., Giacomazzo, B., & Colpi, M. 2024, Phys. Rev. D, 109, 103024 [CrossRef] [Google Scholar]
  35. Franchini, A., Sesana, A., & Dotti, M. 2021, MNRAS, 507, 1458 [NASA ADS] [Google Scholar]
  36. Franchini, A., Lupi, A., & Sesana, A. 2022, ApJ, 929, L13 [NASA ADS] [CrossRef] [Google Scholar]
  37. Franchini, A., Lupi, A., Sesana, A., & Haiman, Z. 2023, MNRAS, 522, 1569 [NASA ADS] [CrossRef] [Google Scholar]
  38. Franchini, A., Bonetti, M., Lupi, A., & Sesana, A. 2024a, A&A, 686, A288 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Franchini, A., Prato, A., Longarini, C., & Sesana, A. 2024b, A&A, 688, A174 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Giacomazzo, B., Baker, J. G., Miller, M. C., Reynolds, C. S., & van Meter, J. R. 2012, ApJ, 752, L15 [Google Scholar]
  41. Gold, R. 2019, Galaxies, 7, 63 [Google Scholar]
  42. Gold, R., Paschalidis, V., Etienne, Z. B., Shapiro, S. L., & Pfeiffer, H. P. 2014, Phys. Rev. D, 89, 064060 [NASA ADS] [CrossRef] [Google Scholar]
  43. Gutiérrez, E. M., Combi, L., Noble, S. C., et al. 2022, ApJ, 928, 137 [CrossRef] [Google Scholar]
  44. Hénon, M. 1965, Ann. Astrophys., 28, 992 [Google Scholar]
  45. Hopkins, P. F. 2015, MNRAS, 450, 53 [Google Scholar]
  46. Kaaz, N., Liska, M., Tchekhovskoy, A., Hopkins, P. F., & Jacquemin-Ide, J. 2025, ApJ, 979, 248 [Google Scholar]
  47. Kocsis, B., Frei, Z., Haiman, Z., & Menou, K. 2006, ApJ, 637, 27 [NASA ADS] [CrossRef] [Google Scholar]
  48. Krauth, L. M., Davelaar, J., Haiman, Z., et al. 2023, MNRAS, 526, 5441 [NASA ADS] [CrossRef] [Google Scholar]
  49. Leahy, D. A., & Leahy, J. C. 2015, Comput. Astrophys. Cosmol., 2, 4 [NASA ADS] [CrossRef] [Google Scholar]
  50. Lopez Armengol, F. G., Combi, L., Campanelli, M., et al. 2021, ApJ, 913, 16 [NASA ADS] [CrossRef] [Google Scholar]
  51. Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] [Google Scholar]
  52. Lupi, A. 2023, MNRAS, 519, 1115 [Google Scholar]
  53. MacFadyen, A. I., & Milosavljević, M. 2008, ApJ, 672, 83 [Google Scholar]
  54. Manikantan, V., Paschalidis, V., & Bozzola, G. 2025, ApJ, 984, L47 [Google Scholar]
  55. Mayer, L., Kazantzidis, S., Mastropietro, C., & Wadsley, J. 2007, Nature, 445, 738 [NASA ADS] [CrossRef] [Google Scholar]
  56. Mignon-Risse, R., Varniere, P., & Casse, F. 2023, MNRAS, 520, 1285 [NASA ADS] [CrossRef] [Google Scholar]
  57. Mignon-Risse, R., Varniere, P., & Casse, F. 2025, MNRAS, 544, 4226 [Google Scholar]
  58. Muñoz, D. J., Lai, D., Kratter, K., & Miranda, R. 2020, ApJ, 889, 114 [Google Scholar]
  59. Noble, S. C., Mundim, B. C., Nakano, H., et al. 2012, ApJ, 755, 51 [NASA ADS] [CrossRef] [Google Scholar]
  60. Noble, S. C., Krolik, J. H., Campanelli, M., et al. 2021, ApJ, 922, 175 [NASA ADS] [CrossRef] [Google Scholar]
  61. Paczynski, B. 1977, ApJ, 216, 822 [CrossRef] [Google Scholar]
  62. Paczyńsky, B., & Wiita, P. J. 1980, A&A, 88, 23 [NASA ADS] [Google Scholar]
  63. Paschalidis, V., Bright, J., Ruiz, M., & Gold, R. 2021, ApJ, 910, L26 [Google Scholar]
  64. Price, D. J., Wurster, J., Nixon, C., et al. 2017, Astrophysics Source Code Library [record ascl:1709.002] [Google Scholar]
  65. Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6 [NASA ADS] [CrossRef] [Google Scholar]
  66. Rigamonti, F., Severgnini, P., Sottocorno, E., et al. 2025, A&A, 693, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Rudak, B., & Paczynski, B. 1981, Acta Astron., 31, 13 [NASA ADS] [Google Scholar]
  68. Ruiz, M., Tsokaros, A., & Shapiro, S. L. 2023, Phys. Rev. D, 108, 124043 [Google Scholar]
  69. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  70. Shi, J., Krolik, J. H., Lubow, S. H., & Hawley, J. F. 2012, ApJ, 749, 118 [Google Scholar]
  71. Xu, H., Chen, S., Guo, Y., et al. 2023, RAA, 23, 075024 [NASA ADS] [Google Scholar]
1

Note that the code can also work with pre-computed trajectories such as those provided by the CBWAVES tool (Csizmadia et al. 2012).

2

As an internal-consistency check, we also verified that the SKS metric recovers the single Kerr metric when one of the two BH masses is set to zero and the massive BH is set to rest, by integrating circular orbits of collisionless particles.

3

During the testing, we found that the splitting procedure in the GR scheme is slightly more diffusive than the Newtonian one. However, in our runs, this difference results in a Ṁ increase of 10−2⟨Ṁ⟩.

4

We let the relativistic run evolve until t = 600T0; a larger cavity with radius r ≃ 500M is carved around the BHs.

5

We also tested whether the initial conditions had any effect on the late-time dynamics of the gas in the disc, finding none.

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Mass accretion rate as a function of time in the Bondi-like test. The solution from Cattorini et al. (2021) is shown with the orange line.
In the text
Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Density map of the central part of the circumbinary disc for the GR run (left) and the Newtonian run (right), at t = 200T0. Density is volume-averaged along the z direction. The white circles mark the BHs’ event horizons, scaled up by a factor of 10 for visualisation purposes.
In the text
Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Surface density (Σ) radial profile for the GR (blue), Newtonian (red), and PW (lime) runs at t = 200T0. The initial Σ profile is shown with a dashed black line.
In the text
Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Mass accretion rate as a function of time for the GR (blue), Newtonian (red), and PW (lime) binaries. The transient phase, caused by the initial data relaxation, is marked with grey shading.
In the text
Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Density-weighted gas temperature distribution at t = 200T0 for the GR (top) and Newtonian (bottom) runs. The inner cavity (r ≲ 225M) is masked with a white patch. The temperature is normalised for a binary of Mb = 106 M⊙.
In the text
Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

Density maps of the minidiscs during a filling phase (after ≃60 BH orbits from the start of the run) for the GR (top) and Newtonian (bottom) runs. The black dots mark the event horizons, the solid lines the innermost stable circular orbit radii, and the dashed lines the Roche lobe radii. Snapshots are aligned with the binary axis.
In the text
Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Average gas eccentricity within the minidiscs as a function of time for the GR (blue), Newtonian (red), and PW (lime) runs.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.