Hydrodynamic Torques in Circumbinary Accretion Disks

TL;DR

This study uses 3D hydrodynamic simulations to analyze how circumbinary disks influence the orbital evolution of equal-mass binary systems, finding that disks tend to cause the binaries to grow in semi-major axis.

Contribution

First detailed 3D hydrodynamic simulations of circumbinary disks showing they cause binary orbital expansion regardless of disk alignment.

Findings

Total torque on the binary is positive, leading to orbital expansion.

Binary growth occurs in both aligned and misaligned disk configurations.

Disk realignment timescale matches the viscous timescale of the disk.

Abstract

Gaseous disks have been proposed as a mechanism for facilitating mergers of binary black holes. We explore circumbinary disk systems to determine the evolution of the central binary. To do so, we perform 3D, hydrodynamic, locally isothermal simulations of circumbinary disks on a Cartesian grid. We focus on binaries of equal mass ratios on fixed circular orbits. To investigate the orbital evolution of the binary, we examine the various torques exerted on the system. For the case where the disk plane and binary orbital plane are aligned, we find that the total torque is positive so that the semi-major axis of the binary increases. For the misaligned case, we run simulations with the binary orbital plane and disk midplane misaligned by 45\degree and find the same results - the binary grows. The timescale for the circumbinary disk to realign to the plane of the binary is consistent with the global viscous timescale of the disk.

Figures

Tables

Table 1. Table 1 Simulation Parameters

Run	Dim	Resolution	Inclination	Sink time ($t_{\text{rm}}$)
A	2D	$256\times 256$	-	29.814
B	2D	$256\times 256$	-	5.0
C	2D	$256\times 256$	-	0.1
D	2D	$128\times 128$	-	29.814
E	3D	$256\times 256\times 128$	0°	29.814
F	3D	$256\times 256\times 128$	45°	29.814

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0,$$

$$\frac{\partial \rho \mathbf{v}}{\partial t}+\nabla \cdot \left(\rho \mathbf{vv}+P\mathbf{I}+{\mathbf{T}}_{\mathrm{vis}}\right)=-\rho \nabla \varphi ,$$

$$c_{\mathrm{s}}^2=h^2|\varphi |,$$

$$\varphi =-\frac{G{M}_1}{|\mathbf{r}-{\mathbf{r}}_1|}-\frac{G{M}_2}{|\mathbf{r}-{\mathbf{r}}_2|}$$

$${\mathbf{T}}_{\mathrm{vis},ij}=-\rho \nu \left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}-\frac{2}{3}\frac{\partial v_k}{\partial x_k}{\delta }_{ij}\right)$$

$$\rho (R,z){|}_{t=0}=\rho (R,z=0)\exp \left[\frac{GM}{c_s^2}\left(\frac{1}{\sqrt{R^2+z^2}}-\frac{1}{R}\right)\right]$$

$$\rho (R,z=0)={\rho }_0(R/a{)}^{-p}\exp \left[-(R/R_s{)}^{-2}\right]$$

$$du/dt=-[u-u(t=0)]\text{sl}F$$

$$\text{sl}F=\frac{R^2-R_{i}R}{R_o^2-R_i{R}_o}{t}_{\mathrm{orb}}(R_i{)}^{-1}$$

$$d\rho /dt=-\rho /t_{\mathrm{rm}},$$

$$t_{\mathrm{rm}}=t_{\mathrm{visc}}(|\mathbf{r}-{\mathbf{r}}_{\mathbf{i}}|)=\frac{2|\mathbf{r}-{\mathbf{r}}_{\mathbf{i}}{|}^2}{3\nu }$$

$${\dot{\mathbf{J}}}_{adv}=\int \rho \mathbf{v}\times \mathbf{r}(\mathbf{v}\cdot \mathbf{r})dS$$

$${\dot{\mathbf{J}}}_{visc}=-\int ({\mathbf{T}}_{visc}\cdot \widehat{{\mathbf{e}}_{\mathbf{r}}})\times \mathbf{r}dS$$

$${\mathbf{T}}_{grav}=\int \rho (\mathbf{r}\times \nabla \varphi )dS$$

$${\dot{\mathbf{J}}}_{tot}={\dot{\mathbf{J}}}_{adv}-{\dot{\mathbf{J}}}_{visc}-{\mathbf{T}}_{grav}$$

$$\frac{\dot{a}}{a}=8\Big{(}\frac{l_{0}}{l_{B}}-\frac{3}{8}\Big{)}\frac{\dot{M}}{M}$$

$$\cos \delta =\frac{{\mathbf{J}}_{\mathbf{B}}\cdot \mathbf{L}}{\mid {\mathbf{J}}_{\mathbf{B}}\mid \mid \mathbf{L}\mid }$$

$$T_{\nu }=\frac{R^2}{\nu }=\frac{1}{\alpha h^2}\sqrt{\frac{R^3}{GM}}.$$

$$R_{\mathrm{break}}\lesssim 50\mu^{1/2}{|\sin{2\theta}|}^{1/2}{\Big{(}\frac{H/R}{10^{-3}}\Big{)}}^{-1/2}{\Big{(}\frac{\alpha}{0.1}\Big{)}}^{-1/2}a.$$

$$\frac{d}{dt}\Big{(}\frac{1}{a}\Big{)}_{\textrm{stars}}=-14.3\frac{G\rho}{\sigma}$$

$$\Big{(}\frac{da}{dt}\Big{)}_{\textrm{stars}}=14.3\frac{G\rho}{\sigma}a^{2}\ .$$

Full text

††slugcomment: Accepted to ApJ.

1. Introduction

