FORMATION OF CIRCUMBINARY PLANETS IN A DEAD ZONE

Previous work on circumbinary accretion disks has assumed the disk to be fully ionized and turbulent. We follow Alexander (2012) with the use of a one-dimensional circumbinary disk model but incorporate a layered disk model (Martin & Lubow 2011, 2013; Lubow & Martin 2012, 2013). The disk is thermally ionized, and MRI active throughout, if the midplane temperature is greater than some critical value, T_c > T_crit. The value of the critical temperature is thought to be around 800 K (Umebayashi 1983), but we find that this temperature is not reached within the circumbinary disk. Below this temperature, the disk may become layered. The surface layers of the disk (with maximum surface density Σ_crit) are ionized by cosmic rays or X-rays from the central stars. If the total surface density is greater than the critical value, Σ > Σ_crit, then the active layers have surface density Σ_m = Σ_crit and the midplane contains a dead zone with mass Σ_d = Σ  −  Σ_crit. If Σ < Σ_crit, then the whole disk is externally ionized and MRI active, Σ_m = Σ. This dead zone may become self gravitating if sufficient mass builds up such that the Toomre (1964) parameter, Q, is less than a critical value that we take to be 2. However, in this work, we find that there is not sufficient build up in a circumbinary disk for self-gravity.

The structure of the MRI active layer and the dead zone has not yet been well constrained. The surface density that is ionized by cosmic rays or X-rays, Σ_crit, has been calculated from the ionization balance of the external effects and internal effects such as Ohmic and ambipolar diffusion (e.g., Bai 2011; Simon et al. 2013). However, these calculations find accretion rates that are much lower than those observed in T Tauri stars, which require Σ_crit > 10 g cm⁻² (e.g., Perez-Becker & Chiang 2011). With these uncertainties in mind, we allow Σ_crit to be a free parameter (e.g., Armitage et al. 2001; Zhu et al. 2009, 2010) and consider different values.

Material orbits the central binary of mass M = M₁ + M₂, at radius R, at Keplerian angular velocity, $\Omega =\sqrt{GM/R^3}$. The governing accretion disk equation, which includes a binary torque term, for the evolution of the total surface density Σ = Σ_m + Σ_d, and time t is

$$\begin{eqnarray} \frac{\partial \Sigma }{\partial t} &=& \frac{1}{R}\frac{\partial }{\partial R} \left \lbrace 3R^{1/2}\frac{\partial }{\partial R}[R^{1/2}(\nu _{\rm m}\Sigma _{\rm m}+\nu _{\rm d}\Sigma _{\rm d})]\right. \nonumber\\ &&\left. -\frac{2\Lambda \Sigma R^{3/2}}{(GM)^{1/2}}\right \rbrace \end{eqnarray} \tag{ 1 }$$

(Pringle 1981; Lin & Papaloizou 1986). The viscosity in the MRI active surface layers is

$$\begin{equation} \nu _{\rm m}=\alpha _{\rm m}\frac{c_{\rm m}^2}{\Omega }, \end{equation} \tag{ 2 }$$

where α_m = 0.01 (e.g., Hartmann et al. 1998) is the Shakura & Sunyaev (1973) viscosity parameter, $c_{\rm m}=\sqrt{{\cal R} T_{\rm m}/\mu }$ is the sound speed, T_m is the temperature in the layer, ${\cal R}$ is the gas constant, and μ is the gas mean molecular weight. Similarly, the viscosity in the dead zone layer is

$$\begin{equation} \nu _{\rm d}=(\alpha _{\rm d}+\alpha _{\rm g})\frac{c_{\rm s}^2}{\Omega }, \end{equation} \tag{ 3 }$$

where $c_{\rm s}=\sqrt{{\cal R}T_{\rm c}/\mu }$ is the sound speed at the midplane. The α_g term due to self-gravity is zero unless Q < Q_crit (see Martin & Lubow 2011, for more details). MHD simulations suggest that there may be some residual viscosity in the dead zone and this is parameterized with α_d but its value is still undetermined (e.g., Fleming & Stone 2003; Simon et al. 2011). We take α_d = 0 in most of our models but consider one with a high value of α_d = 10⁻⁴ for comparison.

The tidal torque from the binary is

$$\begin{equation} \Lambda (R,a)=\frac{q^2 GM}{2R}\left(\frac{a}{\Delta _{\rm p}}\right)^4 \end{equation} \tag{ 4 }$$

