EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 691 (November 2024) A&A, 691 (2024) A250 Full HTML
Open Access
Issue
A&A
Volume 691, November 2024
Article Number A250
Number of page(s) 15
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202449598
Published online 18 November 2024

© The Authors 2024

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. Subscribe to A&A to support open access publication.

1. Introduction

The existence of a massive black hole (MBH) at the centre of most, if not all, galaxies is a well-established observational fact (Kormendy & Ho 2013, and references therein). In the hierarchical clustering framework of structure formation, galaxies grow by accreting material from the intergalactic medium and by merging with other galaxies. Galaxy mergers are expected to result in the formation of MBH binaries (MBHBs). Indeed, during these merging events, the MBHs hosted at the centre of the parent galaxies migrate towards the centre of the merger remnant due to the dynamical friction mechanism (Chandrasekhar 1943). At parsec separation, where the two MBHs bind together to form a binary (Begelman et al. 1980), dynamical friction becomes inefficient in bringing the binary components closer together (Mayer et al. 2007; Dotti et al. 2012). Since gravitational waves (GWs) are efficient only at binary separations below ∼10−3 pc, it is necessary to invoke one or more astrophysical mechanisms to shrink the binary further down, to the point where GWs can efficiently drive the merger. The most explored mechanisms involve the interaction of the binary with either its stellar (Quinlan 1996; Khan et al. 2011; Preto et al. 2011) or gaseous (Armitage & Natarajan 2002; Escala et al. 2005; Lodato et al. 2009; Cuadra et al. 2009) environment. While the evolution driven by stars is always found to shrink the binary due to the energy and angular momentum exchange with stars in three-body stellar encounters (e.g. Sesana et al. 2007; Bortolas et al. 2021), the effect of the interaction with the gaseous environment has recently been widely debated. The effect of the circumbinary disc has been investigated extensively using 2D (Muñoz et al. 2019; Duffell et al. 2020; Tiede et al. 2020; Siwek et al. 2023) and 3D (Heath & Nixon 2020; Franchini et al. 2021, 2022, 2023) hydrodynamical simulations. The general consensus is that, although relatively warm discs with an aspect ratio of around 0.1 have been found to expand the binary, shrinking is usually promoted by cold, less viscous circumbinary discs.

As a result of the interaction with the gaseous environment in which they reside, MBHBs are expected to produce distinctive observational electromagnetic (EM) signatures. Currently known candidates at sub-parsec separations have been inferred from their EM emission (Bogdanović et al. 2022) through either photometric measurements of quasi-periodic variability both in the optical (Graham et al. 2015; Charisi et al. 2016; Liu et al. 2019; Chen et al. 2020) and γ rays (e.g. Sandrinelli et al. 2018; Peñil et al. 2022) or spectroscopic measurements of offset broad emission lines (Bogdanović et al. 2009; Dotti et al. 2009; Boroson & Lauer 2009; Tang & Grindlay 2009; Decarli et al. 2010; Barrows et al. 2011; Tsalmantza et al. 2011; Eracleous et al. 2012; Tsai et al. 2013). Convincing evidence that these sources are indeed MBHBs is still missing, both due to the lack of a firm theoretical understanding of the expected distinctive emission signatures of these systems and because the observed features can be accommodated by alternative explanations (see e.g. Eracleous & Halpern 1994; Vaughan et al. 2016).

These observationally elusive MBHBs are also the main targets of current and future GW experiments. The future space-based Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2023), to be launched in 2035, will survey the milli-hertz GW window, where it is expected to observe the late inspiral and mergers of MBHBs with masses in the range of 104 − 107 M⊙ everywhere in the Universe. At even lower frequencies, in the nano-hertz regime, precision timing of an ensemble of millisecond pulsars within a so-called pulsar timing array (PTA, Foster & Backer 1990) is sensitive to GWs emitted by MBHBs of billions of solar masses at milli-parsec separations (see e.g. Sesana et al. 2008). Although GWs might be eventually needed for a firm detection of MBHBs, the identification of convincing EM counterparts will help with pinning down the properties of the emitting systems, and thus realising the promises of low-frequency multi-messenger astronomy. The importance of understanding the EM emission of MBHBs has became even more compelling following the recent detection of a signal compatible with a GW background (GWB) reported by the European and Indian PTA (EPTA+InPTA, EPTA Collaboration 2023; EPTA Collaboration & InPTA Collaboration 2023a,b, 2024a; EPTA Collaboration & InPTA Collaboration 2024b; Smarra et al. 2023), NANOGrav (Agazie et al. 2023a,b,c; Afzal et al. 2023), Parkes PTA (PPTA, Reardon et al. 2023), and Chinese PTA (CPTA, Xu et al. 2023). Although the properties of the signal are compatible with a cosmic population of MBHBs, its origin cannot yet be observationally established from the GW data alone, and the identification of EM counterparts might aid in confirming the origin of this signal.

Attempts to characterise the EM signatures from MBHBs at very small separations (i.e. tens to a hundred gravitational radii) have been made by several authors (e.g. Tanaka et al. 2012; Tang et al. 2018; d’Ascoli et al. 2018; Gutiérrez et al. 2022; Combi et al. 2022; Krauth et al. 2023; Franchini et al. 2024). The same regime of separations has recently been investigated in a series of 2D numerical simulations, indicating that the lump periodicity coming from the streams and the cavity wall is more noticeable in circular binaries being a possible smoking gun signature of circularity (Westernacher-Schneider et al. 2022). Moreover, Westernacher-Schneider et al. (2024) proposes that EM signatures produced by MBHBs could be different from single MBH quasi-periodic sources as a consequence of the formation of eccentric mini-discs, prompted by periodic mass exchange between mini-discs. The characterisation of the EM signatures at intermediate scales – that is, separations ∼10−4 − 10−2 pc from full 3D hydrodynamical simulations that resolve the thermal evolution of the circumbinary disc – is currently missing and is of fundamental importance for aiding the identification of the origin of the GWB signal found with PTA experiments. Furthermore, the modelling of these signatures is of paramount importance for the identification of MBHB candidates in time-domain surveys, which can constitute the cosmic population of precursors of the merging binaries that will be detected by LISA.

In this work, we have computed for the first time the spectral energy distributions (SEDs) and the multi-wavelength light curves (LCs) from 3D numerical simulations with hyper-Lagrangian refinement of milli-parsec scale binaries. We have taken into account the thermodynamic evolution of the gas using a radiative cooling prescription in the form of black-body radiation. We have also included the disc’s self-gravity, which is usually neglected, and the Shakura-Sunyaev prescription for viscosity (Shakura & Sunyaev 1973). We modelled binaries with an eccentricity of e = 0,  0.9 and mass ratio of q = 0.1,  1 and then computed their emitted spectra and LCs. We also added to the spectrum the contribution from a non-thermal corona, based on the assumption that its emission is proportional to the radiation emitted by the discs around each binary component (also called mini-discs). We placed the simulated binary at different redshifts of z =  0.1,  0.4,  0.7, and analysed the observed flux in different frequency bands, mainly focussing on the optical band that will be probed by the upcoming Vera Rubin Observatory (VRO; LSST Science Collaboration 2009). As an exercise, we then assumed the source to be two orders of magnitude brighter, but still preserving the same LC, and we investigated the observability with the VRO at higher redshifts of z =  1,  2,  3. Finally, we analysed the time evolution of the disc and binary properties such as the aspect ratio, H/R (where H is the disc vertical height and R is the distance from the central MBH in the disc plane), the semi-major axis, a0, and the eccentricity, e, of the binaries.

The paper is organised as follows. In Sect. 2, we describe the numerical details of the simulations and the physical parameters we use, the circumbinary disc physics assumptions, and the thermal emission calculation. 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 draw our conclusions and discuss possible observational implications in Sect. 4.

2. Physical and numerical set-up

2.1. Binary and circumbinary disc model

We performed 3D hyper-Lagrangian resolution hydrodynamics simulations of the evolution of the thermal circumbinary disc around an accreting MBHB. We used the 3D meshless finite mass (MFM) version of the code GIZMO (Hopkins 2015) combined with the adaptive particle splitting technique (see Franchini et al. 2022, for details) to increase the resolution inside the cavity carved by the binary. We ran four simulations combining two values of the binary mass ratio, q = 0.1,  1, and eccentricity, e = 0,  0.9, each spanning over 1300 binary orbital periods. The only exception is the circular binary with a mass ratio of q = 0.1, which we evolved for just over 1060 orbital periods since the disc reached a steady state after fewer orbits in this case. The MBHBs have, in code units, a total mass of MB = M1 + M2 = 1, where M1 and M2 are the masses of the primary and secondary black hole, respectively, and the initial semi-major axis, a0 = 1. Each MBH is modelled as a sink particle with an accretion radius of Rsink = 0.05a. As the binary orbit is allowed to evolve with time (Franchini et al. 2023), conservation of mass, linear and angular momentum are ensured during each accretion event, in the same way as is done in the PHANTOM code (Bate et al. 1995). We initially modelled the circumbinary disc with N = 106 gas particles for a total mass of MD = 0.01 MB. The mass is distributed with an initial surface density profile of Σ ∝ R−1 and a radial extent between Rin = 2a and Rout = 10a. The disc is co-planar with the binary orbit and has an initial aspect ratio of H/R = 0.1 in all our simulations. We generated the disc initial conditions using the SPH code PHANTOM (Price et al. 2017).