Circumbinary accretion disks are a common occurrence in nature, in systems such as stellar binaries (e.g. Dutrey et al., 1994; Mathieu et al., 1997; Orosz et al., 2012) and supermassive black hole binaries (e.g. Begelman et al., 1980; Milosavljević & Phinney, 2005; Dotti et al., 2012). The interactions between the binary and the disk are of great interest for understanding the evolution of the binary, and of the system as a whole. For example, the existence of circumbinary planets indicate the importance of understanding these systems, in order to understand the formation of such planets. Many such systems have been detected by the Kepler mission (Doyle et al., 2011; Welsh et al., 2014).

Black hole binaries are of particular interest, as the tools now exist to observe black hole - black hole mergers using LIGO (Abbott et al., 2017). It is an ongoing question what the mechanism is by which supermassive black hole binaries shrink to reach semi-major axes smaller than about a parsec, known as the final parsec problem (Begelman et al., 1980). A gaseous disk might pslay an important role in the evolution of the system, and may provide a mechanism for facilitating the shrinkage of the binary to these separations (Escala et al., 2005). Though this is a popular mechanism, the gas disk is not a definitive solution. The work of del Valle & Volonteri (2018) and Souza Lima et al. (2017) highlight that AGN feedback in a gas disk causes the shrinkage rate of the binary to significantly decrease. Souza Lima et al. (2017) also explore radiative cooling and find that it causes fragmentation of the disk, which leads to additional slowing of the binary decay rate.

To answer the question of the validity of gas disks as a shrinkage mechanism, it is necessary to understand the angular momentum transport within the circumbinary disk. If there were no accretion of matter on to the binary, then the binary would only lose angular momentum to the circumbinary disk, potentially presenting a solution to the final parsec problem. However, taking accretion into consideration complicates the situation, as accretion causes the binary to gain angular momentum. Therefore, the overall sign of the angular momentum transport is uncertain, and is key in understanding the evolution of the binary.

Additionally, we expect some fraction of circumbinary disks to be misaligned with the binary orbital plane. There is observational evidence that shows this misalignment in binary star systems such as HK Tau (Stapelfeldt et al., 1998; Jensen & Akeson, 2014), KH 15D (Winn et al., 2004), and IRS 43 (Brinch et al., 2016). In binary black hole systems, Nixon et al. (2013) note that accretion events on to supermassive black hole binaries are likely to be chaotic and randomly oriented, as the separation of the binary is on a much smaller scale than the scale of the galaxy feeding the accretion events. Also, internal galactic processes can randomize angular momentum of the gas involved in these accretion events, further leading to the conclusion that we should not expect accretion to be aligned with the binary and therefore can expect some circumbinary disks to form misaligned. Dunhill et al. (2014) have simulated such accretion events onto supermassive black hole binaries, and found that this does indeed results in disks misaligned with the binary orbital plane.

There has been much work done on hydrodynamic simulations of binary systems with a circumbinary accretion disk. In particular, both Miranda et al. (2017) (hereafter MML17) and Tang et al. (2017) have done 2D simulations using a polar grid and examine the torques on the system. MML17 use the grid based code PLUTO, excise a central region and do not resolve the binary, and examine a range of mass ratios. They find that for circular binaries the binary semi-major axis will expand, though this changes with viscosity and binary eccentricity. Tang et al. (2017) perform similar simulations using the moving mesh code DISCO, though they resolve the whole simulation domain, and focus on the case of an equal mass binary. They examine a range of removal rates for the sink cells at the center and find that for faster sink times the net torque is positive, and for slower sink times the net torque is negative which drives the binary towards a merger. Their work builds on the work of Farris et al. (2014) who use the same moving mesh code and examine a range of mass ratios, but do not examine the torques on the system. In their simulations they find strong periodicity of the accretion rate which is associated with a high density lump, and for unequal mass ratios find that accretion causes the mass ratio to increase. The work of Nixon et al. (2013) uses the smooth particle hydrodynamics code PHANTOM to perform simulations of misaligned disks for equal mass binaries, focusing on examining the tearing of the disk. For more inclined disks, they find significant tearing of the entire disk and enhanced accretion rates. Fragner & Nelson (2010) use the grid based code NIRVANA to perform simulations of a binary system with a slightly different configuration - there is a central primary with a secondary that is far from the disk, and misaligned from the disk midplane. They examine the disk inclination and precession evolution closely, and find that thick disks with low viscosity align with the binary orbit on the viscous evolution timescale. Thin disks with higher viscosity also align with the binary orbit on similar timescales but develop significant twists. Most recently, the work of Muñoz et al. (2018) builds on the work done previously in MML17. They use the moving mesh code AREPO to examine equal mass binaries with circular and eccentric orbits, performing a close study of the angular momentum transfer to the binary and the resulting binary orbital evolution. They find in all of their simulations that the binary expands, similar to the results of MML17.

In this paper we present 2D and 3D viscous hydrodynamic simulations of circumbinary disk systems, using the finite-volume MHD code athena++. We adopt a simple prescription for the disk, using a thin ($h\equiv H/R=0.1$) $\alpha$ disk model with no magnetic fields or radiative transport. The system is evolved for 1000 binary orbits, to relax to a quasi-steady state. Our main goal is to examine the evolution of the binary in systems where the binary is aligned with the disk plane, and systems where the binary and disk are misaligned. To do so, we examine the different contributions to the torque, and how this affects the long-term evolution of the binary. Our study is the first 3D hydrodynamic grid-based calculation of the entire circumbinary disk system, including the binary. Shi et al. (2012) have performed 3D MHD simulations of circumbinary disks but have excised the central region which includes the binary. As mentioned earlier others have previously studied misaligned disks but ours is the first study of the torques in such systems. We find that as in some previous studies, in the 2D case the binary gains angular momentum and the semi-major axis expands. More interestingly, in the both the aligned and misaligned 3D simulations the same is true - the binary grows, contrary to the common idea that the binary always loses angular momentum to the disk.

The paper is organized as follows. In Section 2, we summarize our numerical methods and analysis procedures. In Section 3, we present our results on disk properties, including the accretion rate and torques, for both aligned and misaligned disks. Finally, in Section 4 we summarize our findings and discuss the physical implications of the results.

