Gas pile-up, gap overflow and Type 1.5 migration in circumbinary discs: general theory

Abstract

1 Introduction

Understanding the co-evolution of binaries and accretion discs is fundamental in several fields of astrophysics, including planet formation and migration (Goldreich & Tremaine 1980; Ward 1997), patterns in planetary rings (Goldreich & Tremaine 1982), stellar binaries, compact object and binaries involving supermassive black holes (SMBHs).

Despite the long history of the subject, there are no self-consistent analytical models for the co-evolution of binaries and accretion discs, incorporating the fundamental physical effects over the long time-scales on which the binary separation evolves. The standard α-model of radiatively efficient turbulent thin accretion discs (Shakura & Sunyaev 1973) relates the effective kinematic viscosity of the disc to the pressure ν ∝ αp. The viscous evolution of the disc, however, is often modelled without considering the pressure dependence of the viscosity (Lynden-Bell & Pringle 1974). Similarly, models of the gravitational interaction between the disc, which describe the launching of spiral density waves in the disc that remove angular momentum from the binary, also do not account for the tidal heating of the disc and the corresponding feedback on the torque cut-off phenomenon (Goldreich & Tremaine 1980).

The evolution of the circumbinary disc is sensitive to the above-mentioned assumptions, especially when the mass of the secondary is large, and can strongly perturb the disc. For a massive secondary, the tidal torque clears a gap in the disc, and the viscous radial inflow of the gas pushes the object inwards on the viscous time-scale (Type II migration). If the secondary mass, m_s, is larger than the local disc mass, , where Σ is the surface density, then the migration slows down, as the spiral density waves cannot remove angular momentum away from the binary at a rate on which the gas flows in. This leads to the pile-up of gas outside the secondary's orbit, in which the gas density increases by up to a factor B^−3/8, where B = m_d0/m_s < 1 and m_d0 is the unperturbed local disc mass (Syer & Clarke 1995). Once this steady-state level is reached, the viscous gas inflow velocity matches the inward migration of the object (secondary-dominated Type II migration).

In this paper, we focus on such systems, with m_s > m_d0, and point out that the Syer & Clarke (1995) steady-state level of gas pile-up cannot be reached for sufficiently large secondary masses m_s ≫ m_d0, where B ≪ 1. The enhanced viscous dissipation rate (D_ν ∝ B^−5/8) can increase the disc temperature such that it becomes radiation pressure dominated. The enhanced pressure makes the disc puff up (H ∝ B^−5/8), and reduces the relative gap size. Once the gas approaches within a distance less than the scaleheight from the secondary, the torque that the disc exerts on the binary has a cutoff (Goldreich & Tremaine 1980; Artymowicz 1993b, 1993a; Goodman & Rafikov 2001) which limits the migration rate of the secondary. Once the gas enters the Hill radius, it can furthermore flow across the secondary's orbit along horseshoe orbits or accrete on to the secondary. We derive an analytical quasi-steady-state model for the co-evolution of the disc and the orbital migration of the secondary, in which we combine a Shakura & Sunyaev (1973) disc with the theory of the binary--disc interaction by Goldreich & Tremaine (1980) self-consistently. In particular, we adopt the viscosity prescription of standard thin accretion discs proportional to pressure,¹ calculate the sound speed and vertical balance including both gas and radiation pressure (p_gas and p_rad), adopt the simple analytical approximation to the angular momentum exchange between the binary and the disc of Armitage & Natarajan (2002), consider the standard viscous and tidal heating of the disc (Lodato et al. 2009), and self-consistently account for the feedback on the pressure, viscosity, scaleheight and the torque cutoff near the secondary's orbit. We generalize the steady-state model of Hourigan & Ward (1984), Ward & Hourigan (1989) and Liu & Shapiro (2010) by self-consistently including variations in the viscosity and pressure caused by the pile-up. We derive azimuthally averaged steady-state analytical disc models which recover the Goodman & Tan (2004) solution for arbitrary β = p_gas/(p_gas + p_rad) in the limit that the secondary mass m_s approaches zero, but the disc structure is significantly modified by the secondary over multiple accretion time-scales for larger m_s.