In our model, the thermodynamic evolution of the gas follows an adiabatic equation of state with index γ = 5/3, so that the disc is allowed to heat and cool and we can capture the effect of shocks. In the Shakura-Sunyaev accretion disc model, the angular momentum transport within the disc is modelled with the α viscosity parameter, which constitutes a simple parameterisation of the turbulence within the disc. This turbulence might be driven by the magneto-rotational instability in highly ionised discs (Balbus & Hawley 1991). We therefore included viscosity using the Shakura-Sunyaev parameterisation with a kinematic viscosity of ν = αcsH, where cs is the gas sound speed, and set α = 0.1 in all our simulations. Since we do also include the disc self-gravity together with a cooling prescription, our discs may develop gravitational instabilities that will eventually result in an additional source of angular momentum transport (Lodato & Rice 2004, 2005; Lodato 2007). To ensure the gravitational stability of the disc, we set the initial Toomre parameter, Q >  1 (Toomre 1964).

The vast majority of existing numerical simulations that consider the disc self-gravity (Gammie 2001; Lodato & Rice 2005; Cuadra et al. 2009; Roedig & Sesana 2014; Franchini et al. 2021) make use of a simple cooling prescription that assumes gas to cool on a timescale, tcool, that is a multiple, βcool, of the dynamical time, tdyn = Ω−1; that is, βcool = Ωtcool. In this work, we have used a more realistic radiative cooling in the form of black-body radiation. The cooling rate is given by

Q ˙ = 8 3 σ SB T i 4 κ Σ , $$ \begin{aligned} \dot{Q} = \frac{8}{3} \frac{\sigma _{\rm SB} T_{i}^4}{\kappa \Sigma }, \end{aligned} $$(1)

where σSB is the Stephan-Boltzmann constant, κ is the opacity, Σ is the disc surface density, and Ti is the temperature of each element. We assume the opacity, κ, to be a combination of the Kramer opacity, κ Kramer ∝ ρ T i − 7 / 2 $ \kappa_{\mathrm{Kramer}} \propto \rho T_{i}^{-7/2} $, and the electron-scattering opacity, κes = 0.2(1 + X) g cm2, with a hydrogen mass fraction of X = 0.59. When scaled to physical units, all our simulations have a binary total mass of MB = 106 M⊙ and initial separation of a = 4.8 ⋅ 10−4 pc ≃ 1.2 ⋅ 104 Rg, where Rg = GMB/c2 is the gravitational radius of the binary. This implies that we can follow the gas only down to Rsink = 0.05a ∼ 600 Rg ∼ 100 RISCO. Therefore, we should bear in mind that in the following we are neglecting the emission in the region between 100 RISCO and RISCO and the total luminosity of the system is likely going to be underestimated, especially in the UV and X-ray bands. Choosing this particular semi-major axis value sets the orbital period of the binaries to 1 year. According to the observing time of a VRO survey of 10 years, our choice of semi-major axis implies that we can observe the flux emitted by the source over 10 orbital periods. In order to investigate the thermodynamics that governs our circumbinary discs, we measured the rate of change in the internal energy due to both heating (i.e. shocks, PdV work, and viscosity) and cooling processes, and we then compared it to the cooling rate in Eq. (1). We defined an effective βeff as

β eff = u u ˙ Ω , $$ \begin{aligned} \beta _{\rm eff} = \frac{u}{\dot{u}} \Omega , \end{aligned} $$(2)

where u is the particle internal energy. We find βeff to depend on the radius and to reach a stable constant value of βeff ∼ 5 at r = 2a and βeff ∼ 1.7 × 103 at r = 10a after an initial transient phase of about ∼800 PB in all our simulations. We compared βeff with the βcool we measured from the black-body cooling and we found that βcool ≪ βeff. We can therefore conclude that the heating mechanisms ensure the disc stability against the black-body cooling, which initially causes the transient phase.

2.2. Emission model

For each gas particle, i, we calculated the temperature, Ti, assuming that both gas and radiation pressure contribute to the hydrostatic equilibrium of the disc, and thus numerically solving the following implicit equation for Ti at each resolution element:

P tot = P gas + P rad = ρ k B T i m p μ + 4 3 σ SB T i 4 c · $$ \begin{aligned} P_{\rm tot} = P_{\rm gas} + P_{\rm rad} = \frac{\rho k_{\rm B} T_{i}}{m_{\rm p}\mu } + \frac{4}{3}\frac{\sigma _{\rm SB} T_{i}^4}{c}\cdot \end{aligned} $$(3)

We divided the disc temperature domain into a 2D matrix with dimensions of 512  ×  512 pixels in the x − y plane, which coincides with the MBHB orbital plane. For each pixel, we computed the midplane temperature as the average temperature of all the particles within the z co-ordinate −0.05a <  z <  0.05a, obtaining the midplane temperature matrix, T. For a more detailed analysis of the mini-discs, we further divided the spatial domain from the sink radius of each MBH out to r = 3a into a 2D matrix of 512 × 512 pixels. For each pixel in this grid, we computed the midplane temperature, T, using the same method outlined above. We then computed the effective temperature in optically thick approximation (κΣ >  1) in each element of both matrices as

T eff 4 = 4 3 T 4 κ Σ , $$ \begin{aligned} T_{\rm eff}^4 = \frac{4}{3}\frac{T^4}{\kappa \Sigma }, \end{aligned} $$(4)

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

We obtained the flux emitted by each pixel using Planck’s formula:

d L ν ≡ B ν d ν d A = 2 h ν 3 c 2 d ν exp ( h ν k B T eff ) − 1 d A , $$ \begin{aligned} \mathrm{d}L_{\nu } \equiv B_{\nu } \, \mathrm{d}\nu \mathrm{d}A = \frac{2h \nu ^3}{c^2}\, \frac{\mathrm{d}\nu }{\exp \left({\frac{h\nu }{k_{\rm B}T_{\rm eff}}}\right)-1} \mathrm{d}A, \end{aligned} $$(5)

where h is the Planck constant, kB is the Boltzmann constant, and dA is the area of each element. In order to analyse the different contributions to the emission from each part of the disc, we have divided the simulated domain into five different regions: the two mini-discs that extend from the sink radius of each component to the Roche Lobe size, the streams region that extends from outside the Roche Lobe out to r = 3a, an inner part of the disc, 3a <  r <  5a, and an outer part of the disc, 5a <  r <  10a. For each region, we computed the SED through the sum of each pixel flux obtained with Eq. (5).

We further added the non-thermal emission that is expected to originate from the corona using the correction to the bolometric luminosity in Duras et al. (2020) in the hard X-ray (i.e. 0.2 − 10 keV) band. We assume that the corona emission is proportional to the total luminosity of the mini-discs through the correction factor, Kband, which in this case is defined as the ratio between the mini-discs’ luminosity, LMDs, and the luminosity in a given spectral band, Lband; that is, Kband = LMDs/Lband. We assume the mini-discs to behave similarly to single MBH AGN discs, where the X-ray emission is produced by the corona, as there is no clear description available to date of the properties of the non-thermal corona in the accretion discs that form around the binary components. Following Duras et al. (2020), the X-ray correction factor, KX, is

K X ( L MDs ) = a [ 1 + ( log ( L MDs / L ⊙ ) b ) c ] , $$ \begin{aligned} K_{\rm X}(L_{\rm MDs}) = a \left[1 + \left(\frac{\log (L_{\rm MDs}/L_{\odot })}{b}\right)^c\right], \end{aligned} $$(6)

where a, b, c are best-fit parameters shown in Table 1 in Duras et al. (2020). Moreover, we assumed that the hot corona emission produces a power law spectrum with νLν ∝ ν−0.7 (Regan et al. 2019).

Since we ultimately want to obtain LCs to compare with future observations, we placed the source at different distances from the observer and, assuming isotropic emission, we computed the total flux in different bands integrating the observed flux, fν(νo), over a specific range of frequencies:

F = ∫ ν 1 ν 2 f ν ( ν o ) d ν o = ∫ ν 1 ν 2 l ν ( ν e ) 4 π d L 2 ( 1 + z ) d ν o , $$ \begin{aligned} F = \int _{\nu _1}^{\nu _2} f_{\nu }(\nu _{\rm o}) \, d \nu _{\rm o} = \int _{\nu _1}^{\nu _2} \frac{l_{\nu }(\nu _e)}{4\pi d_{\rm L}^2}(1+z) \, d \nu _{\rm o}, \end{aligned} $$(7)

where lν(νe) is the luminosity as a function of the frequency, νe, and dL is the luminosity distance (Hogg 1999; Hogg et al. 2002).

We then analysed the changes in the observed flux over time in order to identify modulations that might signal the presence of a MBHB.

In order to understand whether these changes in the observed flux are detectable, we need to simulate a realistic detection scenario, taking into account the sensitivity limit of the employed survey and its observation time span. In particular, in this study we focus on MBHB observability with the VRO. To infer the observability of the systems, we adopted a similar method to the one used by Kelley et al. (2019). We added to the simulated fluxes a Gaussian noise with variance given by the telescope’s 5σ sensitivity in the considered band. We then performed a fast Fourier transform (FFT) onto a limited number of cycles (commensurate to a survey time span of ≈10 yr) and qualitatively identified whether prominent peaks in the spectrum appear and what their origin is. We defer a more thorough statistical study of the detectability of these features to a follow-up study. The VRO 5σ sensitivities in different bands are listed in Table 1.

Table 1.

Frequency bands and 5-σ flux sensitivities of the VRO telescope’s optical filters.

3. Results

We here present the results that we obtained from our numerical simulations. In particular, we show the time evolution of the disc aspect ratio, H/R, and briefly discuss the evolution of the binary orbital eccentricity and semi-major axis. We then focus on the disc emission and on the LCs in various frequency bands and on the effect that the explored parameters (i.e. mass ratio, eccentricity, and redshift) have on the detectability of the EM emission.

3.1. Disc and binary evolution

Since our numerical simulations are 3D and the gas is allowed to change its temperature with time, we can resolve gas shocks and investigate the effect they have on the final temperature of the disc. As showed by previous simulations in the literature (Artymowicz & Lubow 1996; Hayasaki et al. 2007; MacFadyen & Milosavljević 2008; Roedig et al. 2011; Shi et al. 2012; Farris et al. 2014; Franchini et al. 2022; Westernacher-Schneider et al. 2022, 2024), part of the streams of gas that enter the cavity is flown back to the disc, creating shocks at the inner edge of the cavity wall that affect the disc aspect ratio. Figure 1 shows the evolution of the aspect ratio profile at t = 0, 500, 700, 900, 1000, 1300 PB for the circular cases (upper panels) and eccentric cases (lower panels) with mass ratio q = 1 (left panels) and q = 0.1 (right panels). Since Q >  1, the disc cools down until Q ∼ 1. Indeed, the aspect ratio decreases to H/R ∼ MD/MB ∼ 0.01 within the first few orbits. In all our simulations, the disc then re-expands in the vertical direction and the aspect ratio increases, eventually reaching a stable profile. This is consistent with the disc temperature increase from an initial value at the inner edge of Teff ∼ 103 K to a final value of Teff ∼ 104 K. Since gas particles are continuously flung back into the cavity wall, the shocks that occur in this region maintain the disc aspect ratio to H/R ∼ 0.08, regardless of the choice of binary mass ratio and eccentricity. The radiative cooling in the outer parts of the disc is more efficient and the aspect ratio reaches eventually a value H/R ∼ 0.04. Equal-mass binaries have a bigger quadrupole mass moment that exerts larger torques on the inner edge of the disc cavity, eventually producing more prominent streams and stronger shocks. As a consequence, discs around equal-mass binaries reach a stable configuration after ∼1000 PB (e = 0) and ∼900 PB (e = 0.9), while it takes only ∼700 PB for discs around unequal-mass binaries, regardless of eccentricity.

The evolution of the binary semi-major axis depends on the contributions to the net torque from two distinct components: the gravitational torque exerted by the circumbinary disc on each binary component and the gas accretion torque (see Roedig et al. 2011, for the detailed computation of the two contributions). Early numerical simulations of circumbinary discs showed the binary inspiral to be aided by the presence of the disc (Artymowicz & Lubow 1994, 1996; Armitage & Natarajan 2002). The findings of more recent 2D, fixed binary orbit, hydrodynamic simulations show that the angular momentum transfer onto the binary is positive and this leads the expansion of the binary (Muñoz et al. 2019, 2020; Duffell et al. 2020; Moody et al. 2019). Using a similar numerical scheme, Tiede et al. (2020) found that the fate of the binary depends on the disc temperature, i.e. on the disc aspect ratio. They find the threshold for expansion to be H/R >  0.04, while 3D smoothed particle hydrodynamics (SPH) simulations of locally isothermal discs find the threshold to be much higher, i.e. H/R = 0.2 Heath & Nixon (2020). In a recent study, Franchini et al. (2022) show, employing 3D MFM simulations where hyper-Lagrangian resolution was achieved via adaptive particle splitting (the same method employed here), that the binary inspiral-outspiral depends also on the disc viscosity. Numerical simulations that study the regime where the disc self-gravity cannot be neglected, have found that the interaction between the binary and its gaseous disc leads the binary to shrink as a consequence of the time evolution of the disc temperature, regardless of its initial value (Cuadra et al. 2009; Roedig et al. 2012; Franchini et al. 2021).

The interaction between the binary and its circumbinary disc changes the binary eccentricity as well (Goldreich & Tremaine 1980). Based on analytical arguments, Artymowicz et al. (1991) shown that an equilibrium eccentricity should exist, which was later confirmed by the SPH simulations of Roedig et al. (2011), who derived an equilibrium eccentricity 0.5 <  e <  0.8 for comparable mass binaries, linking the precise value to the disc cavity size. More recent 2D, fixed orbit, hydrodynamic simulations found equal-mass binaries to reach an equilibrium eccentricity value around ∼0.45 (Zrake et al. 2021; D’Orazio & Duffell 2021). Using a very similar numerical scheme, Siwek et al. (2023) finds that binaries with mass ratios q >  0.2 evolve towards an equilibrium eccentricity of e ∼ 0.5.

While these latter works (Franchini et al. 2022; Zrake et al. 2021; D’Orazio & Duffell 2021) assumed the disc to be locally isothermal, here we have allowed the disc temperature to change with time due to PdV work heating and radiative cooling. We find the interaction between the binary and the circumbinary disc to cause the binary to shrink, regardless of the initial conditions. In particular, a decreases by ∼1% over ∼1000 orbits in all cases, except for the binary with q = 0.1 and e = 0.9. In this case, the binary initially experiences an expansion of 0.05% over the first 600 orbits, likely because of the very high eccentricity that brings the lower mass MBH very close to the initial cavity edge. This is only a transient phase, as the binary then carves a larger cavity compared to the circular case and starts to shrink, transferring angular momentum to the circumbinary disc. We find the eccentricity value to remain relatively constant during the whole simulation for initially eccentric binaries, while circular ones feature a mild eccentricity increase, reaching e ∼ 0.04 in the equal-mass case and e ∼ 0.06 in the unequal-mass case.

3.2. Spectral energy distributions

We computed the EM emission from our numerical simulations by using Planck’s Law (see Eq. (5)), taking into account both the gas and the radiation pressure contribution when computing the disc temperature, as is mentioned in Sect. 2.2. We divided our simulated domain into five different regions: the mini-disc region that extends from the sink radius of each binary component out to the Roche Lobe size, the stream region that extends from outside the Roche Lobe out to r = 3a, and the inner (3a <  r <  5a) and outer (5a <  r <  10a) parts of the disc. For each region, the total SED was obtained by integrating the flux emitted by each pixel over the whole spatial domain.

Figure 2 shows the surface density maps, the effective temperature maps, and the SED obtained for each region at t = 1298 PB for the circular equal-mass binary. Each panel line shows a different orbital phase of the binary. Periodically, a small fraction of gas enters the cavity and feeds the mini-discs around each binary component through the streams. We assume an initial disc temperature profile that decreases with radius as R−0.5. The inflow of gas into the cavity combined with the exchange of material among the mini-discs generates shocks that increase the gas temperature, resulting in an effective temperature of Teff ∼ 104 K, warmer than the outer parts of disc.

thumbnail Fig. 1.

Time evolution of the disc aspect ratio, H/R, as a function of radius for circular binaries (top panel) and eccentric binaries (bottom panel) with a mass ratio of q = 1 (left column) and q = 0.1 (right column). For all the simulations, we report the thickness profile at different times with different colours. All the simulations start with H/R = 0.1. The aspect ratio decreases during a transient phase and increases again, reaching a constant value, similar for all the simulations: in the inner part of the disc, H/R ∼ 0.08 − 0.09, while in the outer regions of the disc, H/R ∼ 0.04 − 0.06. Unequal-mass binaries stabilise their disc temperature faster compared to equal-mass binaries.