2. Numerical Methods

2.1. Numerical Setup

We use athena++ (Stone et al., 2008) as the numerical tool to solve the equations of viscous hydrodynamics:

[TABLE]

where $\rho$ is the density, $\mathbf{v}$ the velocity, and $P\equiv\rho c_{\rm s}^{2}$ the pressure. We enforce a local isothermal equation of state 111We achieve the local isothermal with adiabatic equation of state of $\gamma=1.001$, together with instantaneous cooling enforced at each timestep. of the following form (Farris et al., 2014)

[TABLE]

where the disk aspect ratio $h\equiv H/R=0.1$ is fixed in this study. The binary potential is

[TABLE]

The total binary mass is $GM=GM_{1}+GM_{2}=1$, where subscript ‘1’ and ‘2’ denote the primary and the secondary respectively. The location of the binary members are $\mathbf{r_{1}}=(x_{1},y_{1},z_{1})$ and $\mathbf{r_{2}}=(x_{2},y_{2},z_{2})$ in Cartesian coordinates. For disk locations far from the binary $c_{\rm s}\propto R^{-1/2}$ where $R=\sqrt{x^{2}+y^{2}}$ is the cylindrical radius, while in the vicinity of the individual holes $c_{\rm s}\propto|\mathbf{r}-\mathbf{r_{i}}|^{-1/2}$ where $r=(x,y,z)$ denotes the spherical radius. In our 2D runs, the disk density $\rho$ is replaced by the surface density $\Sigma$, and the pressure is therefore the vertically integrated value.

In this work, we choose isotropic viscosity to account for the angular momentum transport within the disk. $\mathbf{T}_{\rm vis}$ is the viscous stress tensor, defined as

[TABLE]

in its component form, where $x_{i},x_{j},x_{k}\in[x,y,z]$. The isotropic viscous coefficient is defined following the $\alpha$-disk prescription (Shakura & Sunyaev, 1973) $\nu=\alpha c_{s}^{2}/\Omega$ with $\alpha=0.1$ and $\Omega=\Omega_{K}\equiv(GM/R^{3})^{1/2}$ where $R=\sqrt{x^{2}+y^{2}}$ the cylindrical radius. MML17 have carried out simulations using both $\alpha=0.1$ and $\alpha=0.05$ and find the same result for circular binaries, so we do not vary $\alpha$ in our study. Numerically, the stress tensor is calculated explicitly and added directly to the hydro momentum flux, while the binary forcing term is added to the right hand side of equation 2 as a source term.

The initial disk density profile is

[TABLE]

where

[TABLE]

is the midplane density. We set density unit $\rho_{0}=1$, binary separation $a=1$ as the length unit, and the time unit $\Omega_{\rm bin}\equiv\sqrt{GM/a^{3}}=1$. For our 3D simulations $p=1.5$ so that it matches $\Sigma\propto R^{-1/2}$ as in the 2D runs, where $p=0.5$. To facilitate a faster steady solution, we start with a gap interior $R\lesssim R_{s}=2.5a$ similar to the work of Tang et al. (2017). We correct the initial azimuthal velocity $v_{\phi}$ of the disk with the pressure effects from the local pressure profile and the binary’s quadruple potential term. The initial radial velocity $v_{R}$ takes into account of the small radial drift due to viscous accretion.

The binary is assumed to be of equal mass, and orbits on a fixed circular orbit with separation $a=1$. In both the 2D and 3D in-plane cases, the binary and the disk share the same orbital plane. In the 3D misaligned case, the binary is out of the $z=0$ disk plane with a fixed inclination angle $i=\pi/4$.

The computational domain extends to $[-10a,10a]\times[-10a,10a]\times[-5a,5a]$ in $x$, $y$, and $z$ respectively. We choose our boundary conditions along $x=\pm 10a$ and $y=\pm 10a$ direction to be outflow, and the vertical boundary conditions to be diode so that inflow is not allowed. A damping zone between $R=R_{i}=8a$ and $R_{o}=10a$ is imposed to quench the $m=4$ modes according to the following formula:

[TABLE]

where $u$ represents the hydro variables and $t_{\rm orb}=2\pi/\Omega_{K}$ is the orbital period at $R$.

We also simulate the accretion process with sink cells centered around each individual binary members with a spherical radius $r_{\rm sink}=0.1a$. Within this radius, we remove the mass according to the $\alpha$-disk model

[TABLE]

with the mass remove timescale

[TABLE]

This is the same sink prescription used by Tang et al. (2017). Our fiducial runs make use of the static mesh refinement of athena++. At the root level, the resolution is $256\times 256\times 128$ on the aforementioned $[-10a,10a]\times[-10a,10a]\times[-5a,5a]$ computational domain. We further make 2-level refinement for the regions $[-5a,5a]\times[-5a,5a]\times[-2.5a,2.5a]$ and $[-2.5a,2.5a]\times[-2.5a,2.5a]\times[-1.25a,1.25a]$. We therefore resolve the sink cell regions with roughly $10$ cells across.

2.2. Analysis Procedure

As in the work of MML17 and Tang et al. (2017), we would like to examine the effects of the torque on the evolution of the binary. To do so, we calculate the advective, viscous, gravitational, and total torques on the system. The definitions of each are as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The torques are calculated within the code at every timestep. Each component of the integral ($x,y,z$) is calculated separately over spherical shells throughout the disk. Each spherical shell has a width of two cells ($2\Delta x$), in order to get a smooth average. The torque data is output instantaneously every tenth of a binary period. In our calculations we use the data from the last 100 orbits in order to ensure we are analyzing the disk in a steady state.

Once calculated, these are used to examine if the binary is expanding or contracting. Though we have all three components of the torque, we consider the z-component of the torque only, as this is the contributing factor to the binary evolution. We examine the binary evolution using the following expression as described by MML17, which assumes a fixed circular orbit:

[TABLE]

Here, $l_{B}=[GMa(1-e_{B}^{2})]^{1/2}$ is the specific angular momentum of the binary, and $l_{0}$ is ${\langle\dot{J}\rangle}/{\langle\dot{M}\rangle}$ calculated from the average total torque described in equation 15. The accretion rate is calculated in spherical shells, described for the torque components, and then averaged over radius. $GM$ and $a$ are 1 in our simulations and the binary is on a fixed circular orbit so $e_{B}=0$, giving us a value for $l_{B}$ of 1. Here we are concerned with whether the binary is shrinking or growing, not the rate at which this is happening. This can be determined from the sign of $\dot{a}$, so the $\dot{M}/{M}$ term does not need to be calculated. From this, we see that the sign of $\dot{a}$ is determined by whether $l_{0}$ is less than or greater than $3l_{B}/8=3/8$.

In our misaligned disk simulations we calculate the inclination angle through the disk, to determine how well the disk is aligned with the binary orbital plane. We use the definition of inclination angle through the disk, $\delta$, as described by Fragner & Nelson (2010):

[TABLE]

where $\mathbf{J_{B}}$ is the angular momentum vector of the binary and $\mathbf{L}$ is the total angular momentum vector at some radius in the disk. As with the torques, $\mathbf{L}$ is calculated over spherical shells. We then use the inclination angle to determine the realignment time of the disk, i.e. the time it will take for the circumbinary disk to realign to the binary orbital plane. This is compared to the viscous time, which is calculated using the following:

[TABLE]

3. Results

In table 1 we show the different simulations we have run with their different properties. All simulations have the same disk aspect ratio ($h\equiv H/R=0.1$) and $\alpha$ ($\alpha=0.1$). The fiducial runs are runs A, E, and F. In runs B and C we test faster sink times for the sink cells representing the black holes. In run D we test a lower resolution to test convergence.

3.1. 2D Simulations

First we run 2D simulations to check our numerical setup and the reliability of the Cartesian grid. Many similar simulations have been done already, as seen in the work of MML17, Farris et al. (2014), Tang et al. (2017), and Muñoz et al. (2018). None use a Cartesian grid in their simulations. The work of Farris et al. (2014) do not examine the torques on the system though the other three do. In addition, the work of MML17 excise a central region so they do not resolve the binary itself.

3.1.1 General Properties

Our simulation (simulation A as described in table 1) has a resolution of $256\times 256$ at each level of mesh refinement, which gives us a resolution of about 102 cells per $a$ in the innermost region of $[-2.5a,2.5a]\times[-2.5a,2.5a]$. The viscous time at $r=a$ is 119 binary orbits, and the global viscous time of the disk is 5020 binary orbits. We run the simulation for 1050 binary orbits, about 1/5 of the global viscous time of the disk, in order to reach a steady state. This requires approximately 140 CPU hours to complete the simulation.

A snapshot of the surface density of the system at late time is shown in Figure 1. Only the inner part of the disk is shown. We can see clearly the circumbinary disk with a gap region cleared out around the central binary as expected, with a gap radius of roughly $r=2a-3a$. MML17 have the same size of the central gap region in their simulation with the same parameters as we have used. The gap is generally not centered about the central binary in both our simulation and those of MML17. Additionally, we see towards the top of the cavity in our simulation a further asymmetry in the gap region. A high density lump has formed in the disk at the edge of the gap. This lump is created from the material that is flung back towards the disk away from the binary. The lump orbits along with the disk, and is periodically torn apart and recreated as it orbits.

Around each black hole, we can also see the formation of mini disks, with accretion streams feeding the mini disks from the circumbinary disk. MML17 see the same features in their simulations (spiral density waves, accretion streams, and lump), though they are unable to see mini disks as they have excised the central region. The simulations of Tang et al. (2017) have the same disk aspect ratio and $\alpha$ values as our simulation, though they use a polar grid and a slightly different sink time. Their simulations are also consistent with ours as they see the same size asymmetric gap region, density lump, accretion streams, and spiral density waves. Tang et al. (2017) do resolve the inner portion of the disk, and also see mini disks that form around each black hole with accretion streams feeding them.

To analyze if the system has indeed reached a steady state we plot the integrated surface density profile of the disk as a function of radius and time, seen in Figure 2. This figure shows the last 100 binary orbits of the simulation. The main structure seen here is the spiral density waves, which are excited by the binary and travel outwards. The lack of secular trends is consistent with the disk having reached a steady state by the end of our simulation time.

3.1.2 Accretion Rate

The accretion rate is another indicator of whether we have reached a steady state. In Figure 3 we show the accretion rate through the disk. The upper plot shows the accretion rate through the disk for the last 100 orbits of the simulation, with each line representing an average over 10 orbits. The red dashed line indicates the average accretion rate, averaged over time and radius for these last 100 orbits. The lower plot is an average over the entire period of the last 100 orbits. The red dashed line is the same as in the upper plot - the average accretion rate averaged over both time and radius. The average rate is 0.012 $\Sigma_{0}\sqrt{GMa}$. This value is larger than for an $\alpha$-disk around a single black hole by a factor of 1.27. Farris et al. (2014) also find their disks to have larger accretion rates than single black hole disks. They note that this does not indicate that binaries enhance accretion, but rather that binaries are unable to fully clear the gap region and suppress accretion.

Although there is larger oscillation on shorter timescales as shown in the upper plot of Figure 3, the lower plot indicates that the average accretion rate is generally steady for the majority of the circumbinary disk. This behavior is consistent with the simulations in MML17. The most variation occurs around about $r=3a$, as seen in more detail in the upper plot. This region is near where the edge of the gap region is located, so there is interplay between material being accreted to the mini disk and some material being flung outwards due to the torque. This variation is the same as we see in Figure 2, which indicates this is where the inner edge of the disk oscillates in and out. We conclude that mass is flowing through the disk at a roughly steady rate.