(Armitage et al. 2002), where Δ_p = max(H, |R − a|), H = c_s/Ω is the disk scale height, a is the binary separation, and q = M₂/M₁ is the mass ratio of the stars. Coupled with Equation (1) we solve a simplified energy equation

$$\begin{equation} \frac{\partial T_{\rm c}}{\partial t}=\frac{2(Q_+-Q_-)}{c_{\rm p}\Sigma } \end{equation} \tag{ 5 }$$

(Pringle et al. 1986; Cannizzo 1993). The disk specific heat is $c_{\rm p}=2.7 {\cal R}/\mu$. The local heating is

$$\begin{equation} Q_+=Q_{\nu }+Q_{\rm tid}, \end{equation} \tag{ 6 }$$

where the viscous heating is

$$\begin{equation} Q_{\nu }=\frac{9}{8}\Omega ^2(\nu _{\rm m}\Sigma _{\rm m}+\nu _{\rm d}\Sigma _{\rm d}) \end{equation} \tag{ 7 }$$

and the tidal heating term is

$$\begin{equation} Q_{\rm tid}=(\Omega _{\rm b} -\Omega) \Lambda \Sigma \end{equation} \tag{ 8 }$$

(e.g., Lodato et al. 2009; Alexander et al. 2011), where the orbital frequency of the binary is $\Omega _{\rm b}=\sqrt{GM/a^3}$. We also include a simple prescription for the effects of irradiation from the central stars. We approximate the irradiation temperature of the disk with

$$\begin{equation} T_{\rm irr} = T_\star \left(\frac{2}{3 \pi }\right)^\frac{1}{4}\left(\frac{R_\star }{R}\right)^\frac{3}{4} \end{equation} \tag{ 9 }$$

(Chiang & Goldreich 1997), where we take T_⋆ = 4000 K and R_⋆ = 2 R_☉. If the midplane temperature of the disk drops below the irradiation temperature, then we assume that the disk is isothermal and T_c = T_m = T_e = T_irr (see also Martin et al. 2012b). The cooling rate is

$$\begin{equation} Q_-=\sigma T_{\rm e}^4, \end{equation} \tag{ 10 }$$

where T_e is the surface temperature of the disk and σ is the Stefan–Boltzmann constant. The temperature at the midplane, T_c, is related to the temperatures in the active layers, T_m, and at the surface, T_e, by considering the energy balance in the layered model above the midplane and is described with Equations (6)–(12) in Martin et al. (2012b).

The binary torque clears the inner regions of the binary orbit (Artymowicz & Lubow 1994, 1996). It is strong enough to prevent all accretion flow into the binary region. Hence, the mass of the disk remains constant in time. However, some accretion may occur on to the binary through tidal streams. The maximum average of this is thought to be about 10% of the steady accretion flow (MacFadyen & Milosavljević 2008). Hence, we also include a model with accretion on to the binary stars. Material is removed from the inner regions of the disk at a fraction of the steady accretion rate 3π(ν_mΣ_m + ν_gΣ).

The initial surface density is chosen to be

$$\begin{equation} \Sigma =\Sigma _0 \left(\frac{R_0}{R}\right) \exp \left(-\frac{R}{R_0}\right),\quad R>10a, \end{equation} \tag{ 11 }$$

where we take R₀ = 15 a and scale the constant Σ₀ such that the total mass of the disk, M_disk, is equal to some constant in the range 0.01–0.1 M_☉ (see also Alexander 2012). This initial condition has little effect on the structure of the disk at later times. The total mass of the disk, not the initial distribution of mass, is the important feature.

We take a grid of 200 points distributed evenly in log R from R_in = 2a up to R_out = 1000 AU. At the inner edge, we have a zero surface density inner boundary condition, but because the binary torque is strong, there is no material here. At the outer edge, we take a zero radial velocity boundary condition to prevent mass loss, but we have chosen a radius large enough that no significant amount of material spreads this far in a reasonable timescale. In the following section, we describe the results of our simulations.

Figure 1 shows the disk structure for two circumbinary disk models with central binary star masses of M₁ = 1 M_☉ and M₂ = 1 M_☉ with a separation a = 0.2 AU, accretion efficiency = 0 and total disk mass M_disk = 0.05 M_☉. Model R1 is fully turbulent and model R2 has a dead zone defined by Σ_crit = 20 g cm⁻². In both cases, the inner parts of binary region are cleared out to a radius of greater than 1 AU. In model R2, the dead zone structure is set up in less than 10⁵ yr and then is fairly constant in time. The disk is fully MRI turbulent in the inner and outer parts where the surface density is less than the critical that is ionized by external sources, Σ < Σ_crit. The dead zone extends from 1.91 AU to 18.5 AU at a time of 3 Myr, for example, and is not self gravitating at any time. The mass of the dead zone decreases in time as the disk spreads out in radius. The peak in the surface density distribution is much farther from the binary star in the disk with a dead zone, than in the fully turbulent case where the peak is close to the inner edge of the disk. The radius of the peak surface density moves inward in time. We discuss this feature further in Section 2.3.