thumbnail Fig. 2.

Surface density map, effective temperature map and SEDs for a circular equal-mass binary. From the top to the bottom row: the binary at four different orbital phases at time t = 1298 Pb. The left, middle, and right panels in each row show the surface density map in the x − y plane, the effective temperature map in the x − y plane and the SEDs. The tick spacing on the x and y axes in the left and central panels is 2a0, where a0 is the initial binary semi-major axis. The surface density upper limit is set at 1.8 × 105 g cm−2. In the SEDs, the contribution of the mini-discs around the primary and the secondary black hole is shown by the dash-dotted blue and orange lines, respectively, the stream region is represented by the dashed green line while the inner and outer part of the circumbinary disc are represented by the dash-dotted red and purple lines, respectively. The corona contribution is shown by the dotted brown line. The solid grey line shows the total emitted spectrum.

As is shown by the four different orbital phases in Fig. 2, the mini-disc and stream temperature variations occur within an orbital period, producing EM emission variations. The spectra obtained by analysing the emission from the mini-discs and the streams region (blue, orange, and green lines) do indeed exhibit more variability than the inner and outer regions of the disc (red and purple lines). The emission peak of the mini-discs and the streams region changes frequency between the optical and UV band (log (ν/Hz)∼14.6 − 15.2) spanning one order of magnitude in luminosity, while the emission peak of the circumbinary disc remains between the IR and optical band during one orbital period.

In the highly eccentric e = 0.9 case (see Fig. A.1), the emission from the cavity shows more variability than in the circular case, due to the geometry of the orbit of the components. Indeed, the emission peak of the mini-discs and the stream component changes between the optical and UV band (log (ν/Hz)∼14.8 − 15.4) spanning two orders of magnitude in luminosity (log (νLν/erg s−1) ∼ 39 − 41).

In Appendix A, we show the surface density maps, the effective temperature maps, and the SEDs for unequal-mass binaries (q = 0.1) for both the circular (Fig. A.2) and the eccentric case (Fig. A.3). We find the mini-discs to have a lower temperature in the circular unequal-mass case with respect to the circular equal-mass case shown in Fig. 2. Indeed, the exchange of material between the mini-discs does not produce significant shocks and the temperature does not increase as much as in the circular case. Therefore, the emission peak of the mini-discs is shifted to lower frequencies, log (ν/Hz)∼14.8 − 15.2, and spans ∼2 orders of magnitude in flux within one orbital period.

In the eccentric unequal-mass case (Fig. A.3), the mini-discs are significantly depleted as the high eccentricity of the binary causes them to strongly interact with each other and the material is either promptly accreted or flung back to the cavity wall. We note that in this case, the variability of the emission from the mini-discs is not totally hidden by the emission from the streams region but does change the spectrum at high frequencies within one orbital period.

In all of the simulations, by construction, the emission of the X-ray corona tracks that of the mini-discs, as per Eq. (6); therefore, depending on the simulation, it can vary by up to two orders of magnitude in luminosity.

3.3. Light curves

We integrated the luminosity emitted in each disc region in order to produce the LCs, showing the flux variation in the system as a function of time. We assume the source to be at redshift z = 0.1 and investigate the effect of redshift on source observability in the next subsection. We have calculated the bolometric flux across the whole frequency range we used to produce the spectra – 108 − 2.8 × 1019 Hz (4.13 × 10−10 − 41.3 keV) – and in three different regions of the EM spectrum: in the optical, using the VRO filters frequency bands (see Table 1), in the near-UV (1.0 − 1.5 × 1015 Hz or 4.13 − 6.20 eV), and in the soft X-ray (7.25 − 48.3 × 1016 Hz or 0.3 − 2 keV). The upcoming VRO will perform a 10-year survey of the sky in the southern hemisphere and it will potentially be able to capture the EM emission from the accretion disc of MBHBs. We therefore focus most of our attention on the optical emission of our simulated systems, and we consider the VRO Z, I, R, and G filters (see Table 1 for details), and the instrument sensitivity, as is described in Sect. 2.2.

The results are shown in Figs. 3 (e = 0) and 4 (e = 0.9). In each panel of the two figures, the left column shows the LCs computed in all the considered bands together with the accretion rate, (grey line in the top panel), while the right column shows the FFT of the accretion rate and of the fluxes computed over 300 orbits, t = 1000 − 1300 PB (t = 760 − 1060 PB for the case q = 0.1, e = 0), normalised to unity. In the VRO filters, the flux was computed including extra Gaussian stochastic fluctuations mimicking the effect of the VRO sensitivity, as is described in Sect. 2.2. In all cases, we see that the emission is brighter in the optical and UV band (except in the eccentric unequal case), while it is dimmer in the soft X-ray, consistent with our SEDs (see Figs. 2–A.3).

thumbnail Fig. 3.

Light curves and their FFT for circular mass binaries at z = 0.1. The left and right panels are for binaries with q = 1 and q = 0.1, respectively. In each panel, the first row shows the accretion rate (grey line) and the flux (orange line) integrated over the whole frequency range we consider (108 − 2.8 × 1019 Hz). The left column shows the flux evolution in time while the right column shows the FFT of the accretion rate and flux over 300 orbits in the time window t = 1000 − 1300 PB (left panel) or t = 760 − 1060 PB (right panel), normalised to unity, in function of f/fK with fK the Keplerian frequency of the binary. The second, third, and fourth rows show the flux and FFT in the optical Z band, UV band, and soft X-ray band, respectively. The optical flux was computed taking into account an extra Gaussian noise component, as is described in Sect. 2.2.

thumbnail Fig. 4.

Same as Fig. 3 but for eccentric e = 0 binaries. Here, the Fourier transform of the accretion rate and fluxes is computed over 400 orbits in the time window t = 900 − 1300 PB.

Indeed, in our model the only contribution that produces luminosity at high frequencies is the corona. Another clear trend shown in all simulations is that variability tends to increase with frequency. In fact, while the flux in the VRO bands changes within a factor of 2−3, the UV and X-ray fluxes can experience oscillations of more than one order of magnitude, in particular in eccentric cases. This is consistent with the physics of the emission from the disc. The optical mostly comes from the circumbinary disc, which is relatively steady and only mildly affected by the binary. Conversely the UV emission is dominated by the streams and mini-discs, which are strongly impacted by the dynamics of the binary, and thus highly variable. Finally, the X-ray corona is directly connected to the UV emission of the mini-discs, which is the component showing the highest variability. The large amplitude variation in the UV flux that we find suggests that these systems might be identifiable by future wide-field UV facilities, such as ULTRASAT (scheduled for launch in 2027, Shvartzvald et al. 2024) and UVEX (scheduled for launch in 2030, Kulkarni et al. 2023).

The left panel in Fig. 3 displays a number of interesting features. It is clear that the main modulation pattern is not related to the binary period, but occurs on longer timescales, and this is true both for the LCs and for the accretion rate. This is confirmed by the FFTs, which show two clear peaks at 0.2 fK and 0.4 fK (second harmonic), where fK is the Keplerian frequency of the binary. This periodicity is associated with the ‘lump’, an over-density region that obits at the cavity edge with an orbital period a few times the binary orbital period, as reported in previous works (MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Krolik 2010; Roedig et al. 2011; Noble et al. 2012; Shi et al. 2012; D’Orazio et al. 2013; Farris et al. 2014; Bowen et al. 2018; Tang et al. 2018; Westernacher-Schneider et al. 2022, 2024). It is also interesting to notice that, while the accretion rate shows an intricate structure of harmonics (Farris et al. 2014; Muñoz et al. 2019; Franchini et al. 2023), this is much less evident in the LCs, where there is significant power only at 2 fK, which corresponds to one half of the orbital period. The situation is strikingly different when the mass ratio of the binary is small (right panel in Fig. 3). In this case, both the accretion rate and the LCs show a clear periodicity on the binary period, which is confirmed by the FFT, where clear peaks are visible at 1 fK and 2 fK (second harmonic). We also note that there is no clear power at f <  fK. This is because no significant lump forms in this case, since perturbations induced by the binary are not sufficient to excite an m = 1 mode at the inner edge of the circumbinary disc cavity. This is in line with results of 2D simulations in the literature (see e.g. Farris et al. 2014).