We also record the amount of mass removed over time by the sink cells representing the central black holes. In doing so, we find that a single black hole has an accretion rate of 0.00559 $\Sigma_{0}\sqrt{GMa}$ giving a total accretion rate of 0.0112 $\Sigma_{0}\sqrt{GMa}$. Comparing to our calculation of accretion rate through the disk we find the two calculations agree to within 7%. The accretion rate calculated by the sink cells is shown in Figure 4. In the left-hand panel the most obvious feature is the burst-like variability on the timescale of about 5 binary orbits. MML17 see the same variability with the same frequency. They note that their power spectrum peaks at about $\Omega_{B}/5$ ($\omega=0.2$). We plot the power spectrum of the accretion rates in the right-hand panel of Figure 4. Our largest peak occurs at 0.2075, nearly the same value. This burst-like variation is associated with presence of the high density lump at the edge of the disk gap. There are other periodic variations with higher frequencies in the left-hand panel, reflected in the many other peaks in the power spectrum. These peaks are seen by MML17 in their simulations of circular binaries. Farris et al. (2014) have also performed similar 2D simulations for circular binaries (though unlike MML17 they do resolve the inner cavity) and find very similar behavior in the accretion rate. Their results for an equal mass ratio also show that the dominant peak occurs around the same frequency of 0.2.

3.1.3 Torques

In Figure 5, we show the time-averaged advective, viscous, and gravitational torques, as well as the total net torque. All torque components are normalized by the average accretion rate, 0.0112 $\Sigma_{0}\sqrt{GMa}$. Our results are very similar to what is seen in Figure 17 of MML17. The main difference occurs after $r=8a$ as our disk has a buffer zone whereas their disk extends much further out. We focus our comparison on our results outside $r=1a$ as their work has excised the central region. All torque components in our simulation have the same profile as in MML17, and the features in the profiles occur at the same values of $r$ and $\dot{J}/\dot{M}$. Our total net torque (the red dashed line in Figure 5) is remarkably flat with the largest variation at about $r=0.5a$. This is due to this radius being the location of the binary. The overall average value of the normalized net torque for our 2D simulation is 0.723 $\sqrt{GMa}$, which is close to the value of about 0.8 stated in MML17. The simulation of a circular equal mass binary in the work of Muñoz et al. (2018) (shown in their Figure 3) shows almost identical behavior to both our simulation and that of their own previous work in MML17. The value of their normalized torque is 0.68, in good agreement with our result.

As we see the same positive net torque as MML17 and Muñoz et al. (2018) with the same shape and value, evolution of the binary orbit will result in binary expansion. By taking the average of the total torque curve, we now have our $l_{0}$ value needed to utilize equation 16. As $0.723>3l_{B}/8$ we have a positive $\dot{a}$, indicating that the binary will be expanding, the same result as found in MML17 and Muñoz et al. (2018).

As we have noted in Section 3.1.1 the work of Tang et al. (2017) is very similar to the 2D simulations we have run. They find a dependence of the torque on the sink time, with rapid sink times resulting in positive torques and slower sink times resulting in negative torques. Though we run a slow sink time, our simulation therefore seems to correspond with their rapid sink time results as we find positive torques. We have varied our sink time to explore the effect on the torques and have not found the same result as Tang et al. (2017), we instead find that the torques in our simulations are not sensitive to the sink time. This is described in more detail in Section 3.1.5. The differences between our work and theirs may be related to the numerical methods used (we use a grid-based code and they use a moving-mesh code) or the initial density profile, which is slightly different. We have carried out a convergence study to make sure our results are not affected by numerical diffusion, described in the following section (Section 3.1.4).

3.1.4 Resolution Test

We have also performed a 2D simulation with resolution of $128\times 128$, keeping all other simulation properties identical to the $256\times 256$ simulation (simulation D as described in table 1). The surface density for the same time in the simulation as in Figure 1 is shown in Figure 6. From this figure we can see the same general structure of the system: the gap region within $r=2.5a$, spiral density waves, mini disks around the individual black holes, and accretion streams. We calculate accretion rate and normalized torque for this lower resolution simulation. The two 2D simulations show the same overall accretion rate through the disk. This is shown in Figure 7, which is the equivalent of the lower panel of Figure 3 but for this lower resolution run. Though there is variation in the accretion rate over time in both simulations, the shape of the accretion rate as a function of radius is the same, with lower variation in the $256\times 256$ simulation than the $128\times 128$ simulation. We find the same trend when we analyze the torque in the lower resolution simulation. The torques have the same numerical value, and the curves as a function of radius have the same behavior as the higher resolution run. Again, we do see slightly larger variation in each torque component curve in the lower resolution run. We conclude that the simulation is converged based on the comparison between the two resolutions showing the same results with very minor differences, which are small variations on the same results.

3.1.5 Comparing Sink Times

As stated in section 2.1 for this 2D simulation we have used the viscous time at the sink radius as the mass removal timescale (the sink time), which is $t_{rm}=29.814$. We now test different sink times to see if there is a relationship between the resulting torques and the mass removal rate. In Figure 8 we show the torque components and net torque for three different sink times, $t_{rm}=29.814$, $t_{rm}=5.0$, and$t_{rm}=0.1$ (these are simulations A, B, and C as described in table 1). The overall behavior is the same in all the curves and the average net torque is within 4% of the simulation with $t_{rm}=29.814$ for both of the other sink times. As we see no discernible difference in the sink time we use the viscous time ($t_{rm}=29.814$) for the 3D simulations that follow.

3.2. 3D Simulations - aligned disk