The disc structure in this overflowing state with a pile-up is intermediate between the weakly perturbed case without a secondary and the case with a gap. Not surprisingly, the migration rate in such an intermediate state, which we label Type 1.5, is significantly different from the corresponding limiting cases of Type I migration and the secondary-dominated Type II migration. The transition between Type I and Type II migration as a function of the secondary mass was previously typically investigated by considering only the change in the surface density due to gap formation, but without investigating the feedback from the changes in viscosity and pressure (Hourigan & Ward 1984; Ward & Hourigan 1989; Korycansky & Papaloizou 1996; Ward 1997; Bate et al. 2003; Crida & Morbidelli 2007). However, simulations show that migration is sensitive to temperature variations and radiation pressure (D'Angelo, Henning & Kley 2003; Paardekooper & Mellema 2006, 2008; Kley & Crida 2008). We derive the Type 1.5 migration rate for the self-consistent radial profile including these effects when the pile-up is significant in an overflowing steady-state disc. As the secondary migrates inwards across the increasingly hotter inner regions of the disc, the gap opening conditions and the migration rate change even if one neglects the feedback on viscosity and temperature due to gas pile-up (Haiman et al. 2009; Kocsis, Yunes & Loeb 2011), but here we show that the changes are significant over a much wider range of masses and radii in the self-consistent model. We discuss migration and gap opening for SMBH binaries in more detail in Kocsis, Haiman & Loeb (2012, hereafter Paper II).

The remainder of the paper is organized as follows. In Section 2, we lay out the basic equations governing the hydrodynamical and thermal evolution of the disc, as well as the migration of the secondary. We solve the equations numerically in Section 3. We then derive an analytical solution in Section 4. We summarize the results, and discuss how they depend on the most important physical parameters, in Section 5. We offer our conclusions in Section 6. A more detailed discussion and the implications for SMBH binary systems are presented in Paper II.

We use geometrical units G = c = 1, and suppress factors of G/c² and G/c³ to convert between mass, length and time units. Our basic notation for the disc and secondary parameters is depicted in Fig. 1.

Figure 1

Gas pile-up and overflow in a circumbinary accretion disc with component masses M_• and m_s, binary separation r_s, and accretion rate . We distinguish five distinct radial zones: an inner and an outer far zone where the effects of the secondary are negligible, an interior and an exterior near zone where the tidal effects are significant, and an extended middle zone with a significant gas pile-up (see Section 4).

Open in new tabDownload slide

2 Thermo-Hydrodynamical Interaction Between A Disc and A Secondary

We examine the evolution of the secondary and an azimuthally and vertically averaged Shakura–Sunyaev disc (i.e. axisymmetric one-zone disc) in local thermal equilibrium. Here we review the basic equations. First, we write down the continuity and angular momentum transport equations including the viscous torque and the gravitational tidal torque of the secondary. The back-reaction of the tidal torque changes the angular momentum of the secondary. The viscous and tidal torques depend on the disc surface density, viscosity and pressure gradient (or scaleheight). We derive the vertically averaged disc structure assuming that: (i) the local viscous plus tidal heating equals the radiative cooling with photon diffusion limited vertical radiative flux (i.e. negligible radial heat transport); (ii) the viscosity is proportional to either gas+radiation pressure (α-disc) or just the gas pressure (β-disc) and (iii) gas plus radiation pressure supports the disc against the vertical gravity. This yields a closed set of non-linear partial differential equations for the disc and the location of the secondary in 1+1 dimensions. We present solutions in subsequent sections below.

2.1 Angular momentum transport

We denote the masses of the primary and the secondary objects by M_• and m_s, the surface density of the disc by Σ (assuming axisymmetry) and the radial bulk velocity of the disc by v_r, which is negative if gas accretes towards r = 0. The continuity and angular momentum equations for the disc are²

1

2

where the total torque T = −T_ν + T_d is due to viscosity and the gravity of the secondary, given by

3

4

Here Λ is the torque per unit mass in the disc, approximately given by

5

where

6

q ≡ m_s/M_•, H ≪ r is the scaleheight of the disc and f is a constant calibrated with simulations. This approximate formula for Λ, introduced by Armitage & Natarajan (2002), accounts for the net contribution of all Lindblad resonances as well as the torque cutoff within r_s ± H (Goldreich & Tremaine 1980; Ward 1997), and guarantees that the torque vanishes at r ≫ r_s. Here outside the torque cutoff in Goldreich & Tremaine (1980), in Lin & Papaloizou (1986), and f = 10⁻² calibrated to match the gap opening conditions in Armitage & Natarajan (2002).³ We adopt a conservative value f₋₂ ≡ f/10⁻² ∼ 1 in our numerical calculations, but keep the f₋₂ terms general in all of our analytical formulae. Note that practically equation (6) assumes that the tidal torque density ‘saturates’ instead of having a true cutoff near the secondary as long as the gas density is non-vanishing there (Artymowicz 1993b), which accounts for the effects of shocks near the secondary (Goodman & Rafikov 2001; Dong, Rafikov & Stone 2011b; Duffell & MacFadyen 2012).^4,5 However, this prescription might be inaccurate for a high-mass secondary forming a gap in the disc, where the tidal torques are due to spiral streams passing near the secondary on horseshoe orbits (MacFadyen & Milosavljević 2008; Baruteau, Ramirez-Ruiz & Masset 2012; Petrovich & Rafikov 2012; Roedig et al. 2012; Shi et al. 2012). We do not consider the torques inside the Hill radius, |r − r_s| < r_H ≡ (q/3)^1/3r_s, assuming that gas reaching this region flows in across the secondary's orbit. Outside this region, we use (4), assume that gravity is dominated by M_• and the orbital velocity is nearly Keplerian, , where .

After some algebra (Frank, King & Raine 2002), equations (1)–(2) simplify to

7

The total mass flux across a ring of radius r is defined as

8

Equation equation (7) along with the definition of the total torque T in equations (3)–(5) describes the evolution of the axisymmetric disc surface density and radial velocity as a function of radius and time.

The evolution of the secondary's orbital radius, r_s, is driven by the tidal torques of the gas and gravitational wave (GW) losses. The angular momentum of the secondary is so that

9

where −T_d is the recoil due to the torque exerted on the disc, equation (4), and with , the torque from the GWs is given by

10

Given ν(r, t) and H(r, t), equations (7) and (9) provide three equations for the three unknowns: Σ(r, t), v_r(r, t) and v_sr(t).

We examine steady-state solutions to these equations where and so that is a constant.⁶ Note that in general the disc need not be in steady state. However, in many cases the inflow rate of gas may be much faster than the radial migration speed of the secondary |v_sr| ≪ |v_r|. If this is satisfied in a wide range of radii up to the outer edge of the disc, then the secondary is effectively stationary in the azimuthally averaged picture, and the radial profile of the disc might be expected to relax to a steady state, independent of the initial condition of the disc. We propose that the secondary then migrates slowly through a sequence of quasi-steady-state configurations of the disc with a fixed . Then, equation (8) becomes

11

This is a first-order ordinary differential equation for T_ν(r), once is specified for a specific disc model.

2.2 Boundary conditions

We distinguish two types of inner boundary conditions corresponding to whether or not the perturbation is strong enough to lead to a truncated disc with a wide hollow circular cavity. Here, ‘wide’ means wider than the Hill radius (see below).

• (I)
  If Σ(r) ≠ 0 all the way to the innermost stable circular orbit r_ISCO of M_• (i.e. the disc does not have a cavity), we require a zero-torque boundary condition (Novikov & Thorne 1973; Penna et al. 2010; Tanaka 2011; Zhu et al. 2012),

  

  12

  Starting with this boundary condition, we obtain, among other properties of the steady-state disc, the gas velocity profile v_r(r). As stated above, if this inflow velocity is much faster than the migration velocity of the secondary over a large range of radii, then one might expect that the disc approaches this steady-state configuration, independent of the initial condition. In the opposite case, the steady-state assumption may be violated by the time-dependent migration of the secondary. As we will show, the steady-state solution with a fixed requires a large build-up of gas outside the secondary for the viscosity to overcome the tidal barrier of the secondary. We refer to these solutions, in which the disc is not truncated outside the secondary, as ‘overflowing’.
• (II)
  If the tidal torques dominate over the viscous torques near the secondary, gas is expelled from the region near the secondary and a wide gap forms. Assuming that the characteristic radius r_g, where the tidal torque is exerted on the disc near the edge of the gap, tracks the inward migration of the secondary with r_g = λr_s where 1 < λ ≲ 3 is a constant, we require that the gas velocity at that radius satisfies (Syer & Clarke 1995; Ivanov, Papaloizou & Polnarev 1999)

  

  13

  Note that λ is not specified by hand ab initio; it is found by assuming steady state in our solutions below. This condition can be understood intuitively, since the secondary cannot ‘run away’ and leave the outer disc behind (if it did, it would cease to be able to torque the disc and would have to slow down). Likewise, the gap edge cannot get closer to the secondary (if it did, gas would pile-up and the gap would eventually close). Although the disc is not in steady state near its boundary, we assume at r > r_g (see discussion in Section 4.1.4.).

Based on equations (8) and (9), equation (13) is equivalent to

14

Note that here and throughout the paper, by ‘gap’ we refer to situations where the gas density becomes effectively zero outside the secondary, such that the inflow of gas from the outside pushes the secondary inwards according to (13). In these cases, we assume that inflow across the orbit is insignificant, and in particular, we neglect torques from the gas interior to the orbit. In our calculations, a gap is effectively a hollow circular cavity in the disc, which is supported by the tidal torques of the secondary. However we emphasize that we do not rule out the presence of a local density decrement, resembling an annular gap, with a significant mass flux across the gap.

In practice, we attempt to find a solution with either of the above two boundary conditions, and then check whether the solution is self-consistent. By construction, only one of the two boundary conditions will lead to a self-consistent solution as confirmed below.

2.3 Physical conditions in the disc

Next we derive H(r) and ν(r) which appear in the tidal and viscous torques in equations (3)–(5).

2.3.1 Vertical balance

Let us first derive the scaleheight, H. If the vertical gravity is dominated by M_• (i.e. |r − r_s| > r_H), then in vertical hydrostatic equilibrium H = c_s/Ω where is the local midplane sound speed and p = p_gas + p_rad is the pressure due to the gas and radiation,⁷p_gas = ρkT_c/(μm_p), , where T_c is the central temperature, a = 4σ/c is the radiation constant, σ is the Stefan–Boltzmann constant, m_p is the proton mass and μ = 0.615 is the mean particle mass in units of m_p. Since ρ = Σ/(2H), so that c_s = 2p/(ΣΩ). The pressure can be expressed as p = p_rad/(1 − β), where β = p_gas/p. If photons are transported to the surface by diffusion, then the mean radiation flux is

15

Here T_s is the surface temperature, τ = κΣ/2 is the optical depth from the midplane to the surface, where κ = 0.35 cm² g⁻¹ is the opacity assumed to be dominated by electron-scattering. We do not investigate changes caused by free–free opacity at large radii for simplicity, and neglect deviations from blackbody radiation (see e.g. Tanaka & Menou 2010 for a more detailed model). Thus, , so that c_s = κc⁻¹F Ω⁻¹(1 − β)⁻¹, and we have

16

Note that equation (16) is valid in general for radiation flux limited, geometrically thin discs, independent of the source of dissipation and viscosity.⁸

2.3.2 Viscosity

In the standard Shakura–Sunyaev α- and β-disc models, the viscous stress tensor, t_ij = ρν∇_iv_j, satisfies , where b = 0 or 1, respectively, and α is a constant parameter (Shakura & Sunyaev 1973; Sakimoto & Coroniti 1981), implying that

17

In the second equality, we have substituted equation (16).

2.3.3 Local thermal equilibrium

We assume steady-state thermal equilibrium in which heat generated by viscosity and the dissipation of the spiral density wave escape the optically thick disc in the vertical direction by photon diffusion. The vertical radiation flux is F = D_ν + D_d. The viscous dissipation rate per disc face element is

18

We assume that the density waves generated by the tidal torque are dissipated locally in the disc and turned into heat, yielding the rate D_d. This is expected to be an adequate approximation based on analytical arguments (Goldreich & Tremaine 1980, equation (97) therein) and numerical studies (Dong et al. 2011b; Duffell & MacFadyen 2012; Rafikov & Petrovich 2012), especially in the regime where the disc is strongly perturbed.

Following Goodman & Rafikov (2001) and Lodato et al. (2009),⁹

19

The total vertical flux or total dissipation rate is

20

Using the above equations we derive Σ and T_c for a given D_ν and F at each radius (see Appendix A1).

2.3.4 Summary

Combining the previous expressions, we obtain

21

22

where

23

All other disc parameters can be derived from these relations. For example, the scaleheight H and the quantity νΣ that determine the torque (equation (11)) are given by equations (16) and (18). In particular, the limiting cases for H are

24

In the limit that the only source of heat is viscosity in a Keplerian disc, , we recover the solution of Goodman (2003) up to a constant of the order of unity.¹⁰

More generally, equations (16), (21) and (23), along with the definition of D_ν and F in equations (18) and (20), and the angular momentum flow equation (11), provide a closed set of equations for the stationary disc, valid throughout the gas and radiation-pressure dominated regions for α and β discs. The solution is self-consistent if for all r, the disc is thin (H < r), the radiation flux is sub-Eddington (), the radial accretion velocity is subsonic (), radial heat transport is negligible, the self-gravity of the disc is negligible and the disc is stable against fragmentation (Q = c_sΩ/(πGΣ) ≥ 1), the disc is optically thick (τ = κΣ/2 ≥ 1) and the boundary conditions are satisfied (implying in particular that v_sr ≪ v_r across a wide range of radii for overflowing solutions; see Section 2.2).¹¹ We verify that these conditions are indeed satisfied for the overflowing solutions given below.

3 Disc Structure - Numerical Solutions

First we generate numerical steady-state solutions for tidally and viscously heated discs assuming that the migration rate is much smaller than the radial accretion velocity in the disc. These numerical solutions are useful to verify the detailed analytical estimates presented in the following section.

We proceed along the following steps.

• (i)

  Obtain the ratio of gas to total pressure, β = β(r, D_ν, F), by inverting equation (23). A unique solution is guaranteed by the intermediate value theorem, since the left-hand side is a monotonic function of β, mapping 0 < β < 1 to all positive real numbers, while the right-hand side is positive and independent of β.

• (ii)

  Substitute the solution for β in equations (16) and (21) to obtain H(r, D_ν, F) and Σ(r, D_ν, F).

• (iii)

  Substitute β, H and Σ in the definition of F in equation (20) to get an equation between F and D_ν for fixed r and r_s. Invert this relation to find F(r, r_s, D_ν). Similar to step (i), one can show that the solution exists and is unique.

• (iv)

  Using equations (18) and (19), obtain the function .

• (v)

  Substitute into equation (11) to obtain an expression for a fixed . Solve this differential equation for T_ν(r, r_s).

• (vi)

  Substituting back into D_ν and F, equations (18) and (20) and the formulae of step (ii), to get Σ(r, r_s), T_c(r, r_s) and H(r, r_s).

The complexity is related to the non-linearities in steps (i), (iii) and (v). Nevertheless, the solution exists and is unique in steps (i) and (iii). However, step (v) is a boundary value problem of a non-linear first-order differential equation, which can have many solutions. We solve the differential equation numerically upstream from an initial value T_ν(r_ISCO) = 0. Without the secondary, the solution is simply , which leads to the Shakura–Sunyaev disc. If q ≪ 1, then the secondary creates a small dip in T_ν(r) in its neighbourhood, where the depth of the minimum increases with q. For larger q, T_ν(r) becomes very small positive approaching the secondary from downstream, and the surface density approaches zero. In this regime, tidal heating dominates over viscous heating, and H > |r − r_s|, implying that the pressure gradient shifts the torques out of resonance, and the torque is suppressed according to equation (5). Since the adopted torque model is valid only outside the secondary's Hill radius, we stop the calculation at r_s − r_H, and restart it at r_i = r_s + r_H assuming that¹²

25

This has a similar effect to smoothing the torque interior to the Hill radius as done previously in Lin & Papaloizou (1986), Syer & Clarke (1995) and Lodato et al. (2009).

The solution is approximately self-consistent if the migration rate is slower than the radial gas velocity outside the secondary. However, if this is not satisfied, a cavity opens and the disc becomes truncated. In this case, we seek a different solution in step (v), which satisfies the boundary condition in equation (13). This is possible by increasing r_i in equation (25) where T_ν(r_i) ≈ 0, until equation (37) is satisfied. Here r_i can be identified as the truncation radius at the inner edge of the disc. We distinguish the characteristic truncation or gap radius to reside at r_g where the tidal effect is exerted on the disc, more specifically the boundary where the tidal torque density becomes subdominant and use r_g in the boundary condition, equation (37).¹³ The surface density increases rapidly within r_i < r ≲ r_g, has a maximum and decreases thereafter. We assume that the disc is truncated interior to r_i if a gap forms with r_i > r_s + r_H.

4 DISC STRUCTURE - ANALYTICAL SOLUTIONS

Here we derive an analytical solution to the non-linear equations in Section 2. Such solutions can be derived asymptotically far from the secondary or near the secondary, where either the tidal torque or the viscous torque dominates, or where the angular momentum flux is negligible. We therefore distinguish the corresponding far, middle and near zones (see Fig. 1). The far zones are well inside and well outside the secondary, where the effects of the secondary are negligible. The middle zone is the region outside the secondary where the tidal effects (i.e. torque and heating) are locally negligible compared to the viscous effects, but where the gas pile-up is significant and the disc profile is modified. The near zones are just inside and just outside the secondary's orbit, where the tidal effects of the secondary dominate over the viscous effects. We restrict the near zone to outside the Hill radius, where the adopted tidal torque formula is valid. In addition to providing a basic understanding of the disc structure, the approximate analytical solutions allow us to infer the migration rate of the secondary.

To keep track of the approximations and notations introduced for the various zones below, we provide Table 1 for convenience. Note that the far/middle/near zones divide the disc into five radial slices, and the asymptotic behaviour further depends on whether the disc becomes truncated forming a wide gap (in the middle zone) and whether the torque saturated by the condition on the radial distance from the secondary is δr ≡ r − r_s < H (in the outer near-zone). Each row in the table corresponds to one of these disc regimes, discussed in a corresponding subsection below, and shows which terms are relevant in equation (11). The subdominant terms are marked with a ‘0’. The column with T_ν shows functions we introduced related to the viscous torque, and shows the scaling of the specific tidal torque in equation (5).

Table 1

Approximations and notations for the various radial zones in the disc used in Section 4

Open in new tabDownload slide

In the following we mostly focus on β-discs (i.e. b = 1) and examine both radiation and gas pressure-dominated discs, but it is straightforward to derive analogous formulae for α-discs in the same way. We also note that in the radiation-pressure-dominated regime, the viscosity of α-discs is larger by a factor of p_gas/(p_gas + p_rad) = β⁻¹. This would generally lead to stronger overflows for a smaller gas pile-up, and the cavity would close for a wider range of parameters than we find below for β-discs.

4.1 Far and middle zones

First we examine the region sufficiently far from the secondary, either inside or outside of its orbit, where

26

In this region, equation (11) can be integrated and substituted in (18)

27

28

where T_bc is an integration constant determined by the boundary condition near the secondary. For a fixed T_bc, equation (28) gives both F and D_ν, from which the surface density and central temperature follow from (21) and (22),

29

30

where

31

Thus, solving the disc structure in these zones amounts to finding the torque at the boundary, T_bc.

If T_ν(r_min ) = 0 then equation (11) shows that

32

In practice, r_min = r_ISCO for a disc without a cavity, and it is the inner edge of the disc if it has a cavity. Depending on which term dominates in equation (27), we distinguish the far zone () and the middle zone (). The far zone can be either well inside or far outside the secondary's orbit, but the middle zone is always outside. Well inside the secondary, the second term can be neglected in equation (32), and well outside of it, the second term dominates and the integration domain can be extended to ∞. In both cases, T_bc is independent of r.

Equations equations (27) and (32) show that in the region outside the secondary, T_bc represents a torque barrier due to the secondary's tidal effects. This parameter can also be used to obtain the migration rate of the secondary. Indeed, combining equations (9) and (32) give

33

4.1.1 Far zone - unperturbed disc

Without the secondary equation (12) implies that . Substituting into equations (27) and (28) gives D_ν and F. Plugging in equations (8), (16) and (21)–(22) leads to the standard Shakura & Sunyaev (1973) solution

34

35

36

37

38

Here and below, the subscript ‘0’ denotes quantities related to the unperturbed disc, α₋₁ = α/0.1, , is the Eddington accretion rate for 10 per cent radiative efficiency, q₋₃ = q/10⁻³, M₇ = M_•/10⁷ M_⊙, r_s2 = r_s/10²M_•, and we introduced

39

Without the secondary, in the radiation pressure-dominated regime (β ≪ 1) the scaleheight is approximately constant, and increases approximately linearly further out where gas pressure dominates (β ∼ 1).

The viscous torque, for future reference:

40

where r₂ = r/10²M_•. Sufficiently far from the secondary, the disc is independent of the secondary and follows equations (34)–(38) with ϕ ≈ 1. However, the disc structure depends on the rate at which gas is allowed to flow in through .

4.1.2 Middle zone

Now let us consider the opposite limit, , where the steady-state perturbation to the torque is significant. In terms of the dimensionless torque barrier,

41

the formulae describing the unperturbed disc, (34), get modified by replacing the boundary term with

42

The disc quantities change to

43

44

45

46

47

and the migration rate follows from equation (33)

48

where we have assumed T_GW ≪ T_bc. Here and below, the superscript ‘m’ labels the middle zone. Note that the dimensionless angular momentum flux k can be interpreted as a brightening factor in the middle zone relative to the unperturbed disc; k_s is representative of the maximum brightening, if k(r) is extrapolated to r_s. In practice, the maximum brightening is even larger than k_s in the near zone due to tidal heating (see Section 4.3).

The disc is modified within a radial range where the dimensionless angular momentum flux satisfies k > 1. This sets the outer boundary of the middle zone, where the disc transitions to the far zone. From equation (41),

49

Equations equations (43)–(47) represent a disc with negligible inflow of angular momentum but an inner boundary condition with a large viscous torque, corresponding to the torque barrier. Such solutions are often (somewhat misleadingly) referred to as a decretion disc (Pringle 1991; Lodato et al. 2009). To avoid confusion, we emphasize that there is accretion (i.e. inflow) in this region, too, with a fixed . However, the radial accretion velocity is greatly reduced, while the surface density, temperature and scaleheight are all greatly increased, relative to an accretion disc around a single compact object.

So far in this subsection, we have derived a solution for an arbitrary torque barrier or k, without specifying its value. In general, k is given by equation (32), which depends on the tidal torque in the near zone. Thus, to complete the derivation of the disc structure in the middle zone, we are first required to obtain the disc structure in the near zone (which we will do in Section 4.2). However, in the case of the steady-state cavity, the particular form of the boundary condition allows us to directly infer k, independently of the near zone, up to a factor λ of the order of unity, which we show next.

4.1.3 Middle zone - steady-state disc with a cavity

When the tidal torque is sufficiently strong to clear a gap so that the secondary and the nearby gas move with a similar velocity, T_ν ∼ T_bc can be substantial over a large range of radii. From equations (13) and (14), this requires

50

Here and below, the superscript ‘g’ refers to solutions with a gap, and r_g = λr_s is the outer radius of the gap. For this value of T_bc, equation (28) gives F and D_ν, and Σ follows from (21). However, since the right hand side (RHS) of equation (50) depends on Σ itself, this gives an algebraic equation for T_bc. The solution is

51

Note that this is independent of the tidal torque model (i.e. the Λ or f in equation (5)), since here the tidal torque is set by the boundary condition of the gap. This solution breaks down, and becomes tidal torque dependent, if the gap closes, which we discuss in Section 4.2.

The dimensionless angular momentum flux from equation (41) is

52

In particular, near the secondary is the ratio of secondary mass to the accumulated local gas mass.¹⁴ The only free parameter in this zone is λ, which we determine explicitly in Section 4.3.2..

In the range , k^mg ≫ 1 and equations (43)–(47) give

53

54

55

56

where Σ(r) ∼ 0 at r ≤ λr_s. Outside , k^mg ≈ 0, and the disc approaches the unperturbed solution given by equations (34)–(38). In the transition zone, between the middle and far zones, , one needs to use equations (43)–(47) with k = k^mg (equation (52)).

When a cavity is present, the migration speed of the secondary follows from equations (33) and (51):

57

This expression is consistent with the secondary-dominated Type II migration rate of Syer & Clarke (1995) who assumed λ = 1. Note that the migration speed is slower than the gas inflow velocity without the secondary, , in equation (37). This is referred to as disc-dominated Type II migration, which is appropriate if the secondary mass is smaller than the unperturbed local disc mass (or equivalently k^mg ≥ 1), but large enough to open a gap.

It is interesting to note that the structure of the middle zone does not depend explicitly on the tidal torque model, (in particular Λ or the f₂ parameter in equation (5)); the dependence is implicit and arises only by fixing the value of λ. Physically, while the tidal torques are negligible in this region, the effects of the tidal torques are still communicated to the region by setting an effective hydrodynamical boundary condition. We determine λ in Section 4.3.2. and find that, in fact, it only weakly depends on .

4.1.4 Consistency of steady state

A basic assumption of our model is that the radial structure of the disc is in a quasi-steady state as the secondary migrates slowly inwards. To check the consistency of these steady-state solutions, we must verify that the implicit time-dependence in the surface density profile through r_s(t) does not violate the continuity equation (1) significantly, so that

58

Integrating over radius, the relative error in the accretion rate

59

which is ∼20 per cent near the gap edge. The error in T_bc based on equation (51) is ∼13 per cent. However, the error in the accretion rate exceeds unity at large radii, outside 4.2λr_s. The steady-state assumption breaks down because as the secondary migrates inwards, the steady-state gas density near the edge of the gap continuously increases with time. If is fixed near the gap edge to be a constant fraction of the Eddington value, the true accretion rate at larger radii must be larger, to supply material for the increasing gas density. Conversely, if is fixed at large radii, then it becomes smaller approaching the gap edge. Such non-steady-state solutions have been derived by Pringle (1991) and Ivanov et al. (1999) by solving the non-linear diffusion equation (7)for a fixed outer boundary condition, assuming that the viscosity can be expressed as ν = kΣ^ar^b where k, a and b are constants. In particular, Ivanov et al. (1999) derived a non-steady, but self-similar solution. In that solution, the migration is slower, and the angular momentum flux is lower, compared to the Syer & Clarke (1995) steady-state solutions with a fixed for the same binary and disc parameters (Syer & Clarke 1995). The quasi-steady migration rate and brightening factors for a truncated disc with λ > 1 are intermediate between the Syer & Clarke (1995) and Ivanov et al. (1999) solutions.

The steady-state condition is typically not violated in the overflowing solution over a wide radial range. For global-steady state, a necessary condition is

60

where we have used equations (29), (48) and (49) and defined and . This sets a maximum limit for the dimensionless angular momentum flux k_s. For a β-disc, this implies

61

and k_smax is larger (i.e. less restrictive) for radiation pressure-dominated α-discs.

An important qualitative difference between the overflowing model presented here and the Syer & Clarke (1995) model for a truncated disc is that γ_Σs > 0 for the former as we show below. In contrast, in the overflowing case, the excess surface density and the dimensionless angular momentum flux k_s in the middle zone both gradually decrease during the inward migration of the secondary. Thus, the excess surface density diffuses radially outwards. If k_s < k_smax then the diffusion is sufficiently fast to reach a global quasi-steady state throughout the middle zone. If this is not satisfied, then the outer parts of the middle zone cannot respond as quickly as the object moves inwards and the structure of the disc in these regions will depend on its previous history. However, since the viscous time-scale is always much smaller than the migration time-scale in at least the inner parts of the middle zone, the local disc structure of the overflowing solution might approach an approximate steady state there with a constant even if k_s ≳k_smax. The migration rate of the secondary depends on the near zone of the disc, which is expected to remain insensitive to perturbations in the outer parts of the middle zone in an overflowing disc.¹⁵ We leave a detailed investigation of the time-dependent overflowing solutions to future work (Salem et al., in preparation).

4.2 Near zone

Now let us consider the regions near the secondary where the tidal torque and heating are important. We discuss steady-state solutions inside and outside of the secondary's orbit, in turn, without and then with a circular cavity. Deriving the physical properties of the disc in this region is useful to provide an estimate of the torque barrier, T_bc, at the interface between the near zone and the middle zone. As explained previously, the torque barrier sets the overall scale of the physical properties in the middle zone, as well as the migration rate of the secondary. We therefore first compute the value of the torque barrier for an overflowing disc, as well as for a disc with a wide gap. In the latter case, we then compare the value with the torque barrier in the middle zone derived above (equation (51)). By equating the two, we can estimate the gap size (i.e. λ), and obtain the conditions for gap opening and closing.

4.2.1 Inside the secondary orbit

Consider the region just downstream the secondary, outside the torque cut-off Δ = r_s − r > H in equation (6), assuming a steady-state overflow (i.e. no hollow circular cavity). Based on (11), the viscous torque decreases in the vicinity of the secondary and for a sufficiently large secondary mass, the angular momentum exchange is dominated by the tidal torque. In this regime,

62

implying that

63

for a Keplerian disc (here and below, the superscript ni refers to the solutions in the inner near zone). Equation equation (20) shows that

64

The asymptotic limit corresponds to Δ ≪ r. From equations (8) and (16)

65

66

The latter equation shows that the secondary makes the disc thinner downstream if the disc is radiation pressure dominated. Combining equations (63)–(64) with (21) gives . The viscous torque then follows from (18). To first beyond leading order, for b = 1,

67

exhibits a sharp cut-off near the secondary.

We can verify that the working assumptions hold in this region. Equation equation (66) shows that Δ > H holds for all Δ, since the unperturbed disc is thin, H₀ < r. Since which implies that is indeed satisfied for sufficiently small Δ. Coincidentally, the assumption in equation (62) is satisfied within a distance Δ_ni from the secondary, where

68

and

69

The disc parameters in the region r_s − Δⁿⁱ ≲ r ≲ r_s − r_H are

70

71

72

73

74

During the inward migration of the secondary, r_s decreases, and the surface density evolves in a self-similar way. The surface density, midplane temperature, surface brightness and scaleheight all decrease significantly near the secondary with large q₋₃. The radial flow velocity becomes very large in the close vicinity, and the flow may become advection dominated there. However, we do not extrapolate this solution inside the Hill radius of the secondary because the torque model is invalid there.

It is remarkable that for a fixed , the fractional perturbation to the surface brightness and the radiation-pressure-dominated scaleheight are universal in this region, independent of the binary and disc parameters. This property is general for an arbitrary torque or viscosity model in radiatively efficient steady-state discs. The surface density in this regime is also independent of the viscosity model but it is sensitive to the torque model: it is set to ensure that the tidal torque matches the angular momentum flow associated with . The original value of the surface density is suppressed by a factor proportional to q⁻²Δ. These solutions are valid only for discs with relatively large secondary masses, such that Δ_ni > r_H, but in which there is gas inflow across the secondary orbit.

4.2.2 Outside the secondary - unsaturated torque

Next consider the region just outside the secondary. Here we examine the case where the secondary is massive enough for the tidal torques to be important. After a significant amount of gas pile-up the viscous torque eventually counteracts the tidal torque and creates a stationary inflow. In this regime the tidal and viscous torques counteract one another and both greatly exceed the momentum flux in equation (11) such that

75

The tidal heating rate is much larger than the viscous heating rate of a disc without a satellite, making the disc much hotter and thicker. Let us assume D_d ≫ D_ν here, equation (20) implying that

76

In this subsection we examine the case where the torque is not saturated, H(r) < r − r_s, so that Δ = r − r_s in equation (6). This is most relevant for relatively small mass ratios (i.e. typically q ≲ 10⁻³; see equation (95) and Paper II), where the banking up of the stream in this region is modest. Substitute equation (21) for Σ with b = 1, and use equation (18),

77

where a_Σ is the constant coefficient in equation (21). Solve this for F and plug back into equation (76)

78

Integrating both sides between r_i, the inner edge of this region (see further below for a discussion of the value of r_i) and r,

79

We are most interested in the case where the tidal torques increase T_ν substantially in this region. If so, we approximate T_ν(r_i) = 0. Substituting equation (5) and rearranging give

80

where the superscript neu refers to the case of unsaturated torque in the external near zone, and

81

In the second line, the approximation is accurate to within 6 per cent for Δ_i ≡ r_i − r_s ≤ 0.3 r_s. Note that depends on radius only through ζ; it increases monotonically and approaches a constant value, which depends very sensitively on Δ_i/r_s. In practice, one might expect

82

for a disc without a cavity because the tidal torque model is valid only outside this region, and within this distance the gas may flow across the secondary orbit along radial streams or horseshoe orbits. In the following we keep Δ_i/r_H general. We incorporate a factor of (r_H/r_s)^−5/2 in the prefactor of equation (80) and introduce a renormalized ζ, as

83

Then equation (80) becomes

84

where

85

86

The outer edge of this region is where equation (75) is first violated, i.e. where the viscous torque density¹⁶ becomes comparable to the accretion term . We substitute from equation (78) utilizing (80) and get

87

We label the radial distance of this interface from the secondary as . While this equation of a single variable can be easily solved numerically for for any fixed r_s and Δ_i, we may derive an analytical approximate solution as follows. Assuming , ζ is close to its asymptotic maximum, which implies

88

To obtain the dependence of the physical parameters on radius, we first derive D_ν and F by substituting equation (84) into (18) and (77). Then Σ and T_c follow from equations (21) and (22). In the range ,

89

90

91

92

93

Here r = r_s + Δ, and these equations are formally correct to first beyond leading order in Δ/r_s, but we find them to be a good approximation typically within 15 per cent even for Δ/r_s ≳ 1. Interestingly, all of these physical parameters have a local extremum in this zone. We label the distance corresponding to the maximum local disc luminosity as . We find that is a slowly decreasing function of Δ_i, which varies between 1.55 and 1.4 for 0 < Δ_i ≲ r_s.

To match T_ν at a radius , the interface between this region and the middle zone, we must set . Based on equations (84) we may use

94

This equation is typically valid in the overflowing case (see equation (88)), as long as q is large enough that (strongly perturbed solution), but not too large so that H < r − r_s. We discuss the solutions if the latter condition is violated in Section 4.2.3. below. Comparing equations (40) and (86) shows that the minimum mass ratio to cause a significant gas buildup with unsaturated torques:

95

For smaller masses, the disc structure is not modified significantly outside the secondary, and one has to use equation (67) for the torque at the inner boundary in equation (79) with . We do not show these more general but more complicated expressions here.

We note that the results in this section are sensitive to Δ_i if different from r_H, which sets the distance at which the gas can flow in across the secondary's orbit without significant resistance. While Δ_i ∼ r_H is reasonable based on the horseshoe orbits in the restricted three-body problem, we keep it as a free parameter in the following.

4.2.3 Outside the secondary - saturated torque

Here we again assume that (75)–(76) hold, but now examine the case of much higher secondary masses, where the scaleheight is increased so much that H(r) > r − r_s and the tidal perturbation enters the torque cut-off regime. We make the simplifying assumption here that H(r) > r − r_s holds throughout the near zone so that the tidal torque in equation (5) does not alternate between saturated and unsaturated. It is straightforward to obtain more general solutions, but we find this exclusively saturated or unsaturated assumption to be an excellent approximation in most cases.

Due to the large torque barrier and tidal heating, the disc in this region is typically radiation pressure dominated, and we accordingly assume β ≈ 0 for the analytical solutions below.

Equations equations (16), (18) and (21) show that and for b = 1, where a_Σ is a constant. Substituting into equation (76),

96

This equation can be solved for F as a function of r and T_ν. Plugging back into equation (75) leads to a separable first-order differential equation for T_ν

97

where a₁ is a constant independent of r and r_s. Now use¹⁷ Ω/Ω_s ≈ (r/r_s)^−3/2 and integrate both sides assuming T_ν(r_i) ≈ 0 for some r_i ≳ r_s. Here r_i is the radius where the torque model breaks down, for which we adopt the Hill radius around the secondary r_i = [1 + (q/3)^1/3]r_s if a gap does not form. Thus,

98

where we introduced a dimensionless function

99

100

101

Here is the incomplete beta function, and the last two lines are simple approximations, typically accurate to within 15 per cent. We can now use equations (18) to get D_ν. For T_ν(r_i) ≈ 0, after substituting the value of a₁, equation (98) yields

102

103

Here the superscript nes refers to the case of saturated torque in the external near zone. Note that depends on radius only through ψ. Close to the inner boundary of this region r_i, it grows quickly with δr/r and saturates to a constant at δr/r ∼ 1. This can be understood, since this solution neglects angular momentum flow, the viscous torque is equal to the integrated tidal torque density, and the latter has a cutoff at δr/r ∼ 1. We can verify that is indeed satisfied in this region (cf. equation (40)) and so the first assumption, equation (75), and T_ν(r_i) ≈ 0 are well justified.

Can we use the asymptotic maximum of as an estimate of the torque at the outer boundary of this region to estimate in the middle zone of an overflowing disc? In many cases no, because the disc transitions to the middle zone much closer δr ≪ r_s, implying that the torque at the outer boundary of this region can be much less than its asymptotic maximum. The outer boundary of this region is where H = r − r_s. To figure out exactly where this happens, we proceed to determine the disc structure in this region.

Now equations (103) and (96) give D_ν and F; all other disc parameters then follow from equation (8), (16) and (21)–(22). The result within r_i ≤ r ≲ r_s + H^nes is

104

105

106

107

108

We can now confirm the consistency of the second assumption, equation (76), using equations (103) and (106). Indeed, near the secondary since scales with a higher power of ψ. Σ and T_c have a maximum, v_r decreases and becomes practically constant, while H slowly increases in this regime.

The outer boundary of this region, , is where¹⁸

109

We use as an approximation to the transition radius to the middle zone. After substituting equation (108), equation (109) is a non-linear algebraic equation for . While it is easy to solve it numerically for any choice of parameters, it is still useful to derive approximate analytical solutions. We find the following method yields results that are accurate within 20 per cent for a wide range of parameters. Use (100), expand H^nes to second order in , use (1 + ax) ≈ (1 + x)^a for small x and 1 − x^a ≈ −aln x for small a. This gives an approximate relation

110

where δr_i ≡ r_i − r_s. Equation equation (110) can be solved analytically using the Lambert W-function (Corless et al. 1996)¹⁹

111

as

112

where the two forms are equivalent, and we have introduced

113

Finally, we substitute in equation (102) and use equation (101),

114

Note that depends logarithmically weakly on the disc parameters (Corless et al. 1996), in practice for 10⁻⁴ ≲ a ≤ 1/e = 0.368. Here a ≤ 1/e is required for this solution to exist, implying that q and r_s have to be sufficiently large and small, respectively. In the opposite case, the torque is unsaturated (Section 4.2.2).

4.3 Transition between near and middle zones

4.3.1 The case with overflow

The value of T_ν at the outer edge of the near zone is to be matched with that in the middle zone, T_bc. If the tidal torque is unsaturated in the near zone, we approximate T_bc with the asymptotic maximum value, . Otherwise, if it is saturated, then we set . Matching the middle and near zones at assumes that T_ν does not grow substantially in the transition region between the saturated near zone and the middle zone, i.e. outside of but within a radius where the tidal effects are still non-negligible. We find this approximation to be better than 10 per cent. Thus, to match the torque at the outer boundary of the near zone and the inner boundary of the overflowing middle zone, we combine the saturated and unsaturated cases as

115

where and are given by equations (86) and (114). If this satisfies for δr_i = r_H (see equation (51)), then the satellite migration velocity is less than the gas bulk local inflow velocity, and the overflowing steady-state solution is self-consistent. In the opposite case the disc forms a gap with δr_i > r_H.

The dimensionless angular momentum flux in the middle zone (equation (41)) is

116

117

Gap overflow causes the torque level to decrease in the middle zone, which suppresses k. We discuss gap closing in Section 4.4.

The disc flux in the near zone is larger, due to tidal heating, than in the middle zone. To show this, we next compare the luminosity of the near zone to the middle zone explicitly. In the case of unsaturated torques in the near zone, we find that the local disc luminosity, (see discussion following equation (91)), peaks sharply near . We find that the integrated flux from within a ring of width 0.5 r_H is approximately

118

In the unsaturated case, r_H/r_s = (q/3)^1/3 ≪ 1, and so the net luminosity of the middle zone exceeds that of the near zone. For saturated torques, the maximum brightness corresponds to the outer boundary of the near zone, . We assume an effective radial width given by equation (110). From equations (16) and (109), assuming a radiation pressure dominated near zone, the dimensionless angular momentum flux of the near zone relative to the unperturbed disc is

119

Comparing equations (117) and (119), and recalling that typically , we conclude that the luminosity of the saturated near zone can exceed that of the middle zone by a factor between ∼ 3 and 10 for q ∼ 0.1.

Given the dimensionless angular momentum flux in the middle zone, the disc parameters and the migration rate are given by equations (43)–(48). We substitute equations (116)–(117) to obtain an explicit formula for the disc parameters. The migration rate in the overflowing disc with saturated and unsaturated torques is, respectively,

120

121

4.3.2 Disc with a cavity

Let us next turn to the case with a gap. We determine the radial distance to the outer edge of the gap, λr_s, here by requiring that the migration velocity matches the rescaled gas inflow velocity in the middle zone as stated in the boundary condition, equation (13). With λ in hand, equations (53)–(57) determine the disc parameters in the middle zone. The solution is different when the tidal torque is unsaturated near the inner edge and when it is saturated, which we discuss in turn below.

First, assume that the torque is unsaturated all the way outside of the gap (H ≤ r − r_s). We use equation (87) to obtain ζ(λr_s, r_s, r_i) at the interface between the near and the middle zone. We substitute in (92) to obtain the gas velocity at the interface rescaled by 1/λ:

122

To obtain the migration velocity, we identify in equation (48), and substitute equation (80), and eliminate ζ using equation (87). This gives

123

The boundary condition, equation (13), states that equations (122) and (123) must be equal. This provides a non-linear equation for λ. We solve this equation perturbatively. To first beyond leading order,

124

where

125

The ‘’ subscript is introduced to distinguish the case with unsaturated torques.

Next, consider the case of the torque cut-off. The formulae in Section 4.2.3 are not limited to the overflowing case, as long as the torque is saturated (δr ≡ r − r_s ≤ H). If a gap opens then δr_i marks the distance to the edge of the disc in equations (98) and (99), for which δr_i > r_H. Here, δr_i can be eliminated using the boundary condition equation (13) as follows. We identify in equation (13) where marks the edge of the near zone according to equation (109), so that λ = 1 + δ₁ where . Combine equations (107)–(109) to eliminate ψ(r₁) from the bulk gas velocity at r₁

126

Similarly, assuming in equation (48), eliminate from equations (102) and (109), we get

127

The boundary condition, equation (13), states that equations (126) and (127) are equal. After rearranging, we get

128

The last two terms can be omitted within 10 per cent accuracy for 0 < δ < 3. This gives the characteristic gap scale in the torque cut-off zone

129

provided that λ_s − 1 > r_H/r_s (i.e. r − r_s > r_s + r_H); otherwise no gap is possible.

Thus, the dimensionless angular momentum flux follows after substituting into equation (52). In the unsaturated and saturated cases,

130

131

The migration speed of the secondary in case of a gap with unsaturated and saturated tidal torques, respectively, is

132

133

These estimates depend on the somewhat arbitrary definition of λ that we have adopted in the two cases. In the unsaturated case, we have identified it with the outer edge of the transition region between the near and middle zones, where , while in the saturated case, we considered it to be the inner edge of the transition region, where H = r − r_s. While these conventions could be modified, they do not affect the overflowing solution. They do, however, influence the secondary orbital radius where the gap closes, which we discuss next.

4.4 Gap opening and closing

The previous sections define the disc uniquely, which constitute the solution to the basic equations of Section 2. By looking at the solution, we can identify cases where a cavity is kept empty in steady state or when the disc overflows.

The cavity refills if the inner edge of the disc outside the secondary, r_i, falls within the Hill radius of the secondary. The viscous torque in the near zone, either or , is a monotonically decreasing function of r_i. Thus, if the disc is truncated, then the viscous torque at any radius in the near zone is decreased relative to its value for an overflowing disc, r_i = r_s + r_H (equation (82)). This shows that the state of the disc is uniquely determined by the smallest torque barrier,²⁰T_bc, or equivalently, the smallest dimensionless angular momentum flux:

134

given by equations (116)–(117) and (130)–(131). We therefore must distinguish four different possible cases of migration and disc behaviour. First, corresponds to the standard case with a wide gap, with the secondary exhibiting Type II migration. If , then the gap edge is located within a scaleheight in the near zone so that the torque cutoff limits the tidal torques, but they can nevertheless support a gap against viscosity as the secondary migrates inward. However, if , then the saturated tidal torque becomes smaller than the viscous torque all the way to the Hill radius, and the cavity refills. Finally, if , then the disc reaches the Hill radius and overflows already while the torques in the near zone are still unsaturated.

The migration rate of the secondary is proportional to k_s, ((48)), i.e.

135

given by equations (120)–(121) and (132)–(133). Note that and are given by practically the same formula, up to an order-of-unity factor of λ. This corresponds to the case of secondary-dominated Type II migration (Syer & Clarke 1995). The migration rate in the overflowing case is or for unsaturated or saturated torques. Here the disc is still strongly perturbed, but gas inflow across the orbit limits the efficiency of migration. The disc structure in this new regime is intermediate between a disc with an empty gap (normally associated with Type II migration) and a weakly perturbed disc (Type I migration). Although the migration speed in this regime is slower than either in standard Type II or Type I migration, we refer to this regime as ‘Type 1.5’.

The gap can also close due to three-dimensional overflow for large secondary masses if the tidal heating is substantial to make the disc puff up. Equation equation (56) shows that this happens (H ≳ r) in a radiation pressure-dominated disc if

136

or equivalently if

137

In this case, we do not derive the geometrically thick overflowing disc or the migration rate.

We discuss the gap opening and closing conditions in more detail in Paper II, where we contrast them explicitly with the standard expressions widely used in the literature, and also compare Type 1.5 migration to the previously known Type I and II cases.

5 RESULTS AND DISCUSSION

We have derived analytical solutions to the disc model in different radial regions, where either the tidal, the viscous torques, or the angular momentum flux is negligible relative to the other two terms in equation (11). In particular, we have identified the far zones, either well inside or outside the secondary's orbit, where the effects of the secondary are negligible, the exterior middle zone, where the disc structure is greatly modified but where the tidal torque and heating are locally negligible, and the near zones just inside or outside the secondary, where the tidal effects dominate.²¹ We distinguish two cases in the middle zone, depending on whether the disc has a gap (i.e. the disc is truncated well outside the Hill radius) or if the disc is overflowing across the secondary's orbit. We also distinguish two cases in the near zone outside the secondary's orbit, depending on whether the tidal torque from the binary is saturated (i.e. whether the location of the cavity edge falls within a scaleheight H from the secondary). Furthermore, we have investigated asymptotic results for gas and radiation pressure-dominated cases.

Distinguishing all of the above cases allowed us to adopt separate approximations, each valid in the corresponding regime, and to derive analytical results to the perturbed accretion disc interacting with the secondary. We further used this to estimate the migration speed for the secondary.

In this section, we collect all of the resulting analytical solutions for the most important disc parameters in the various zones, and present them in a form suitable for easy use. We refer the reader to Paper II for physical interpretations, and discussions on possible implications of our results for real binary systems.

5.1 Disc model

Our results in this paper apply to geometrically thin, optically thick accretion discs, and describe vertically and azimuthally averaged properties. All physical parameters can be written as

138

where X denotes any of {Σ, T_c, H, v_r, F, T_ν}; r₂ and r_s2 denote the radius from the primary and the orbital radius of the secondary in units 100 M_•, respectively; and denotes an extra function of the parameters . Note that are the usual parameters of a standard solitary accretion disc; q and f represent the mass ratio and the normalization of the azimuthally averaged tidal torque (see equation (5)). The −N index denotes normalizing with 10^−N, e.g. q₋₃ = q/10⁻³. The C prefactor, c_i exponents and the Φ function in equation (138) are given in Table 2 in the different zones and cases for β-discs (i.e. ν ∝ p_gas). For β-discs, many of the physical parameters, including the surface density, Σ, and the central temperature, T_c, are independent of whether gas or radiation pressure dominates. This, however, is not true for the scaleheight, where we quote results in both regimes, labelled with a ‘gas’ or ‘rad’ subscript.

Table 2

Pre-factors and exponents in the analytical disc model in different zones, C and c_i in equation (138). The third column is in cgs units except where marked by * (where it is in units of GM_•/c²). Columns 4–10 are exponents; the last column is the extra multiplicative function (see text). The last two block of parameters show the migration rate of the secondary and other useful parameters

Open in new tabDownload slide

The last column of Table 2 shows in equation (138). Here we introduced the following notation:

139

140

141

142

143

144

145

146

147

Here B(x; a, b) is the incomplete Beta function, and is the Lambert W-function, the branch defined on the negative real axis, evaluated at −a given in the last row of Table 2. Typically, . The approximations shown for ψ and are better than 20 per cent. Here, r_ISCO is the innermost stable circular orbit near the SMBH [i.e. 6 M_• (1 M_•) for a non-spinning (maximally spinning) SMBH)], r_H = (q/3)^1/3r_s is the Hill radius, r_i marks the radius at which the viscous torque becomes very small outside the secondary's orbit for which we use