Zoom In Zoom Out Reset image size

Similarly, the peak in the temperature distribution of the fully turbulent disk is closer to the binary star than the disk with a dead zone. The snow line radius, R_snow, is the radius beyond which ice condenses from the disk, which occurs at a temperature of around T_snow = 170 K (e.g., Lecar et al. 2006; Martin & Livio 2012, 2013). Thus, the snow line temperature occurs at a radius farther from the binary in the disk with a dead zone. Because the solid mass density is much larger beyond the snow line, planetesimal formation is much easier there and thus giant planets may form more easily. Hence, we expect both gas giants and rocky terrestrial planets to form in circumbinary disks, as they do in circumstellar disks. In the dead zone model, terrestrial planet formation is easier compared with the turbulent model as there is a larger radius out to which they can form.

In Table 2, we show models with various disk parameters and find that the dead zone is extensive for a range of disk masses and critical surface densities. The smaller the critical surface density, Σ_crit, the larger the dead zone and the larger the radius of the peak surface density (comparing models R2 and R4, for example). If the disk mass is small (less than 0.05 M_☉), then the critical surface density must be small, (Σ_crit < 20 g cm⁻²) for a dead zone to exist for the whole disk lifetime, up to 10 Myr (see models R6 and R7). On the other hand, if the disk is massive, then the dead zone is more extensive and the peak surface density is at a larger radius (see model R8). The accretion flow from the disk on to the binary (model R5 with > 0) has little effect on the structure of the disk and the dead zone. In this case, the average accretion rate on to the binary stars is 6.5 × 10⁻¹⁰ M_☉ yr⁻¹.

Table 2. Summary of Circumbinary Disks Models

Model	M1	M2	a	Σcrit	αd		Mdisk	MDZ	Rd1	Rd2	Rpeak	Rpeak	Mdisk	MDZ	Rsnow
	(M☉)	(M☉)	(AU)	(g cm−2)			(M☉)	(M☉)	(AU)	(AU)	(AU)	(AU)	(M☉)	(M☉)	(AU)
Time/Myr	all						0	3				10
R1	1.0	1.0	0.2	>44	0.0	0.0	0.05	...	...	...	2.72	2.94	0.05	...	2.94
R2	1.0	1.0	0.2	20	0.0	0.0	0.05	0.0231	1.91	18.5	5.72	4.70	0.05	0.00575	5.95
R3	1.0	1.0	0.2	20	10-4	0.0	0.05	0.0179	1.77	20.8	4.70	4.35	0.05	0.00349	5.29
R4	1.0	1.0	0.2	10	0.0	0.0	0.05	0.0369	1.91	21.6	7.82	7.52	0.05	0.0249	9.15
R5	1.0	1.0	0.2	20	0.0	0.1	0.05	0.0202	1.91	18.5	5.50	4.35	0.0435	0.00243	5.09
R6	1.0	1.0	0.2	20	0.0	0.0	0.02	0.00208	1.99	11.6	4.35	3.18	0.02	...	...
R7	1.0	1.0	0.2	10	0.0	0.0	0.02	0.00994	1.91	17.1	6.69	5.72	0.02	0.00328	5.72
R8	1.0	1.0	0.2	20	0.0	0.0	0.1	0.0669	1.70	22.5	6.43	5.94	0.1	0.0381	9.15
R9	0.7	0.2	0.2	20	0.0	0.0	0.05	0.0362	1.06	12.0	4.02	3.72	0.05	0.0249	5.72
R10	0.7	0.2	0.2	10	0.0	0.0	0.05	0.0435	1.20	12.0	5.29	5.09	0.05	0.0375	6.43
R11	1.0	1.0	1.0	>13.3	0.0	0.0	0.05		...	...	12.5	12.9	0.05	...	...
R12	1.0	1.0	2.0	>5.8	0.0	0.0	0.05	...	...	...	24.7	25.4	0.05	...	...
R13	1.0	1.0	5.0	>2.0	0.0	0.0	0.05	...	...	...	60.3	63.1	0.05	...	...