To examine if the same effect on binary evolution is seen in 3D we first perform 3D simulations with the binary orbiting in the same plane as the circumbinary disk, to test if the assumptions made in 2D simulations are valid. Some work in 3D has been done by Shi et al. (2012), though they perform full MHD simulations rather than pure hydro simulations. Our simulation is simulation E as described in table 1. As mentioned in Section 2.1 our resolution is $256\times 256\times 128$ on a computational domain of $[-10a,10a]\times[-10a,10a]\times[-5a,5a]$ in $x$, $y$, and $z$, with two levels of mesh refinement. In order to compare with our 2D simulations, we run the simulations for the same amount of time - 1050 binary orbits, equivalent to about 1/5 of the global viscous time of the disk. This requires about 106000 CPU hours to complete.

3.2.1 General Properties

Looking at the general properties of the 3D simulation, we see very similar features to the 2D disk as expected for the aligned case. Figure 9 shows isosurfaces of density, with isosurface chosen to particularly illustrate the orientation of the mini disks and circumbinary disk, as well as the high density features: the mini disks, high density lump near the edge of the gap region, and the beginning part of the accretion stream. Because we use isosurfaces the gap region looks to be devoid of gas. This is simply an artifact of the isosurfaces. The structure of the disk looks very similar to the structure of the 2D disk. We examine the surface density to see if we are in a steady state. The surface density is integrated over spherical shells, and we plot a spacetime diagram of surface density in Figure 10 similar to Figure 2 for the 2D simulation. The structure is nearly identical to the 2D case; we see spiral density waves and not much other structure. From this we expect the accretion rate and torques to behave similarly to the 2D case.

3.2.2 Accretion Rate

The accretion rate through the 3D aligned disk is shown in Figure 11. As in Figure 3 the upper plot is the rate through the disk for the last 100 binary orbits, each line an average over 10 orbits. The lower plot is the average accretion rate over this entire period. The red dashed line in both cases is the overall average of the accretion rate, averaged over radius and time. Again we see larger oscillation on shorter timescales within the gap region and inner edge of the disk, whereas the average over the last time period (the right-hand plot of Figure 11) is steadier. From Figure 10 we see that again this variation is associated with the oscillation of the inner edge of the disk. The value of the accretion rate is lower than the 2D case by roughly a factor of of 5. When normalized by the density of the disk, the normalized accretion rate of the 3D aligned disk agrees with the normalized accretion rate of the 2D disk.

As with the 2D simulation, we compare our calculation of accretion rate through the disk with the accretion rate as measured by mass removal by the sink cells representing the black holes. Each black hole has an accretion rate of 0.001195 $\Sigma_{0}\sqrt{GMa}$, giving a total accretion rate of the binary of 0.00239 $\Sigma_{0}\sqrt{GMa}$. This agrees with the disk accretion rate to 8%, about the same accuracy as we found with the 2D simulation. The individual accretion rates can be seen in the left-hand panel of Figure 12, which shows the accretion rates over ten binary orbits. There is an overall variability on the timescale of roughly 5 or 6 binary orbits, similar to what we saw in the 2D simulation. The right-hand panel of Figure 12 shows the power spectrum of the accretion rate for both the primary and secondary. There is a clear peak at $\omega=0.2245$, which is roughly the same frequency fluctuation discussed in Section 3.1.2 and seen in MML17 ($\omega=\Omega_{B}/5$). Comparing the 2D (Figure 4) and 3D simulations, we see that in the 3D simulation we do see the high frequency fluctuations, but the overall power of the 3D simulations is much lower than the 2D simulations.

3.2.3 Torques

Figure 13 shows the z-component of the torque for the aligned simulation, as we focus on the evolution of the binary which is fixed to be in the same plane as the disk midplane ($z=0$). The gravitational and viscous torque have much of the same shape and behavior that we see in the 2D simulation, shown in Figure 5. The main difference between the 2D case and the 3D aligned case lies in the advective torque, in the outer regions of the disk. There is more variation in the 3D aligned case in the region $r=4a-6a$. Additionally, Figure 5 shows that the advective and viscous torques have the same slope in the region outside roughly $r=4$ for the 2D simulation, whereas Figure 13 shows a slightly shallower slope for the advective torque compared to the viscous torque. In this outer region, the disk is thicker which means movement in the z direction is more likely to occur than in the inner regions of the disk. This increases the x- and y-components of the advective torque, and in turn also perhaps decreasing the z-component of the torque. The disk holds equatorial symmetry, as indicated by the average z-momentum of the disk. This value varies through the simulation but is always less than 0.0006 $\Sigma_{0}\sqrt{GM/a}$, which is roughly 1/1000 of the x- and y-momenta.

The total torque is still very close to constant throughout the disk, with an average value of 0.558 $\sqrt{GMa}$. This value is smaller than those found in 2D simulations, both our own and those of MML17 and Muñoz et al. (2018). We apply the same analysis as in Section 3.1.3 to determine the evolution of the binary. We find a positive $\dot{a}$ in the evolution equation (equation 16). We conclude that in 3D when the binary and the disk are in the same plane, we will see expansion of the binary over time.

3.3. 3D Simulations - misaligned disk

Now that we have examined systems with the binary in the same plane as the disk, we perform simulations of misaligned systems. The disk is initialized in the $z=0$ plane, and the inclination of the binary orbital plane is changed to a fixed value. We focus on the simulation where the binary inclination is 45°(simulation F as described in table 1). The simulation has the same computational domain as the aligned disk case ($[-10a,10a]\times[-10a,10a]\times[-5a,5a]$) with the same resolution ($256\times 256\times 128$). The total run time remains the same (1050 binary orbits), and the required CPU time is 89200 hours.

3.3.1 General Properties