148

and λ is the dimensionless gap or truncation radius scale

149

where²²

150

is the dimensionless angular momentum flux or brightening factor.

We collect the formulae for the transition radii separating different radial zones and physical regimes in Table 3. We label , where marks the radius of the interface between zone a and b. For example, is the transition between the middle zone and the torque-saturated exterior near zone if there is a wide gap, which also sets the truncation radius scale λ_s in equation (140).

Table 3

Transition radii between different zones relative to the secondary orbital radius, r_s. Here is the radius at the interface between zone a and b, , and . Different columns show the constant prefactor and exponents in equation (138) as in Table 2

Open in new tabDownload slide

The state of the disc and the migrate speed of the binary are directly set by the dimensionless angular momentum flux k_s. If k_s =[], then the disc is in the [un]saturated overflowing state, whereas if k_s =[] then it is in the [un]saturated state with a wide gap. These four possibilities are indicated in Table 2. The appropriate choice of k_s also determines which case in Table 2 is to be used in the middle zone for the other parameters (‘mou’, ‘mos’ or ‘mg’), and in the near exterior zone (‘neu’ or ‘nes’).

Thus, the ‘phase space’ of solutions consists of four regions.²³ The transition between two different solutions, with or without wide gaps, corresponds to the parameters for which or . Our hypothesis is that this must represent a physical transition in the disc+binary system, as the secondary migrates inwards from large radius. Initially, a central cavity is created, and the outer edge of the cavity lies far from the secondary's orbit. However, as the secondary migrates inwards, the distance between the cavity edge and the secondary shrinks (at least when measured in units of the Hill radius of the secondary). This may happen both because the viscosity increases as the pressure grows during pile-up, and also because the tidal torque decreases with increasing scaleheight due to the torque cutoff. The cavity finally closes once the cavity wall nudges inside the Hill radius. The dimensionless angular momentum flux or brightening factor, k_s, is largest when this transition occurs.