Notes. Column 2 shows the mass of the primary star, Column 3 shows the mass of the secondary star and Column 4 shows the semi-major axis of the binary orbit. Column 5 shows the critical surface density that is ionized by cosmic rays or X-rays from the central binary and Column 6 shows the viscosity parameter in the dead zone. Column 7 shows the parameter that describes the efficiency of accretion on to the binary and Column 8 shows the initial mass of the disk. Columns 9–12 describe the disk structure at a time t = 3 Myr. Column 9 shows the mass of the dead zone, Columns 10 and 11 show the inner and outer radius of the dead zone, respectively and Column 12 shows the radius of the peak surface density. Columns 13–16 describe the disk structure at t = 10 Myr. Column 13 shows the radius of the peak surface density, Column 14 shows the disk mass, Column 15 shows the dead zone mass and Column 16 shows the radius of the snow line. In rows where there is no snow line radius shown, the snow line is inside the inner edge of the disk.

Download table as:  ASCIITypeset image

For a disk without a dead zone, the maximum surface density occurs close to the inner edge of the disk, as shown in the top plots of Figure 1. However, in a disk that contains a dead zone, the peak surface density may lie far from the inner edge (see Table 2).

The inner parts of the circumbinary disk show little radial movement over the lifetime of the disk (as shown in Figure 1). In such a disk the inward or outward movement of dust particles is determined by the sign of the radial pressure gradient (e.g., Takeuchi & Lin 2002, 2005). The pressure gradient, ∂P/∂R, where $P=c_{\rm s}^2\Sigma /\sqrt{2\pi }H$, is positive inside of the radius of the peak surface density, R_peak, and negative outside. This suggests that particles of all sizes drift toward and collect around this radius. The radius of the peak surface density occurs far from the inner edge of the disk for models with a dead zone. It moves outward with decreasing active layer surface density.

Paardekooper et al. (2012) found that planetesimal accretion can occur in R > 20a = 4.4 AU for Kepler-16, R > 12a = 2.64 AU for Kepler-34, and R > 15a = 2.7 AU for Kepler-35. The results of Rafikov (2013) suggest that these numbers may overstate the difficulty of forming planetesimals, at least while the gas disk remains relatively massive. Models R9 and R10 represent the specific example of Kepler-16. Models R1–R8 represent circumbinary disk models for Kepler-34 and Kepler-35. For reasonable dead zone models, we find that the peak surface density is farther out than the region where planetesimal accretion is inhibited. Thus, the dead zone is a likely formation site for the circumbinary planets. The binary perturbations do not affect planetesimal growth in a circumbinary disk with a dead zone.

A non-zero viscosity in the dead zone is uncertain, and thus we have assumed that the dead zone has zero turbulence, α_d = 0, in most of our models. However, shearing box simulations suggest that the MHD turbulence generated in the disk surface layers may produce some hydrodynamic turbulence in the dead zone layer that may produce a small but non-zero viscosity (e.g., Fleming & Stone 2003; Simon et al. 2011; Gressel et al. 2012).

With this uncertainty in mind, we also performed one simulation with a viscosity in the dead zone, model R3. We choose the parameters of model R2 but include a viscous term in the dead zone with α_d = 10⁻⁴. We note that this is likely to be an upper limit to the viscosity in the dead zone. In a circumstellar accretion disk, the amount of mass flow through the dead zone should be less than that through the active layer because the turbulence in the dead zone is generated by the turbulence in the active layers (see also Bae et al. 2013). Even in this extreme case, we find that the disk structure is not significantly altered. The extent of the dead zone is actually larger, but the mass in the dead zone is slightly smaller. The peak in the surface density occurs at a radius slightly closer to the binary stars, but still farther out than the region that is hostile to planetesimal formation for Kepler-34 and Kepler-35. If the magnitude of the viscosity in the dead zone is reduced by a factor of around 10 or more, then the effect on our results is minimal. Thus, we conclude that a small but non-zero viscosity in the dead zone will not significantly affect the results presented here.

The table also includes some circumbinary disk models for larger separation binaries (models R11–R13). Here, the surface density of the disk (and thus the mass) must be much larger for a dead zone to be present. For an equal mass binary with a separation of 1 AU, the peak surface density is only 13.3 g cm⁻² (at a time of 3 Myr, see model R11) and this decreases with binary separation (see, for example, models R12 and R13). If a dead zone is required for planet formation, then this suggests that circumbinary planet formation is far more likely in close binaries.