In Figure 14, we show density halfway through the simulation to illustrate the structure of the misaligned disk. At the very beginning of the simulation, the disk is aligned in the $z=0$ plane. Very quickly the disk begins to warp, realigning itself with the orbital plane of the binary. The inner regions begin to warp first, with the warp propagating outward in the disk slowly. The very outer edges of the disk remain mostly aligned with the $z=0$ plane for much of the simulation. Mini disks again form around each individual black hole, but now the orientation of the mini disk is also variable. When the mini disks are formed they are aligned with the plane of the disk. Throughout the simulation, the disks’ orientation changes to realign with the binary orbital plane, which occurs much more quickly than the realignment of the circumbinary disk. Figure 14 illustrates all these properties: the outer portions of the disk are aligned with the $z=0$ plane, the mini disks are aligned with the binary orbital plane at 45°, and the inner portion of the disk has started to align with the binary orbital plane but has not yet reached full alignment. As stated in Section 3.2.1 the apparent lack of gas in the gap region is an artifact of the fact that we are using isosurfaces of density, which were chosen to illustrate high density features and the orientation of the circumbinary and mini disks.

Nixon et al. (2013) perform 3D simulations around binary black holes, using a smoothed particle hydrodynamics code. They study a range of inclinations for equal mass binaries on circular orbits. In their simulations they see tearing of the disk for misaligned disks. The inner region of our disks show some disruption, though not as dramatic tearing as what Nixon et al. (2013) see. To estimate the breaking radius of the disk, i.e. where the disk will tear. Nixon et al. (2013) give the condition that the disk will break when the viscous torque in the disk is less than the precession torque. Using the analytic expressions for both (they do not calculate torques within their simulations), they give the disk breaking radius with the following expression:

[TABLE]

They use a disk aspect ratio of $h=0.01$, which indicates that their entire disk (which extends to $r=8a$) will break. Our simulation uses an aspect ratio of $h=0.1$ which gives a breaking radius of about $3.5a$. Within this radius we primarily have the gap region and the very inner edge of the disk, so most of our disk will remain intact and not tear, hence why we see much less dramatic tearing than Nixon et al. (2013).

3.3.2 Accretion Rate

The accretion rate through the disk is shown in Figure 15. In the upper panel we have the accretion rate plotted in averages over 10 orbits for the last 100 orbits of the simulation, and in the lower panel we have the average accretion rate over that entire 100 orbit period. We see similar behavior in this misaligned case that we saw previously in the aligned case. One difference is the value of the average accretion rate throughout the disk, 0.00899 $\Sigma_{0}\sqrt{GMa}$, is larger than the 3D aligned simulation, though the two rates are still within an order of magnitude. Nixon et al. (2013) find that their misaligned disks have much higher accretion rates compared to the aligned disk, i.e. tearing of the disk enhances the accretion rate. Their 45°inclined simulation has an accretion rate roughly 50 times the aligned simulation accretion rate. While our simulations also show an increased accretion rate in the 45°case, it is slightly less than a factor of four larger than the aligned case.

Figure 16 shows the accretion rate as measured by the mass removal of the sink cells. Both plots look very similar to the 3D aligned case - the frequency spectrum peaks at $\omega=0.229$ indicating the existence of a density lump in this misaligned case. The misaligned case shows more variation at higher frequencies than the aligned case.

3.3.3 Torques

In order to examine the torques for the misaligned system, we calculate the torque components in the same way as done previously, and then project the components into the binary orbital frame. We then take the z-components of the torques in the binary frame, and perform the same analysis to determine the binary evolution. The z-components in the binary frame are shown in Figure 17, along with the total torque. We can see that the shapes of the components have vastly different form compared to the in-plane case. As our disk is misaligned, the torque is no longer dominated by the z-component only.

We see from Figure 17 that the total torque is remarkably constant throughout the disk. Applying the analysis described in Section 2.2, we calculate an $l_{0}$ value of 4.883, which is greater than $3l_{B}/8$, again indicating expansion of the binary orbit. We conclude that in all cases, a viscous disk will cause the binary orbit to expand.

3.3.4 Disk Realignment

We expect the disk to realign with the plane of the binary. In Figure 18 we show the inclination angle of the disk as a function of radius, where the angle is the relative angle between the disk and the binary orbital plane. The different curves show inclination as a function of radius for different times in the second half of the simulation. The disk realignment is clearly seen from the movement of these different curves towards lower relative inclination angle. We can also see that the inner portion of the disk is realigning more quickly than the outer portion of the disk. One significant feature is that the buffer zone seems to be acting as an anchor preventing the very outer region of the disk from realigning. This could potentially be affecting the inner portions of the disk, and future work should include a much larger disk or different treatment of the buffer zone to limit this effect.

In Figure 19 we show the inclination angle of the disk relative to the binary orbital plane over time for half of the total simulation time but only in the outer region of the disk. This region is from $r=7a$ to $r=8a$ (before the buffer zone begins) in order for us to examine the global evolution of the disk. We perform a least squares regression fit to this curve to calculate the realignment time, and find a time of $\sim$ 4134 binary orbits. The global viscous time of the disk (calculated for $r=10a$) is $\sim$ 5020 binary orbits, so $t_{realign}/t_{visc}\approx 0.8236$. If we compare the realignment time to the viscous time at $r=7.5a$, we find $t_{realign}/t_{visc7.5}\approx 1.271$. This result is in good agreement with the work of Fragner & Nelson (2010), which finds the realignment time to be roughly on the timescale of the viscous evolution. Foucart & Lai (2014) study the evolution of warp in accretion disks analytically, and apply their findings to circumbinary disks. They find that circumbinary disks should realign on the global precession timescale, which is much shorter than the viscous timescale. Their surface density profile that is steeper than we use, and they note that shallower profiles make it more difficult for torques applied on the smaller amount of material at small radii to affect the global evolution of the disk. Shallower surface density profiles therefore will have longer realignment timescales (than the precession time) which is consistent with our results.

In the inner regions of the disk we perform the same analysis, and find as expected the inner regions realign more quickly than the outer regions, but are still within the same magnitude of the viscous evolution time.

4. Summary and Discussion

4.1. Summary