Results for the eccentric binary simulations are shown in Fig. 4. The equal-mass case (left panel) shows an interesting periodicity structure. The fluxes emitted in the optical bands all exhibit a clear modulation of a factor ≈2, combining periodicities related to both the binary and the lump dynamics. In the FFT, we can clearly see the lump periodicity at f ≈ 0.15fK: this is lower than the circular binary case, as the cavity is larger and the period associated with its inner edge is ∝R3/2. Contrary to the circular equal-mass case, clear peaks are visible also at 1 fK and 2 fK (second harmonic), driven by the binary orbital period. In the eccentric unequal-mass case instead (right panel), the lump periodicity is absent and the peak at f ∼ 2 fK is less prominent than what was found in the circular unequal binary case. We have also explored the effect of placing the binary at different redshifts. For simplicity, we discuss the main results without including plots. Besides the obvious difference in flux, making closer binaries easier to detect, there is also some minor change difference in the displayed periodicities due to the redshifting of the spectra at z <  0.6. At redshifts of z >  0.6, fluxes are very noisy and periodicity peaks are not always distinguishable in all the optical bands. However, all the main features described here for binaries at z = 0.1 remain valid.

3.4. Periodicity identification in the VRO survey

We have so far computed the FFT over a large number of orbits in order to distinguish the different periodicities associated with binaries with different mass ratio and eccentricity. However, VRO might only have a handful of binary orbits at its disposal in its 10-year observational campaign. This is because binaries with periods shorter than a few years are expected to be primarily driven by GW emission, meaning that wide binaries (with long periods) live longer and are more likely to be present in the data (see Kelley et al. 2019). In fact, the binary period of 1 year considered here was chosen to be representative of the typical system that might be detected with VRO.

Since we output ten snapshots per binary orbit, the cadence of our simulated lightcurves is roughly 36 days. This sampling rate is sparser compared to what future VRO surveys are expected to achieve. In fact, the VRO Deep Drilling Fields will be covered with daily cadence, whereas in the most optimal scenario, the main Wide Fast Deep survey will get back on the same target at a roughly weekly cadence. Therefore, our simulated data miss high frequency information that, if present, might be detected in real observations. However, we stress here that our lightcurve sampling rate is good enough to cover periodicities occurring at the binary and at the inner disc cavity orbital frequencies, which are the most prominent signatures proposed in the literature. To get a qualitative idea of whether the periodicities found in the previous section can be identified in VRO data, we have computed the FFT over 5, 10 and 50 binary orbits along a total of 300 orbits for the circular binary cases and 400 orbits for the eccentric binaries, essentially shifting the FFT window within these few hundred orbits. We have then computed the average FFT and its standard deviation (STD) at different redshifts, z = 0.1,  0.4,  0.7.

In Fig. 5 we show the results for each of the four binaries we considered. In each plot, the first row shows the mean FFT both of the flux in the optical G band (blue line) and of the accretion rate (brown line) using a 5-orbit window at redshift z = 0.1,  0.4,  0.7 from the left to the right. The second and the third rows show the results obtained with 10- and 50-orbit windows, respectively. We note that, given a VRO survey of 10 years, 5, 10, and 50 orbits correspond to binaries with a period of two years, one year, and 2.4 months, respectively.

thumbnail Fig. 5.

Fast Fourier transform of LCs from simulations. The top panels are for simulations of circular binaries, while the bottom panels are for simulations of eccentric e = 0.9 binaries. The left and right columns are for binaries with q = 1 and q = 0.1, respectively. The first row of each case shows the FFT of the Optical G band flux (blue line) and the FFT of the accretion rate (brown line) computed over 5-orbit windows within a total of 300/400 orbital periods at redshift z = 0.1,  0.4,  0.7. The second and the third rows show the FFT computed over 10- and 50-orbit windows, respectively.

The main result of this exercise is that periodicities are more easily observed in unequal-mass binaries than in equal-mass ones. In fact, starting from the top left panel of Fig. 5, we see that periodic features in circular equal-mass binaries will be extremely hard to pick. With only 5 orbits, the lump periodicity falls in the lowest frequency bin of the Fourier decomposition and cannot be securely identified. Other periodicities are much weaker and do not show significant power in the FFT over 5 orbits. Things improve with increasing numbers of orbits and the lump periodicity clearly emerges when considering 50 orbits. Conversely, the circular unequal-mass case (top right panel) shows that the peaks at f = fK and f = 2 fK are already prominent after only 5 orbits, which is very promising, and the situation naturally improves if more orbits are sampled.

A similar situation is observed for eccentric binaries, with periodicities that are more prominent in the unequal-mass case. There are noticeable differences though. In the equal-mass case (lower left panel), periodicities at f = fK and f = 2 fK start to emerge already after 5 orbits. However, the variance is large, meaning that these features are not prominent and might be hard to detect. The situation naturally improves with the number of orbits, but a large variance remains, even after 50 orbits. It is also interesting to note that the lump frequency, clearly present in the circular case, does not seem to emerge here. In the eccentric unequal-mass case (lower right panel), the orbital periodicity is clearly identified already after 5 orbits, although the second harmonic is much less prominent than in the circular case.

In all cases, the variance tends to increase with redshift. This is a natural effect due to the inclusion of the VRO sensitivity limit in the computation. In particular, the VRO telescope could encounter substantial difficulties in detecting flux periodicities emanating from binaries with mass MB ∼ 106 M⊙ at redshift z >  0.5. This is clearly illustrated by the z = 0.7 panels of Fig. 5, where the peaks generally observable at z = 0.1,  0.4 tend to be swamped in the variance.

This does not necessarily mean that VRO cannot identify LISA MBHB precursors beyond z ≈ 0.5. In fact, the bolometric luminosity of the systems simulated here is Lbol ≈ 1042 erg s−1 ≈ 0.01 LEdd. The largely sub-Eddington luminosity is mainly due to the relatively low temperature of the gas, reaching only T ≈ 3 × 104 K, which is a factor of a few cooler than the temperature of a standard thin disc around a 106 M⊙ MBH (Frank et al. 2002). This is also expected, due to the fact that Rsink = 0.05a ∼ 100 RISCO, as is mentioned in Sect. 2.1, and the overall luminosity of the system might indeed be higher. We therefore explore here also the detectability of a putative brighter binary by simply multiplying the emission by a factor of 100 – that is, L = 100  ×  Ltrue – thus preserving all the variability properties found in our simulations. Although this is not self-consistent with our simulation, we consider it a useful exercise to assess VRO performances against lightcurves that could be representative of an M = 106 M⊙ MBHBs emitting at the Eddington luminosity, or potentially of a more massive system of M = 108 M⊙ but emitting at ≈0.01 LEdd. We repeated the process described above: we analyse the FFT of the flux over 5, 10 and 50 binary orbits along a total of 300 orbits for the circular binary cases and 400 orbits for the eccentric binaries and we compute the average FFT and its STD for different redshifts z = 1,  2,  3. Results are shown in Fig. 6. Most of the features discussed for the case L = Ltrue are still observed, but now peaks in the FFT can be easily identified up to z = 2, around cosmic noon. A full statistical assessment of VRO capabilities of correctly identifying these peaks will be the subject of future work.

thumbnail Fig. 6.

Same as Fig. 5 but assuming that the emitted luminosity is 100 times higher than the luminosity we obtained from simulations.

4. Conclusions

In this work, we computed SEDs and multi-wavelength LCs from 3D numerical simulations with hyper-Lagrangian refinement of milli-parsec-scale MBHBs embedded in thin circumbinary gaseous discs. The discs in our simulations are described by an adiabatic equation of state. We therefore allowed the gas to heat due to shocks, viscosity, and PdV work, and to cool via black-body radiation. We explored binary eccentricities of e = 0,  0.9 and a mass ratio of q = 1,  0.1, and computed LCs in different bands, placing the sources at different redshifts, z = 0.1,  0.4,  0.7.

We investigated the evolution of the disc aspect ratio, H/R, the binary semi-major axis, a, and the binary eccentricity, e. We found that, after an initial phase in which the black-body cooling dominates the gas thermodynamic evolution, the disc thickens again, reaching H/R ∼ 0.08 in the inner parts and maintaining a lower H/R ∼ 0.04 (corresponding to a lower temperature) in the outer part of the disc, regardless of the initial choice of binary mass ratio and eccentricity. The final equilibrium state is, in fact, mostly driven by the initial disc mass and radial extension, which are the same in all simulations. Therefore, self-regulation (Lodato 2007) drives all discs to reach a similar aspect ratio at the end of the transient phase. We find the interaction between the binary and the circumbinary disc to cause the binary to shrink, regardless of the initial conditions. Since in our model the temperature changes with time, this result further supports previous findings in the literature (Roedig & Sesana 2014; Franchini et al. 2021). We find that circular orbits tend towards higher eccentricity values, in agreement with previous works (e.g. D’Orazio & Duffell 2021; Siwek et al. 2023), whereas very eccentric binaries experience a negligible eccentricity evolution within the timeframe of our simulations. We notice, however, that we followed the evolution of the relaxed disc only for about 400 orbits, corresponding to 400 years. It might be that the eccentricity evolution for very eccentric binaries occurs on a longer timescale than the one simulated here. Indeed, here we simulated the binary evolution over ∼1300 yr, while in both the cited works the timescale is longer.