These gap opening/closing conditions are quite different from those in the literature (e.g. Ward 1986), which state that the gap closes if the disc extends into the region closer than either the local scaleheight or the Hill radius from the secondary. Note that our gap closing condition combines statements on the scaleheight and the Hill radius, but since the scaleheight and viscosity vary significantly near the secondary in the strongly perturbed case with a large pile-up, they depend on the actual perturbed profiles rather than the averaged quantities describing accretion discs around a solitary object. We discuss gap opening/closing, and its physical implications, in more detail in Paper II.

We have verified that the analytical approximate solutions to the disc model match the numerical solutions typically to within tens of per cent for a wide range of disc and binary parameters when the disc is strongly perturbed. Fig. 2 shows two examples when the disc is in the overflowing state with unsaturated (left-hand panel) and saturated (right-hand panel) tidal torques near the secondary. Different regions are indicated with the abbreviations used in Table 2. For further examples and other physical quantities, see Paper II.

Figure 2

Comparison of the numerical solution (thick yellow solid lines) with the asymptotic analytical approximations in the various zones (dotted, dashed and dash–dotted lines, labelled as in Table 2) for the surface density of the disc around a 10⁵ M_⊙ primary. The mass ratio and binary separation are (q, r_s) = (0.01, 500 M) and (0.1, 1000 M) on the left- and right-hand panels, respectively. In both cases the disc is in the overflowing steady state. The tidal torque is unsaturated in the near zone outside the Hill radius on the left-hand panel, but it is saturated on the right-hand panel. The disc structure is significantly modified in both cases in an extended region around the secondary.