We have performed hydrodynamic simulations of circumbinary accretion disks using the grid-based finite volume code athena++. The binary is on a fixed circular orbit with an equal mass ratio, and individual members of the binary are treated as sink cells. The disks are locally isothermal and use an $\alpha$ viscosity. We performed 2D simulations, and 3D simulations for both aligned and misaligned disks and study the mass accretion rate and the torques on the system.

The mass accretion rate in all three simulations shows steady accretion with variation on short timescales which is associated with a high density lump on the inner edge of the disk. The 2D and 3D aligned disk simulations have comparable accretion rates though the 3D rate is roughly a factor of 5 smaller than the 2D rate. This is due to the density profiles through the disk. As mentioned in Section 2.1 the surface density of both disks has the same scaling but the full density profile scales differently, resulting in slightly different overall surface density profiles and therefore slightly different values of the accretion rate. The 3D misaligned disk shows an enhanced accretion rate compared to the aligned disk, possibly due to some tearing of the disk.

The net torque on the system indicates in all simulations that the binary orbit will expand. This confirms some previous results in 2D. The different components of torque show the same behavior in 2D and the 3D aligned disk. The 3D misaligned disk shows very different behavior in the individual components of the torque, but as stated the net torque is still positive and indicates expansion of the binary.

The overall sign of the net torque is very sensitive to the value of the individual torque components, in particular the advective torque. This indicates that a careful examination of the value of $\alpha$ is needed. As mentioned briefly in Section 2.1 MML17 carry out simulations of circular binaries using $\alpha=0.1$ and $\alpha=0.05$ and find the same result for both values of $\alpha$. In recent MHD simulations of accretion disks it is found that the equivalent $\alpha$ that results from turbulence in such disks is around 0.1 so this value is reasonable (Zhu & Stone, 2018).

4.2. Implications for Binary and Disk Evolution

It is a common assumption that gaseous disks are the mechanism by which supermassive black hole binaries are able to shrink their semi-major axes to less than a parsec (Begelman et al., 1980; Escala et al., 2005). As discussed in Section 2.2, we analyze the torques in the simulation to determine if the binary will expand or shrink. In all of our simulations, both 2D and 3D, we find that the torques on the system cause the binary to grow, contrary to some previous results. This suggest that hydrodynamic circumbinary disks cannot be a solution to the final parsec problem for supermassive black hole binaries, at least for equal mass binaries. Perhaps different configurations of the disk may provide a solution, such as different mass ratios, eccentricities, or more massive disks. For equal mass binaries however, it seems that circumbinary accretion disks will prevent the merger of these binaries.

To examine this more closely we compare the rate of change of the semimajor axis due to the gas disk (which we find widens the binary) with the rate of change due to dynamical friction from stars (which shrinks the binary). To calculate the rate of change from the gas disk we use equation 16 using typical values for a binary system. We use a mass of $10^{6}M_{\odot}$ for each black hole and a binary separation ($a$) of one parsec. The $\dot{M}$ and $l_{0}$ values are from our 3D aligned simulation (simulation E in Table 1), which are 0.0026 $\Sigma_{0}\sqrt{GMa}$ and 0.558 $\sqrt{GMa}$ respectively. To calculate the rate of change due to stars we use equation 8.20 from Binney & Tremaine (2008):

[TABLE]

or the following which is in a more directly comparable form to equation 16

[TABLE]

$\rho$ is the density of stars and $\sigma$ is the velocity distribution of the stars. We choose a typical stellar density of $3.386\times 10^{-24}$ g/cm3 which corresponds to 0.1 stars per cubic parsec and average stellar mass of 0.5$M_{\odot}$.The velocity dispersion is set to 75 km/s. Using these parameters, we find that the two rates are comparable if the mass of the disk is around 0.04 times the combined mass of the black holes. For the case where the mass of the disk is 0.1 the mass of the black holes, $({da}/{dt})_{\textrm{gas}}=11.24$ cm/s and $({da}/{dt})_{\textrm{stars}}=4.100$ cm/s. This is significant as it indicates that the circumbinary disk can potentially slow or stall the shrinking of the binary.

We found the length of the disk realignment time to be on the same timescale as the global viscous time of the disk, in good agreement with the results of Fragner & Nelson (2010). Their disk configuration has the primary at the center of the disk, with the secondary at a distance of $30a$ from the primary. This leads us to conclude that the alignment evolution of the disk in any configuration will be on the viscous evolution time.

4.3. Limitations and Future Work

One major limitation of this work is the absence of magnetic fields. We use the $\alpha$-disk prescription for viscosity, and in reality we expect the disk viscosity to be driven by turbulence, resulting from the magnetorotational instability (Balbus & Hawley, 1991). Some magnetohydrodynamic simulations have already been performed, by Shi et al. (2012). They perform 3D simulations of an aligned disk around an equal mass binary on a circular orbit, and find that the binary shrinks. However, they cut out the central region of the simulation which includes the binary. In future work, magnetic fields should be included and full 3D simulations for a range in binary-disk inclinations explored. Given that various results indicate the importance of the advective torque to the evolution of the binary, it is even more important that MHD calculations be performed.

In addition, several assumptions made should be relaxed or explored further in the future. This includes the binary mass ratio, disk aspect ratio, and disk mass. We also use a local isothermal equation of state, so in future work this should change to include a more detailed treatment of the thermodynamics to include the cooling time of the gas. Though we assume the binary orbit to be fixed as the evolution of the binary does not change the eccentricity significantly, a careful treatment would include coupling the evolution of the binary with the evolution of the disk, so the semi-major axis and eccentricity are allowed to change as the simulation runs. Additionally, we do not include the self-gravity of the disk or general relativistic effects, which could significantly affect the dynamics of the binary-disk interaction.

We thank Ryan Miranda, Diego Muñoz, Dong Lai, Zoltán Haiman, Alwin Mao, Goni Halevi, and Cole Holcomb for helpful input and discussions. We also thank the referee for helpful comments that have improved our paper.