We computed the SEDs from the circumbinary discs in our simulations assuming black-body emission. We find that the luminosity emitted by the innermost region of the disc – the mini-discs and the streams – exhibits more variability than the outer parts of the disc. In the circular equal-mass case, the emission peak of the mini-discs and the streams region changes frequency between the optical and UV band (log (ν/Hz)∼14.8 − 15.4), spanning one order of magnitude in luminosity. In the unequal-mass case, the mini-discs have a lower temperature due to the absence of shocks produced by the exchange of material among the binary components. Thus, the emission peak of the mini-discs is shifted to slightly lower frequencies – log (ν/Hz)∼14.8 − 15.2 – and spans up to ∼2 orders of magnitude in luminosity. In both the eccentric binary cases, we find that the emission from the mini-discs is completely (in the equal-mass case) or partially (in the unequal-mass case) covered by the streams emission. The X-ray photons are provided by the corona, which we assume to have an emission proportional to that of the mini-discs. Therefore, the luminosity in the X-ray band displays the highest variability, which can reach up to two orders of magnitude for unequal-mass binaries.

We computed the LCs in different frequency bands, mainly focussing on the optical window that will be probed by VRO. We calculated the (thermal) flux emitted over the whole EM spectrum that we used to produce the SEDs; that is, within the frequency band 108 − 2.8 × 1019 Hz (4.13 × 10−10 − 41.3 keV), in the optical frequency band using the VRO filters, in the near-UV band within 1.0 − 1.5 × 1015 Hz (4.13 − 6.20 eV), and in the soft X-ray band, in the range of 7.25 − 48.3 × 1016 Hz (0.3 − 2 keV).

In almost all cases, the flux is notably higher in the UV band, while in the soft X-ray is dimmer, consistent with the shape of the SEDs. As the frequency increases, the flux variability grows, in particular in the UV and soft X-ray bands. Here, the flux oscillates by more than an order of magnitude, while in the VRO optical bands fluxes vary within a factor of two, in line with the physics of emission from the disc.

In the circular equal-mass case, both the flux and the accretion rate FFTs reveal clear peaks at 0.2 fK and 0.4 fK associated with the lump periodicity. Moreover, in the LCs, a peak at 2 fK is also present, which corresponds to a periodicity of one half the binary orbital period. In the eccentric equal-mass case, the lump periodicity is significant only in the optical and UV fluxes, while its amplitude is negligible in the soft X-ray band. This is probably due to the larger cavity carved by the binary, which causes the lump region to emit in the optical/UV band rather than in the soft X-ray band. We indeed found the lump modulation peak to be shifted to fK ∼ 0.15. This is consistent with the fact that the cavity is larger than in the circular equal-mass case since the period associated with its inner edge is ∝R3/2. We note that we found evidence of lump periodicity only in equal-mass binaries. This is consistent with previous works that show the lump modulation amplitude to decrease with binary mass ratio (D’Orazio et al. 2013). Therefore, a lack of sub-orbital modulation in the presence of a clear orbital modulation might indicate a small binary mass ratio. We also found a prominent flux and accretion rate modulation over the orbital period of the binary and half of it in all the simulations, with the exception of the circular equal-mass case, which shows a a weak periodicity only at 2 fK. By exploring different redshifts, there are some minor change in periodicities at a redshift of z <  0.6, while at higher redshifts fluxes are very noisy and the periodicities are not always well distinguished.

All the aforementioned considerations are valid when the FFT is computed over a large number of binary orbits (i.e. 300 − 400 PB). However, the VRO survey is planned for 10 years and most compact MBHBs are expected to have periods of ≈years. Therefore, we have computed the FFT of the flux and of the accretion rate over 5, 10, and 50 binary obits at different redshifts to make a more realistic assessment of the prominence of these periodicities in the VRO data. In the circular equal-mass case, detecting periodicities is challenging, in particular with only five orbits. We found that, as the number of orbits used to compute the FFT increases, the periodicities become more distinct and the associated variance decreases, as was expected. Still, the identification of equal-mass, circular binaries appears to be the most challenging. Conversely, binaries with q = 0.1 show promising results, with distinguishable peaks at 1 fK and 2 fK even after five orbits. Similar trends are observed in eccentric cases. The lump periodicity is totally absent in all cases except for the circular equal-mass case, which shows a hint of lump periodicity when computing the FFT over 50 orbits. Thus, the chances of detecting it during a 10-year survey by assuming unequal MBHB with an orbital period of 1 year are negligible.

Due to the intrinsic faintness of our system, Lbol ≈ 1042 erg s−1, the detectability of periodicities with VRO are limited to systems at z <  0.5. As an exercise, we increased the luminosity of all our systems by a factor of 100, mimicking a MBHB of M = 106 M⊙ emitting at the Eddington limit. We found that in this case, periodicities can be identified by VRO up to z ≈ 2, opening the possibility of finding these systems in a large fraction of the co-moving volume of the Universe.

In general, our results indicate that periodicities related to unequal-mass binaries will be easier to identify in VRO data, compared to equal-mass ones. Moreover, LCs from binaries with different properties (equal vs. unequal mass, circular vs. eccentric) are characterised by different distinctive modulations, hinting at the possibility (to be further investigated) of constraining the binary properties from their time-domain EM data. These data can inform current and future low-frequency GW searches with LISA and PTAs at different levels. Working at millihertz frequencies, LISA will detect MBHBs with orbital periods of ∼hours, in their late inspiral and final coalescence. Therefore, VRO can provide a census of LISA precursors; binaries that will coalesce in the LISA band within the next few million years. Measurements of their mass ratios, eccentricities, and environment from EM observations can then be used to set priors for LISA searches. Moreover, combined GW and EM information can be used to test our physical understanding of MBHB evolution. For example, we expect binaries evolving in circumbinary gaseous discs to feature high, aligned spins (Dotti et al. 2010). Therefore, a hypothetical large population of gas-rich MBHBs detected by VRO should translate into a number of spin-aligned systems detected by LISA. Conversely, working at nanohertz frequencies, PTAs have the potential to identify MBHBs with orbital periods comparable to those probed by VRO. It should be noted, however, that systems likely to be individually resolved by PTAs have masses ≳109 M⊙, orders of magnitude larger than those simulated here (Rosado et al. 2015). If we assume that the trends we found with mass ratio and eccentricities in the EM signatures hold also at much higher masses, then putative multi-messenger observations of PTA sources will allow us to constrain the parameters of the system. In fact, most PTA resolved sources will be monochromatic, meaning that from the GW signal we shall only be able to extract a wave amplitude. The identification of an EM counterpart will first provide a measurement of the distance to the source, allowing one to estimate the source chirp mass from the amplitude. Moreover, distinctive periodicity features might help in constraining the binary mass ratio, allowing one to estimate the individual masses of the two black holes. It is therefore clear that refining our understanding of EM signals from MBHBs might prove extremely fruitful in the era of multi-messenger astronomy.

In this work, we present advanced 3D simulations that improve the description of accreting MBHBs embedded in gaseous discs, which is fundamental in order to make more realistic predictions of the emission signatures of these sources. We note here that we included the effect that radiation pressure has in determining the gas temperature only a posteriori. The inclusion of radiation pressure on the fly is the subject of a future work (Cocchiararo et al., in prep.). The inclusion of radiation pressure will allow us to obtain a more comprehensive model of these systems and possibly a better characterisation of their EM signatures. Finally, a wider exploration of the binary-disc parameter space using our numerical simulations is needed to make more robust observational predictions and will be the subject of future work.

Acknowledgments

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. AF and AS acknowledge financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691). AL acknowledges support by the PRIN MUR “2022935STW”. FC thanks Nataliya Porayko, Golam Shaifullah and Enrico Panontin for helpful discussions and suggestions.