Open in new tabDownload slide

5.2 Brightening factor

As stated in the previous sections, the k_s parameter sets the angular momentum flux and the brightness of the disc in the middle zone relative to the unperturbed value, and determines the state of the disc.

Remarkably, the disc parameters can differ dramatically from the unperturbed values not only in the near zone, where the tidal effects dominate, but also in the middle zone, where tidal effects are already negligible. The tidal effects of the secondary are short range, but the corresponding effect is communicated to distant regions by setting an effective boundary condition. The size of this region can be times larger than the secondary orbital radius (see Table 3). This long-range behaviour is confirmed in our full numerical solutions and is also present in Lodato et al. (2009) and Chang et al. (2010) for a circumbinary disc with a cavity. The same k_s parameter also sets the increase in the scaleheight, as well as the migration speed of the secondary. We have derived the analytical formulae for k_s in the four relevant regimes summarized in the previous section. Due to its extended radial size, the integrated luminosity of the middle zone is often larger than that of the near zone in the torque-unsaturated case, overflowing state. However, the near zone may be brighter than the middle zone, when the system is in the torque-saturated, overflowing state (see Paper II for further discussions).

5.3 Migration rate

Finally, we have derived the migration speed of the secondary; approximate analytical formulae for our results are listed in Table 2. We find that in the regimes when a cavity forms, the migration in our solution is slower than the standard steady-state secondary-dominated Type II migration rate (Syer & Clarke 1995) by a factor λ^11/16 ∼ 2. The Type II rate may be further reduced in non-steady state models if the accretion rate or the gas mass is limited (Ivanov et al. 1999; Lodato et al. 2009). However, in the new ‘Type 1.5’ regime we identify, with a partial pile-up and overflow, the migration is even slower.

More generally, the migration speed for a strongly perturbed disc follows the minimum of the three possible solutions:

151

In general, Type 1.5 migration is more rapid at larger r_s, in contrast with the Type II rate, which is nearly constant, and the GW inspiral rate, which strongly decreases with r_s. Therefore, as a real system evolves, it will first transition from an initial Type II migration at large radii to Type 1.5 migration at smaller radii, before finally being driven by GWs at still smaller separations.

Regarding the dependence on the primary mass, Type II is nearly constant, while the Type 1.5 speed increases with M_•. This implies that at any fixed orbital separation, Type 1.5 migration is relevant for lower masses (roughly those in the range expected to be detectable by a space-based gravitational wave mission such as LISA, and Type II is relevant for higher-mass binaries (in the sensitivity range of Pulsar Timing Arrays; see Paper II). Type 1.5 migration may also be important for stellar mass binaries in proto-stellar discs or Jupiter-mass planets around M-type dwarf stars of mass (0.1–1) M_⊙ (see Johnson et al. 2012, for a recent discovery of such a system). Finally, we note that Type 1.5 (Type II) migration operates for larger (smaller) accretion rates, if fixing all other parameters.

All of the above conclusions are based on the analytical solutions we obtained; however, we have verified, by numerically solving the equations presented in Section 2, that our solutions are accurate to within tens of per cent for a broad range of parameters. In particular, in Fig. 3 we show the analytical approximation of the migration rates for the case of r_s = 1000 M_• for M_• = 10⁵ M_⊙, together with the rates obtained from a numerical solution. In this case, a disc is in the overflowing state for all values of q ≳ 10⁻³ shown, and the migration rate is significantly different from both Type I and II. As the figure shows, the analytical Type 1.5 formulae give a good approximation over a wide range of q for the migration speed.

Figure 3

The migration speed of the secondary for different mass ratios for M_• = 10⁵ M_⊙ at r_s = 1000M_• (or ). The steady-state numerical solutions (blue) are well represented by our analytical formula for Type 1.5 migration (magenta), but not by those for Type I or II (black dotted). The dashed and solid magenta lines show the unsaturated and saturated Type 1.5 cases, respectively. The discrepancy between the numerical and analytical solutions becomes more significant at q < 0.002 due to the fact that the disc is not strongly perturbed there.

Open in new tabDownload slide

5.4 Caveats

Our findings are subject to many possible caveats.

• (i)

  We assumed a radiatively efficient disc model in which the effective viscosity is proportional to the gas pressure in the disc with a constant α coefficient even in the radiation pressure-dominated regime. Future studies should investigate alternative models in which the viscosity is proportional to the total gas+radiation pressure (Shakura & Sunyaev 1973), or where the viscosity is generated by magneto-rotational instability (MRI; see Turner et al. 2003; Shi et al. 2012; Giacomazzo et al. 2012; Noble et al. 2012 for simulations of circumbinary discs leading to an ‘antigap’).

• (ii)

  We assumed steady-state models where the accretion rate is constant over radius. This is expected to be valid as long as the gas inflow is much faster than the migration rate of the secondary over a large range of radii, if the total gas supply is not limited, and if the accretion rate is set at the inner or outer boundary (i.e. for the new Type 1.5 migration regime we focus on here). However, this assumption may be violated for tidally truncated circumbinary discs (Type II migration; Ivanov et al. 1999) or for models where the gas supply rate is limited at the outer boundary (Lodato et al. 2009; Rafikov 2012).

• (iii)

  We assumed unequal-mass binaries, averaged over the azimuthal angle, assumed that the density waves generated by secondary are dissipated locally, and that the radial tidal torque profile follows the formula given by Armitage & Natarajan (2002). This assumes that the tidal torque saturates near the secondary at a radial distance closer than the scaleheight (). We also assumed that the gas entering a distance comparable to the Hill radius can flow freely across the secondary's orbit. However, accretion on to the secondary may affect the Type 1.5 migration rate. Farther away from the secondary, the assumed torque density has a steep cutoff (); extrapolating beyond r > 2r_s might be inaccurate. The |r − r_s|⁻⁴ scaling may also be inaccurate in the local non-linearly perturbed regime especially for comparable mass-ratio binaries (MacFadyen & Milosavljević 2008; Petrovich & Rafikov 2012; Roedig et al. 2012). These issues should be investigated using simulations, which could also address comparable mass binaries where the disc may be significantly non-axisymmetric (e.g. MacFadyen & Milosavljević 2008; Cuadra et al. 2009) and where the accretion of the secondary is non-negligible (Lubow et al. 1999).

• (iv)

  We neglected non-axisymmetric inflow into the secondary's orbit or on to the secondary if a cavity is formed. Inflow across the gap or accretion of the secondary reduces the amount of pile up outside the secondary, reduces the Type II migration rate and could affect the gap closing transition between the continuously overflowing solutions and the cases with a gap. We have also neglected the corotation torques in the overflowing case.

• (v)

  We found that the enhanced pressure dominates over the increase in surface density outside the secondary's orbit, which makes the overflowing disc stable against gravitational fragmentation (see Paper II). However, the steep pressure gradient in the near zone around a massive secondary may lead to global non-axisymmetric dynamical instabilities (Papaloizou & Pringle 1985; Goldreich, Goodman & Narayan 1986). The corresponding enhancement of the effective viscosity and angular momentum transport in the disc might reduce the pile-up outside the secondary's orbit and further reduce the Type 1.5 migration rate. A detailed stability analysis and an investigation of its implications go beyond the scope of this paper.

• (vi)

  We restricted our attention to circular binaries. However, binary eccentricity may be excited in cases with a cavity, when the masses of the two compact objects are comparable and the gap edge itself becomes significantly non-axisymmetric (Artymowicz 1992; Armitage & Natarajan 2005; Roedig et al. 2011, 2012). We note that such non-axisymmetries excited in discs with a gap diminish once the mass ratio is q ≲ 0.1 (D'Orazio, Haiman & MacFadyen 2012). Nevertheless, it remains to be seen if the binary develops significant eccentricities in the unequal-mass, overflowing state, with a significant pile-up.

• (viii)

  We assumed that the binary is sufficiently widely separated that gravitational wave emission is negligible. For a complete picture, future studies should investigate the gravitational wave driven regime (Chang et al. 2010; Tanaka & Menou 2010; Farris, Liu & Shapiro 2011; Kocsis et al. 2011; Yunes et al. 2011; Bode et al. 2012; Baruteau et al. 2012; Giacomazzo et al. 2012; Noble et al. 2012; Tanaka, Menou & Haiman 2012).

6 Conclusions

In summary, we have presented new analytical solutions to disc properties and migration rates, obtained from self-consistent solutions of a coupled binary--disc system. The evolution equations are solved analytically in the strongly perturbed limit, including the angular momentum exchange between the disc and the binary and the modifications to the density, scaleheight and viscosity self-consistently, including viscous and tidal heating, diffusion limited cooling, radiation pressure and the orbital decay of the binary.

In addition to recovering solutions with a central cavity, similar to previous ‘Type II migration’ scenarios, we have identified a distinct new regime, applicable at smaller separations and masses, larger accretion rates and mass ratios in the range 10⁻³ ≲ q ≲ 0.1. For these systems, gas piles up outside the binary's orbit, but rather than creating a cavity, it continuously overflows as in a porous dam. The disc properties are intermediate between those in an unperturbed disc and a disc with a wide gap. The migration rate of the secondary in this ‘Type 1.5’ regime is typically slower than both Type I and Type II rates.

In this paper, we have presented simple analytical formulae that comprehensively describe binary systems with different parameters in various stages of evolution. The analytical results provide simple scaling relations, which may be useful to scale and interpret the results of numerical simulations to different disc or binary parameters. It allows us to map out the effects of varying α, , or the binary parameters, over a wide range from planetary discs to active galactic nuclei around SMBH binaries.

We discuss the applicability of the new Type 1.5 regime and its physical implications for specific systems, as well as possible observable signatures in a companion paper (Paper II).

Acknowledgments

We thank Re'em Sari, Taka Tanaka, Alberto Sesana and Roman Rafikov for useful discussions. BK acknowledges support from NASA through Einstein Postdoctoral Fellowship Award Number PF9-00063 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. This work was supported in part by NSF grant AST-0907890 and NASA grants NNX08AL43G and NNA09DB30A (to AL) and NASA grant NNX11AE05G (to ZH).

References

Appendix A: Viscously and Tidally Heated Discs

Here we provide the details of the algebraic manipulations that give the disc model for fixed F and D_ν for either α or β-discs (i.e. b = 0 or 1).

For fixed β, we can reduce the problem to three equations and three unknowns Σ, T, ν. From equations (15), (17) and (18),

A1

A2

Solve this for Σ and T,

A3

A4

Finally β is given by

A5

where Σ, T_c and H are to be substituted from equations (154), (155) and (16). From this

A6

Finally equations (154) and (155) can be simplified by eliminating 1 − β using equation (157), which leads to (21)–equations (22).

Author notes