References

  1. Afzal, A., Agazie, G., Anumarlapudi, A., et al. 2023, ApJ, 951, L11 [NASA ADS] [CrossRef] [Google Scholar]
  2. Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023a, ApJ, 951, L8 [NASA ADS] [CrossRef] [Google Scholar]
  3. Agazie, G., Alam, M. F., Anumarlapudi, A., et al. 2023b, ApJ, 951, L9 [CrossRef] [Google Scholar]
  4. Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023c, ApJ, 951, L10 [CrossRef] [Google Scholar]
  5. Amaro-Seoane, P., Andrews, J., Arca Sedda, M., et al. 2023, Liv. Rev. Relat., 26, 2 [NASA ADS] [CrossRef] [Google Scholar]
  6. Armitage, P. J., & Natarajan, P. 2002, ApJ, 567, L9 [NASA ADS] [CrossRef] [Google Scholar]
  7. Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [Google Scholar]
  8. Artymowicz, P., & Lubow, S. H. 1996, ApJ, 467, L77 [NASA ADS] [CrossRef] [Google Scholar]
  9. Artymowicz, P., Clarke, C. J., Lubow, S. H., & Pringle, J. E. 1991, ApJ, 370, L35 [NASA ADS] [CrossRef] [Google Scholar]
  10. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
  11. Barrows, R. S., Stern, D., Madsen, K., et al. 2011, ApJ, 744, 7 [Google Scholar]
  12. Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362 [Google Scholar]
  13. Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [Google Scholar]
  14. Bogdanović, T., Eracleous, M., & Sigurdsson, S. 2009, New Astron. Rev., 53, 113 [CrossRef] [Google Scholar]
  15. Bogdanović, T., Miller, M. C., & Blecha, L. 2022, Liv. Rev. Relat., 25, 3 [Google Scholar]
  16. Boroson, T. A., & Lauer, T. R. 2009, Nature, 458, 53 [Google Scholar]
  17. Bortolas, E., Franchini, A., Bonetti, M., & Sesana, A. 2021, ApJ, 918, L15 [NASA ADS] [CrossRef] [Google Scholar]
  18. Bowen, D. B., Mewes, V., Campanelli, M., et al. 2018, ApJ, 853, L17 [NASA ADS] [CrossRef] [Google Scholar]
  19. Chandrasekhar, S. 1943, ApJ, 97, 255 [Google Scholar]
  20. Charisi, M., Bartos, I., Haiman, Z., et al. 2016, MNRAS, 463, 2145 [Google Scholar]
  21. Chen, Y.-C., Liu, X., Liao, W.-T., et al. 2020, MNRAS, 499, 2245 [Google Scholar]
  22. Combi, L., Lopez Armengol, F. G., Campanelli, M., et al. 2022, ApJ, 928, 187 [NASA ADS] [CrossRef] [Google Scholar]
  23. Cuadra, J., Armitage, P. J., Alexander, R. D., & Begelman, M. C. 2009, MNRAS, 393, 1423 [Google Scholar]
  24. d’Ascoli, S., Noble, S. C., Bowen, D. B., et al. 2018, ApJ, 865, 140 [CrossRef] [Google Scholar]
  25. Decarli, R., Dotti, M., Montuori, C., Liimets, T., & Ederoclite, A. 2010, ApJ, 720, L93 [NASA ADS] [CrossRef] [Google Scholar]
  26. D’Orazio, D. J., & Duffell, P. C. 2021, ApJ, 914, L21 [CrossRef] [Google Scholar]
  27. D’Orazio, D. J., Haiman, Z., & MacFadyen, A. 2013, MNRAS, 436, 2997 [Google Scholar]
  28. Dotti, M., Montuori, C., Decarli, R., et al. 2009, MNRAS, 398, L73 [NASA ADS] [Google Scholar]
  29. Dotti, M., Volonteri, M., Perego, A., et al. 2010, MNRAS, 402, 682 [NASA ADS] [CrossRef] [Google Scholar]
  30. Dotti, M., Sesana, A., & Decarli, R. 2012, Adv. Astron., 2012 [CrossRef] [Google Scholar]
  31. Duffell, P. C., D’Orazio, D., Derdzinski, A., et al. 2020, ApJ, 901, 25 [NASA ADS] [CrossRef] [Google Scholar]
  32. Duras, F., Bongiorno, A., Ricci, F., et al. 2020, A&A, 636, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. EPTA Collaboration (Antoniadis, J., et al.) 2023, A&A, 678, A48 [Google Scholar]
  34. EPTA Collaboration& InPTA Collaboration (Antoniadis, J., et al.) 2023a, A&A, 678, A50 [CrossRef] [EDP Sciences] [Google Scholar]
  35. EPTA Collaboration& InPTA Collaboration (Antoniadis, J., et al.) 2023b, A&A, 678, A49 [Google Scholar]
  36. EPTA Collaboration& InPTA Collaboration (Antoniadis, J., et al.) 2024a, A&A, 690, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. EPTA Collaboration& InPTA Collaboration (Antoniadis, J., et al.) 2024b, A&A, 685, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Eracleous, M., & Halpern, J. P. 1994, ApJS, 90, 1 [NASA ADS] [CrossRef] [Google Scholar]
  39. Eracleous, M., Boroson, T. A., Halpern, J. P., & Liu, J. 2012, ApJS, 201, 23 [Google Scholar]
  40. Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2005, ApJ, 630, 152 [Google Scholar]
  41. Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z. 2014, ApJ, 783, 134 [NASA ADS] [CrossRef] [Google Scholar]
  42. Foster, R. S., & Backer, D. C. 1990, ApJ, 361, 300 [NASA ADS] [CrossRef] [Google Scholar]
  43. Franchini, A., Sesana, A., & Dotti, M. 2021, MNRAS, 507, 1458 [NASA ADS] [Google Scholar]
  44. Franchini, A., Lupi, A., & Sesana, A. 2022, ApJ, 929, L13 [NASA ADS] [CrossRef] [Google Scholar]
  45. Franchini, A., Lupi, A., Sesana, A., & Haiman, Z. 2023, MNRAS, 522, 1569 [NASA ADS] [CrossRef] [Google Scholar]
  46. Franchini, A., Bonetti, M., Lupi, A., & Sesana, A. 2024, A&A, 686, A288 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition (Cambridge: Cambridge University Press) [Google Scholar]
  48. Gammie, C. F. 2001, ApJ, 553, 174 [Google Scholar]
  49. Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [Google Scholar]
  50. Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, Nature, 518, 74 [Google Scholar]
  51. Gutiérrez, E. M., Combi, L., Noble, S. C., et al. 2022, ApJ, 928, 137 [CrossRef] [Google Scholar]
  52. Hayasaki, K., Mineshige, S., & Sudou, H. 2007, PASJ, 59, 427 [NASA ADS] [Google Scholar]
  53. Heath, R. M., & Nixon, C. J. 2020, A&A, 641, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Hogg, D. W. 1999, ArXiv e-prints [arXiv:astro-ph/9905116] [Google Scholar]
  55. Hogg, D. W., Baldry, I. K., Blanton, M. R., & Eisenstein, D. J. 2002, ArXiv e-prints [arXiv:astro-ph/0210394] [Google Scholar]
  56. Hopkins, P. F. 2015, MNRAS, 450, 53 [Google Scholar]
  57. Kelley, L. Z., Haiman, Z., Sesana, A., & Hernquist, L. 2019, MNRAS, 485, 1579 [NASA ADS] [CrossRef] [Google Scholar]
  58. Khan, F. M., Just, A., & Merritt, D. 2011, ApJ, 732, 89 [Google Scholar]
  59. Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [Google Scholar]
  60. Krauth, L. M., Davelaar, J., Haiman, Z., et al. 2023, MNRAS, 526, 5441 [NASA ADS] [CrossRef] [Google Scholar]
  61. Krolik, J. H. 2010, ApJ, 709, 774 [NASA ADS] [CrossRef] [Google Scholar]
  62. Kulkarni, S. R., Harrison, F. A., Grefenstette, B. W., et al. 2023, ArXiv e-prints [arXiv:2111.15608] [Google Scholar]
  63. Liu, T., Gezari, S., Ayers, M., et al. 2019, ApJ, 884, 36 [Google Scholar]
  64. Lodato, G. 2007, Nuovo Cimento Rivista Serie, 30, 293 [NASA ADS] [Google Scholar]
  65. Lodato, G., & Rice, W. K. M. 2004, MNRAS, 351, 630 [Google Scholar]
  66. Lodato, G., & Rice, W. K. M. 2005, MNRAS, 358, 1489 [Google Scholar]
  67. Lodato, G., Nayakshin, S., King, A. R., & Pringle, J. E. 2009, MNRAS, 398, 1392 [NASA ADS] [CrossRef] [Google Scholar]
  68. LSST Science Collaboration (Abell, P. A., et al.) 2009, ArXiv e-prints [arXiv:0912.0201] [Google Scholar]
  69. MacFadyen, A. I., & Milosavljević, M. 2008, ApJ, 672, 83 [Google Scholar]
  70. Mayer, L., Kazantzidis, S., Madau, P., et al. 2007, Science, 316, 1874 [NASA ADS] [CrossRef] [Google Scholar]
  71. Moody, M. S. L., Shi, J.-M., & Stone, J. M. 2019, ApJ, 875, 66 [NASA ADS] [CrossRef] [Google Scholar]
  72. Muñoz, D. J., Miranda, R., & Lai, D. 2019, ApJ, 871, 84 [CrossRef] [Google Scholar]
  73. Muñoz, D. J., Lai, D., Kratter, K., & Miranda, R. 2020, ApJ, 889, 114 [Google Scholar]
  74. Noble, S. C., Mundim, B. C., Nakano, H., et al. 2012, ApJ, 755, 51 [NASA ADS] [CrossRef] [Google Scholar]
  75. Peñil, P., Ajello, M., Buson, S., et al. 2022, ArXiv e-prints [arXiv:2211.01894] [Google Scholar]
  76. Preto, M., Berentzen, I., Berczik, P., & Spurzem, R. 2011, ApJ, 732, L26 [Google Scholar]
  77. Price, D. J. 2007, PASA, 24, 159 [Google Scholar]
  78. Price, D. J., Wurster, J., Nixon, C., et al. 2017, Astrophysics Source Code Library [record ascl:1709.002] [Google Scholar]
  79. Quinlan, G. D. 1996, New Astron., 1, 35 [Google Scholar]
  80. Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6 [NASA ADS] [CrossRef] [Google Scholar]
  81. Regan, J. A., Downes, T. P., Volonteri, M., et al. 2019, MNRAS, 486, 3892 [CrossRef] [Google Scholar]
  82. Roedig, C., & Sesana, A. 2014, MNRAS, 439, 3476 [Google Scholar]
  83. Roedig, C., Dotti, M., Sesana, A., Cuadra, J., & Colpi, M. 2011, MNRAS, 415, 3033 [NASA ADS] [CrossRef] [Google Scholar]
  84. Roedig, C., Sesana, A., Dotti, M., et al. 2012, A&A, 545, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Rosado, P. A., Sesana, A., & Gair, J. 2015, MNRAS, 451, 2417 [NASA ADS] [CrossRef] [Google Scholar]
  86. Sandrinelli, A., Covino, S., Treves, A., et al. 2018, A&A, 615, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Sesana, A., Haardt, F., & Madau, P. 2007, ApJ, 660, 546 [NASA ADS] [CrossRef] [Google Scholar]
  88. Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192 [NASA ADS] [CrossRef] [Google Scholar]
  89. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  90. Shi, J.-M., Krolik, J. H., Lubow, S. H., & Hawley, J. F. 2012, ApJ, 749, 118 [Google Scholar]
  91. Shvartzvald, Y., Waxman, E., Gal-Yam, A., et al. 2024, ApJ, 964, 74 [NASA ADS] [CrossRef] [Google Scholar]
  92. Siwek, M., Weinberger, R., & Hernquist, L. 2023, MNRAS, 522, 2707 [NASA ADS] [CrossRef] [Google Scholar]
  93. Smarra, C., Goncharov, B., Barausse, E., et al. 2023, ArXiv e-prints [arXiv:2306.16228] [Google Scholar]
  94. Tanaka, T., Perna, R., & Haiman, Z. 2012, MNRAS, 425, 2974 [CrossRef] [Google Scholar]
  95. Tang, S., & Grindlay, J. 2009, ApJ, 704, 1189 [NASA ADS] [CrossRef] [Google Scholar]
  96. Tang, Y., Haiman, Z., & MacFadyen, A. 2018, MNRAS, 476, 2249 [NASA ADS] [CrossRef] [Google Scholar]
  97. Tiede, C., Zrake, J., MacFadyen, A., & Haiman, Z. 2020, ApJ, 900, 43 [NASA ADS] [CrossRef] [Google Scholar]
  98. Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
  99. Tsai, C.-W., Jarrett, T. H., Stern, D., et al. 2013, ApJ, 779, 41 [NASA ADS] [CrossRef] [Google Scholar]
  100. Tsalmantza, P., Decarli, R., Dotti, M., & Hogg, D. W. 2011, ApJ, 738, 20 [Google Scholar]
  101. Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, MNRAS, 461, 3145 [Google Scholar]
  102. Westernacher-Schneider, J. R., Zrake, J., MacFadyen, A., & Haiman, Z. 2022, Phys. Rev. D, 106, 103010 [NASA ADS] [CrossRef] [Google Scholar]
  103. Westernacher-Schneider, J. R., Zrake, J., MacFadyen, A., & Haiman, Z. 2024, ApJ, 962, 76 [NASA ADS] [CrossRef] [Google Scholar]
  104. Xu, H., Chen, S., Guo, Y., et al. 2023, RAA, 23, 075024 [NASA ADS] [Google Scholar]
  105. Zrake, J., Tiede, C., MacFadyen, A., & Haiman, Z. 2021, ApJ, 909, L13 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: SEDs from unequal and/or eccentric binaries

thumbnail Fig. A.1.

Same as Fig. 2 but for the highly eccentric e = 0.9 equal-mass binary.

thumbnail Fig. A.2.

Same as Fig. 2 but for a binary with mass ratio q = 0.1 and eccentricity e = 0 and at time t = 1000 Pb.

thumbnail Fig. A.3.

Same as Fig. 2 but for a binary with mass ratio q = 0.1 and eccentricity e = 0.9 and at time t = 1298 Pb.

All Tables

Table 1.

Frequency bands and 5-σ flux sensitivities of the VRO telescope’s optical filters.

In the text

All Figures

thumbnail Fig. 1.

Time evolution of the disc aspect ratio, H/R, as a function of radius for circular binaries (top panel) and eccentric binaries (bottom panel) with a mass ratio of q = 1 (left column) and q = 0.1 (right column). For all the simulations, we report the thickness profile at different times with different colours. All the simulations start with H/R = 0.1. The aspect ratio decreases during a transient phase and increases again, reaching a constant value, similar for all the simulations: in the inner part of the disc, H/R ∼ 0.08 − 0.09, while in the outer regions of the disc, H/R ∼ 0.04 − 0.06. Unequal-mass binaries stabilise their disc temperature faster compared to equal-mass binaries.
In the text
thumbnail Fig. 2.

Surface density map, effective temperature map and SEDs for a circular equal-mass binary. From the top to the bottom row: the binary at four different orbital phases at time t = 1298 Pb. The left, middle, and right panels in each row show the surface density map in the x − y plane, the effective temperature map in the x − y plane and the SEDs. The tick spacing on the x and y axes in the left and central panels is 2a0, where a0 is the initial binary semi-major axis. The surface density upper limit is set at 1.8 × 105 g cm−2. In the SEDs, the contribution of the mini-discs around the primary and the secondary black hole is shown by the dash-dotted blue and orange lines, respectively, the stream region is represented by the dashed green line while the inner and outer part of the circumbinary disc are represented by the dash-dotted red and purple lines, respectively. The corona contribution is shown by the dotted brown line. The solid grey line shows the total emitted spectrum.
In the text
thumbnail Fig. 3.

Light curves and their FFT for circular mass binaries at z = 0.1. The left and right panels are for binaries with q = 1 and q = 0.1, respectively. In each panel, the first row shows the accretion rate (grey line) and the flux (orange line) integrated over the whole frequency range we consider (108 − 2.8 × 1019 Hz). The left column shows the flux evolution in time while the right column shows the FFT of the accretion rate and flux over 300 orbits in the time window t = 1000 − 1300 PB (left panel) or t = 760 − 1060 PB (right panel), normalised to unity, in function of f/fK with fK the Keplerian frequency of the binary. The second, third, and fourth rows show the flux and FFT in the optical Z band, UV band, and soft X-ray band, respectively. The optical flux was computed taking into account an extra Gaussian noise component, as is described in Sect. 2.2.
In the text
thumbnail Fig. 4.

Same as Fig. 3 but for eccentric e = 0 binaries. Here, the Fourier transform of the accretion rate and fluxes is computed over 400 orbits in the time window t = 900 − 1300 PB.
In the text
thumbnail Fig. 5.

Fast Fourier transform of LCs from simulations. The top panels are for simulations of circular binaries, while the bottom panels are for simulations of eccentric e = 0.9 binaries. The left and right columns are for binaries with q = 1 and q = 0.1, respectively. The first row of each case shows the FFT of the Optical G band flux (blue line) and the FFT of the accretion rate (brown line) computed over 5-orbit windows within a total of 300/400 orbital periods at redshift z = 0.1,  0.4,  0.7. The second and the third rows show the FFT computed over 10- and 50-orbit windows, respectively.
In the text
thumbnail Fig. 6.

Same as Fig. 5 but assuming that the emitted luminosity is 100 times higher than the luminosity we obtained from simulations.
In the text
thumbnail Fig. A.1.

Same as Fig. 2 but for the highly eccentric e = 0.9 equal-mass binary.
In the text
thumbnail Fig. A.2.

Same as Fig. 2 but for a binary with mass ratio q = 0.1 and eccentricity e = 0 and at time t = 1000 Pb.
In the text
thumbnail Fig. A.3.

Same as Fig. 2 but for a binary with mass ratio q = 0.1 and eccentricity e = 0.9 and at time t = 1298 Pb.
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.