text:
# Saturation of the corotation resonance in a gaseous disk

## 1 Introduction

The forcing of a gaseous disk by a gravitational perturber at a resonance can result in a strong response and an interchange of energy and angular momentum between the perturber and the disk. A corotation resonance arises where the pattern speed of the forcing matches the local angular speed of disk material. This occurs in the context of galaxies where a stellar bar or spiral arms force motions in the interstellar medium (Binney & Tremaine 1987). Corotation resonances also arise when satellites orbit within a disk, as occurs with young planets in protoplanetary disks, or moons in planetary ring systems (Goldreich & Tremaine 1979, 1980; hereafter GT79 and GT80 respectively).

For planets, two types of corotation resonance need to be distinguished. A coorbital corotation resonance occurs where the orbital period of the planet matches the orbital period of disk material. For a planet with a circular orbit, this is the only type of corotation resonance that can arise. The second type of corotation resonance is non-coorbital and occurs if the planet executes an eccentric orbit. The forcing due to a planet with an eccentric orbit can be decomposed into a series of rotating components of various strengths and pattern speeds (GT80). Bar or spiral galaxies, on the other hand, give rise to non-coorbital resonances.

The analyses of the coorbital and non-coorbital resonances are different. In the coorbital case, which has recently been considered by Masset (2001, 2002) and Balmforth & Korycansky (2001), the forcing is stronger and involves multiple Fourier components. There is more likely to be a reduction in density caused by the tendency of the planet to open a gap in the disk. In this paper, however, we are concerned with non-coorbital corotation resonances, which are critical in determining the eccentricity evolution of planets resulting from planet–disk interactions (GT80), and are therefore likely to be important in explaining the eccentricities of many of the observed extrasolar planets. For a disk in which a gap is opened by a planet having a small orbital eccentricity, corotation resonances tend to damp the planet’s eccentricity, while Lindblad resonances cause it to grow. The effects of each type of resonance are strong, but are nearly equal in magnitude, to the extent that they nearly cancel. The balance is slightly in favor of eccentricity damping, if the corotation resonances operate at maximal efficiency (i.e., are unsaturated). The final outcome of eccentricity evolution depends on the details.

The disturbances of the disk caused by a corotation resonance remain localized to the corotation region (GT79). They are unable to propagate away from the resonance, as occurs in the case of Lindblad resonances. Instead, they act back locally on the disk and may change its density and vorticity in such a manner as to reduce (saturate) the corotation torque (e.g., Ward 1991). However, the effective turbulent viscosity of the disk could lessen the saturation by limiting the back-reaction on the disk (Ward 1992). Full saturation occurs when the torque is reduced to zero, but even a small degree of saturation (5%) could change the sense of eccentricity evolution in the GT80 model from decay to growth (Goldreich & Sari 2002).

In this paper, we develop a detailed model for the saturation of a non-coorbital corotation resonance in a viscous accretion disk. We determine how the strength of the corotation torque varies with the strength of the tidal forcing and the effective viscosity of the disk.

## 2 Physical picture

The interaction between the perturber and the disk in the corotation region can be considered to involve two distinct aspects. In one aspect, the tidal potential excites an evanescent disturbance in the disk, i.e., a disturbance that decays exponentially with distance from the resonance over the corotation region. Associated with this disturbance are evanescent angular momentum fluxes that are oppositely directed on the two sides of the resonance. In the linear, inviscid theory of GT79, there is a jump in the flux at the resonance that corresponds to a torque being injected into the disk exactly at the corotation radius. The torque depends on the radial derivative of the ratio of the surface density to the vorticity at corotation. In the second aspect, the disk responds to the torque by changing the distributions of surface density and angular velocity over the corotation region (Lubow 1990). These changes are such as to cause a reduction in the torque, and, ultimately, it is expected that the torque will be eliminated. This picture is time-dependent because the torque injection continues until full saturation occurs.

In the case of an accretion disk, full saturation may not occur owing to the effects of turbulent viscosity in limiting the back-reaction. We therefore aim to find steady-state solutions in which the torque on the disk is nonzero. We also allow simultaneously for both aspects described above, by formulating a nonlinear problem in which the feedback is included self-consistently. In the process of the analysis, we will determine the size of the region over which the torque is injected into the disk.

## 3 Analysis of the corotation region

### 3.1 Basic equations

We consider a three-dimensional gaseous disk having a barotropic equation of state, so that the pressure $`p`$ depends only on the density $`\rho `$. We assume that the forcing potential and the properties of the unperturbed disk vary smoothly in radius across the corotation region. These assumptions may need to be reconsidered if the resonance is very close to the planet (i.e., at a distance comparable to the vertical scale-height $`H`$ of the disk), or close to a sharp disk edge (where the density varies on a radial scale comparable to $`H`$). However, the torque cut-off (GT79) limits the effectiveness of resonances that are closer than approximately $`H`$, and a gap in the disk about the planet’s orbit may eliminate such close resonances in any case. We also consider the planet to execute a fixed orbit, and ignore any effects of migration.

The disk is governed by the equation of mass conservation,

$$(_t+𝒖)\rho =\rho 𝒖,$$
(1)
and the equation of motion,

$$(_t+𝒖)𝒖=(\mathrm{\Phi }+h)+\frac{1}{\rho }𝐓.$$
(2)
Here $`𝒖`$ is the velocity, $`\mathrm{\Phi }`$ the gravitational potential, $`h=\rho ^1𝑑p`$ the enthalpy, and $`𝐓`$ the stress tensor, which represents a viscous or turbulent stress. For a Navier–Stokes model we have

$$𝐓=\mu \left[𝒖+(𝒖)^\mathrm{T}\right]+(\mu _\mathrm{b}\frac{2}{3}\mu )(𝒖)\mathrm{𝟏},$$
(3)
where $`\mu `$ is the shear viscosity and $`\mu _\mathrm{b}`$ the bulk viscosity. We assume that, as is conventional in accretion disk theory, the vertically integrated shear viscosity $`\nu \mathrm{\Sigma }`$ may be expressed as a function of the radius $`r`$ and the surface density $`\mathrm{\Sigma }`$. The logarithmic derivative $`D_\mu =\mathrm{ln}(\nu \mathrm{\Sigma })/\mathrm{ln}\mathrm{\Sigma }`$ will feature in the analysis below.

The unperturbed disk is a steady, axisymmetric solution $`\{\rho ^{(\mathrm{u})},𝒖^{(\mathrm{u})},h^{(\mathrm{u})},𝐓^{(\mathrm{u})}\}`$ of these equations in the presence of the steady, axisymmetric potential $`\mathrm{\Phi }^{(\mathrm{u})}`$ of the central mass. Let the potential be perturbed by the addition of a nonaxisymmetric external potential $`\mathrm{\Phi }^{}`$ that rotates rigidly with angular pattern speed $`\mathrm{\Omega }_\mathrm{p}`$. The solution in the presence of the perturbed potential $`\mathrm{\Phi }^{(\mathrm{p})}=\mathrm{\Phi }^{(\mathrm{u})}+\mathrm{\Phi }^{}`$ is $`\{\rho ^{(\mathrm{p})},𝒖^{(\mathrm{p})},h^{(\mathrm{p})},𝐓^{(\mathrm{p})}\}`$, and may be assumed to be steady in a frame of reference that rotates at $`\mathrm{\Omega }_\mathrm{p}`$.

We write the nonlinear equations for the Eulerian perturbations $`\rho ^{}=\rho ^{(\mathrm{p})}\rho ^{(\mathrm{u})}`$, etc., in cylindrical polar coordinates $`(r,\varphi ,z)`$ in the rotating frame, using $`(u,v,w)`$ to denote the cylindrical velocity components, and omitting the superscript (u):

$$D\rho ^{}+(u^{}_r+w^{}_z)\rho =\rho ^{}\left[\frac{1}{r}_r(ru)+_zw\right](\rho +\rho ^{})\left[\frac{1}{r}_r(ru^{})+\frac{1}{r}_\varphi v^{}+_zw^{}\right],$$
(4)
$$Du^{}+(u^{}_r+w^{}_z)u2\mathrm{\Omega }v^{}\frac{v_{}^{}{}_{}{}^{2}}{r}=_r(\mathrm{\Phi }^{}+h^{})+\mathrm{},$$
(5)
$$\frac{1}{r}D(rv^{})+2Bu^{}=\frac{1}{r}_\varphi (\mathrm{\Phi }^{}+h^{})+\mathrm{},$$
(6)
$$Dw^{}+(u^{}_r+w^{}_z)w=_z(\mathrm{\Phi }^{}+h^{})+\mathrm{}.$$
(7)
Here

$$D=(u+u^{})_r+\left(\mathrm{\Omega }\mathrm{\Omega }_\mathrm{p}+\frac{v^{}}{r}\right)_\varphi +(w+w^{})_z,$$
(8)
and the dots denote viscous terms that we do not write out owing to their complexity. $`\mathrm{\Omega }=v/r`$ is the unperturbed angular velocity in the inertial frame, which depends only on $`r`$ in a barotropic disk, and

$$2B=\frac{1}{r}\frac{d}{dr}(r^2\mathrm{\Omega })$$
(9)
is the vertical vorticity, or twice Oort’s second parameter. The epicyclic frequency $`\kappa `$ is given by $`\kappa ^2=4\mathrm{\Omega }B`$. Also

$$h^{}=\frac{v_\mathrm{s}^2}{\rho }\rho ^{}+O(\rho _{}^{}{}_{}{}^{2}),$$
(10)
where $`v_\mathrm{s}^2=dp/d\rho `$ is the square of the sound speed.

### 3.2 Expansion about corotation

Let $`r_\mathrm{c}`$ be the corotation radius, defined by $`\mathrm{\Omega }(r_\mathrm{c})=\mathrm{\Omega }_\mathrm{p}`$. There are three characteristic length-scales that are potentially relevant for determining the radial extent of the corotation region,

$$\delta _c=\frac{c}{\kappa },\delta _\nu =\left(\frac{\nu }{md\mathrm{\Omega }/dr}\right)^{1/3},\delta _\mathrm{\Psi }=\left(\frac{\mathrm{\Psi }}{\kappa ^2d\mathrm{ln}\mathrm{\Omega }/d\mathrm{ln}r}\right)^{1/2}.$$
(11)
Here $`c`$ is an appropriate average of the sound speed $`v_\mathrm{s}`$, while $`m`$ and $`\mathrm{\Psi }`$ are the azimuthal wavenumber and amplitude of the forcing potential $`\mathrm{\Phi }^{}`$. We assume that $`d\mathrm{\Omega }/dr<0`$. In a linear, inviscid theory (GT79) the scale is set by the evanescent wavelength $`\delta _c`$ of the two-dimensional mode (density wave) near corotation. However, the solution involves a cusp, which may be expected to be smoothed over a characteristic length $`\delta _\nu `$ in the presence of an effective viscosity. On the other hand, a ballistic treatment of the corotation region suggests that an annulus of librating trajectories is formed, of half-width $`2\sqrt{2}\delta _\mathrm{\Psi }`$ (Goldreich & Tremaine 1981).

In order to allow all these effects to compete, we adopt scalings such that these three lengths are formally of the same order of magnitude. (Subsequently, we will take limits in which one or other of the quantities $`\delta `$ can be distinguished as the relevant scale.) The characteristic width of the corotation region is then $`H`$, the semi-thickness of the disk, as was also found by Lubow (1990). This allows us to resolve the radial structure of the region in the same way as for the vertical structure, by introducing scaled coordinates

$$\xi =\frac{rr_\mathrm{c}}{ϵ},\zeta =\frac{z}{ϵ},$$
(12)
where $`ϵ1`$ is a typical value of the angular semi-thickness $`H/r`$ of the disk. In units such that $`r_\mathrm{c}=1`$, the corotation region is found where $`\xi `$ and $`\zeta `$ are of order unity.

The central potential may be expanded in the corotation region in a double Taylor series,

$$\mathrm{\Phi }=(\mathrm{\Phi }_{00}+ϵ\mathrm{\Phi }_{01}\xi +\mathrm{})+ϵ^2\frac{1}{2}\zeta ^2(\mathrm{\Phi }_{20}+ϵ\mathrm{\Phi }_{21}\xi +\mathrm{})+\mathrm{}.$$
(13)
The perturbing potential may be expanded similarly, but is assumed formally to be smaller by a factor $`O(ϵ^2)`$:

$$\mathrm{\Phi }^{}/ϵ^2=[\mathrm{\Phi }_{00}^{}(\varphi )+ϵ\mathrm{\Phi }_{01}^{}(\varphi )\xi +\mathrm{}]+ϵ^2\frac{1}{2}\zeta ^2[\mathrm{\Phi }_{20}^{}(\varphi )+ϵ\mathrm{\Phi }_{21}^{}(\varphi )\xi +\mathrm{}]+\mathrm{}.$$
(14)
This assumes that both potentials are symmetric about $`z=0`$ and vary smoothly in $`r`$ in the neighborhood of $`r_\mathrm{c}`$. The scaling of $`\mathrm{\Phi }^{}`$ will turn out to provide a critical level of nonlinearity in the solution.

Properties of the unperturbed disk generally have non-trivial structure in $`\zeta `$ but may be expanded in Taylor series in $`r`$,

$`\rho `$ $`=`$ $`\rho _0(\zeta )+ϵ\rho _1(\zeta )\xi +\mathrm{},`$
$`u`$ $`=`$ $`ϵ^3\left[u_0(\zeta )+\mathrm{}\right],`$
$`w`$ $`=`$ $`ϵ^4\left[w_0(\zeta )+\mathrm{}\right],`$
$`\mathrm{\Omega }`$ $`=`$ $`\mathrm{\Omega }_0+ϵ\mathrm{\Omega }_1\xi +\mathrm{},`$
$`B`$ $`=`$ $`B_0+ϵB_1\xi +\mathrm{},`$
$`v_\mathrm{s}`$ $`=`$ $`ϵ\left[v_{\mathrm{s0}}(\zeta )+\mathrm{}\right],`$
$`\mu `$ $`=`$ $`ϵ^3\left[\mu _0(\zeta )+\mathrm{}\right],`$
$`\mu _\mathrm{b}`$ $`=`$ $`ϵ^3\left[\mu _{\mathrm{b0}}(\zeta )+\mathrm{}\right].`$ (15)
Note that $`\mathrm{\Omega }_0=\mathrm{\Omega }_\mathrm{p}`$ by definition of the corotation radius. The scaling of $`\mu `$ implies that the dimensionless viscosity parameter, $`\alpha =\mu \mathrm{\Omega }/p`$, is formally $`O(ϵ)`$, and will turn out to provide a critical level of dissipation in the solution. The scalings of the accretion flow, $`u`$ and $`w`$, follow from this assumption.

Perturbed variables generally have non-trivial structure in $`\xi `$ and $`\zeta `$. We have found the required scalings to be

$`\rho ^{}`$ $`=`$ $`ϵ\left[\rho _0^{}(\xi ,\varphi ,\zeta )+\mathrm{}\right],`$
$`u^{}`$ $`=`$ $`ϵ^2\left[u_0^{}(\varphi )+ϵu_1^{}(\xi ,\varphi ,\zeta )+\mathrm{}\right],`$
$`v^{}`$ $`=`$ $`ϵ^2\left[v_0^{}(\xi ,\varphi )+ϵv_1^{}(\xi ,\varphi ,\zeta )+\mathrm{}\right],`$
$`w^{}`$ $`=`$ $`ϵ^3\left[w_0^{}(\xi ,\varphi ,\zeta )+\mathrm{}\right],`$
$`h^{}`$ $`=`$ $`ϵ^3\left[h_0^{}(\xi ,\varphi )+ϵh_1^{}(\xi ,\varphi ,\zeta )+\mathrm{}\right],`$
$`\mu ^{}`$ $`=`$ $`ϵ^4\left[\mu _0^{}(\xi ,\varphi ,\zeta )+\mathrm{}\right],`$
$`\mu _\mathrm{b}^{}`$ $`=`$ $`ϵ^4\left[\mu _{\mathrm{b0}}^{}(\xi ,\varphi ,\zeta )+\mathrm{}\right].`$ (16)
Equations (4), (5) and (6) at leading order, $`O(ϵ^2)`$, then imply

$$(u_0^{}_\xi +\mathrm{\Omega }_1\xi _\varphi )\rho _0^{}+u_0^{}\rho _1+w_0^{}_\zeta \rho _0=\rho _0\left(\frac{u_0^{}}{r_\mathrm{c}}+_\xi u_1^{}+\frac{1}{r_\mathrm{c}}_\varphi v_0^{}+_\zeta w_0^{}\right),$$
(17)
$$2\mathrm{\Omega }_0v_0^{}=\mathrm{\Phi }_{01}^{}_\xi h_0^{},$$
(18)
$$2B_0u_0^{}=\frac{1}{r_\mathrm{c}}_\varphi \mathrm{\Phi }_{00}^{}.$$
(19)
Also required is equation (6) at the next order, $`O(ϵ^3)`$,

$$(u_0^{}_\xi +\mathrm{\Omega }_1\xi _\varphi )v_0^{}+2B_0u_1^{}+2B_1\xi u_0^{}=\frac{\xi }{r_\mathrm{c}}\left(_\varphi \mathrm{\Phi }_{01}^{}\frac{1}{r_\mathrm{c}}_\varphi \mathrm{\Phi }_{00}^{}\right)\frac{1}{r_\mathrm{c}}_\varphi h_0^{}+\frac{\mu _0}{\rho _0}_\xi ^2v_0^{}+\frac{r_\mathrm{c}\mathrm{\Omega }_1}{\rho _0}_\xi \mu _0^{}.$$
(20)
The last two terms are the only viscous terms to enter our analysis. Finally, equation (7) at leading order, $`O(ϵ^3)`$, requires only a hydrostatic balance,

$$0=\mathrm{\Phi }_{20}^{}\zeta _\zeta h_1^{}.$$
(21)
The enthalpy and density perturbations are related by equation (10) at leading order, $`O(ϵ^3)`$,

$$h_0^{}=\frac{v_{\mathrm{s0}}^2}{\rho _0}\rho _0^{}.$$
(22)
### 3.3 Derivation of the governing equation

We eliminate $`\rho _0^{}`$, $`u_1^{}`$, $`v_0^{}`$, and $`w_0^{}`$ (by integrating eq.  with respect to $`\zeta `$ over the full vertical extent of the disk) to obtain a single equation for $`h_0^{}`$,

$$(u_0^{}_\xi +\mathrm{\Omega }_1\xi _\varphi )(\kappa _0^2I_3I_1_\xi ^2)h_0^{}+(I_4_\xi ^2+2r_\mathrm{c}\mathrm{\Omega }_0\mathrm{\Omega }_1I_5)_\xi ^2h_0^{}=\frac{2\mathrm{\Omega }_0}{r_\mathrm{c}}(_\varphi \mathrm{\Phi }_{00}^{})\left(I_2\frac{B_1}{B_0}I_1\right),$$
(23)
involving the integrals

$$I_1=\rho _0𝑑\zeta ,I_2=\rho _1𝑑\zeta ,I_3=\frac{\rho _0}{v_{\mathrm{s0}}^2}𝑑\zeta ,I_4=\mu _0𝑑\zeta ,I_5=\frac{1}{h_0^{}}\mu _0^{}𝑑\zeta .$$
(24)
Here $`u_0^{}`$ is given directly by equation (19). We recognize $`I_1=\mathrm{\Sigma }_0`$ and $`I_2=\mathrm{\Sigma }_1`$ as coefficients in the Taylor expansion of the unperturbed surface density $`\mathrm{\Sigma }(r)`$. Defining the mean sound speed $`c`$, the mean kinematic viscosity $`\nu `$, and the logarithmic viscosity derivative $`D_\mu `$ in physical terms by

$$\frac{\mathrm{\Sigma }}{c^2}=\frac{\rho }{v_\mathrm{s}^2}𝑑z,\nu \mathrm{\Sigma }=\mu 𝑑z,D_\mu =\frac{\mathrm{ln}(\nu \mathrm{\Sigma })}{\mathrm{ln}\mathrm{\Sigma }},$$
(25)
we find $`I_3=\mathrm{\Sigma }_0/c_0^2`$, $`I_4=\nu _0\mathrm{\Sigma }_0`$, and $`I_5=D_\mu \nu _0\mathrm{\Sigma }_0/c_0^2`$.

When the $`ϵ`$-scalings and subscripts are omitted, the final equation in physical terms reads

$$\left(u^{}_\xi +\frac{d\mathrm{\Omega }}{dr}\xi _\varphi \right)\left(\frac{\kappa ^2}{c^2}_\xi ^2\right)h^{}+\nu \left(_\xi ^2+2r\mathrm{\Omega }\frac{d\mathrm{\Omega }}{dr}\frac{D_\mu }{c^2}\right)_\xi ^2h^{}=\frac{2\mathrm{\Omega }}{r}(_\varphi \mathrm{\Phi }^{})\frac{d}{dr}\mathrm{ln}\left(\frac{\mathrm{\Sigma }}{B}\right),$$
(26)
where

$$u^{}=\frac{1}{2rB}_\varphi \mathrm{\Phi }^{}.$$
(27)
Here, all quantities except $`\xi `$, $`\mathrm{\Phi }^{}(\varphi )`$, $`u^{}(\varphi )`$, and $`h^{}(\xi ,\varphi )`$ are to be understood as constants evaluated at $`r=r_\mathrm{c}`$.

Equation (26) describes, in a reduced but systematically obtained manner, how the enthalpy perturbation is forced by the external potential in the corotation region. If the advective term involving $`u^{}`$ were omitted, this would be a linear inhomogeneous equation, and the forcing would be proportional to the gradient of the vortensity $`\mathrm{\Sigma }/B`$. The $`u^{}`$ term derives from a product of perturbed quantities and is in this sense a nonlinear effect. An alternative way to write equation (26) is

$$\frac{d\mathrm{\Omega }}{dr}\xi _\varphi \left(\frac{\kappa ^2}{c^2}_\xi ^2\right)h^{}+\nu \left(_\xi ^2+2r\mathrm{\Omega }\frac{d\mathrm{\Omega }}{dr}\frac{D_\mu }{c^2}\right)_\xi ^2h^{}=\frac{2\mathrm{\Omega }}{r}(_\varphi \mathrm{\Phi }^{})\left[\frac{d}{dr}\mathrm{ln}\left(\frac{\mathrm{\Sigma }}{B}\right)+_\xi \left(\frac{\mathrm{\Sigma }^{}}{\mathrm{\Sigma }}\frac{B^{}}{B}\right)\right],$$
(28)
where $`\mathrm{\Sigma }^{}=\mathrm{\Sigma }h^{}/c^2`$ and $`2B^{}=_\xi v^{}`$ are the leading-order perturbations of the surface density and the vertical vorticity. This shows that the nonlinear effect can be understood either as an additional advection of the perturbation or as the feedback of the perturbation on the vortensity gradient of the disk.

### 3.4 Reduction to dimensionless form

We assume that the external potential contains a single Fourier component, $`\mathrm{\Phi }^{}=\mathrm{\Psi }\mathrm{cos}m\varphi `$, and rewrite the governing equation (26) in a dimensionless form by means of the transformations

$$\xi =x\frac{c}{\kappa },\varphi =\frac{\theta }{m},h^{}(\xi ,\varphi )=f(x,\theta )\frac{c^3}{\kappa }\frac{d}{dr}\mathrm{ln}\left(\frac{\mathrm{\Sigma }}{B}\right).$$
(29)
Thus

$$(a\mathrm{sin}\theta _x+x_\theta )(1_x^2)f\stackrel{~}{\nu }(_x^2b)_x^2f=a\mathrm{sin}\theta ,$$
(30)
where

$`a`$ $`=`$ $`2\left({\displaystyle \frac{d\mathrm{ln}r}{d\mathrm{ln}\mathrm{\Omega }}}\right){\displaystyle \frac{\mathrm{\Psi }}{c^2}},`$
$`b`$ $`=`$ $`2r\mathrm{\Omega }{\displaystyle \frac{d\mathrm{\Omega }}{dr}}{\displaystyle \frac{D_\mu }{\kappa ^2}},`$
$`\stackrel{~}{\nu }`$ $`=`$ $`\left({\displaystyle \frac{\nu }{md\mathrm{\Omega }/dr}}\right){\displaystyle \frac{\kappa ^3}{c^3}}.`$ (31)
The problem therefore depends on only three dimensionless parameters: the forcing amplitude $`a`$, the viscosity derivative $`b`$, and the dimensionless viscosity $`\stackrel{~}{\nu }`$. For a Keplerian disk, $`b=3D_\mu `$. In terms of the characteristic length-scales defined in equation (11), $`\stackrel{~}{\nu }=(\delta _\nu /\delta _c)^3`$.

Equation (30) is to be solved on $`\mathrm{}<x<\mathrm{}`$, $`0<\theta <2\pi `$ subject to a periodic boundary condition in $`\theta `$, which implies that the physical solution has the same periodicity in $`\varphi `$ as the forcing potential. The solution $`f`$ should be bounded as $`|x|\mathrm{}`$ in order to be a valid localized solution in the corotation region.

We note immediately that one non-localized solution of equation (30) is $`f=x+\mathrm{constant}`$. This represents an axisymmetric perturbation that exactly cancels the vortensity gradient of the unperturbed disk, so nullifying the corotation resonance. Such a solution should be excluded, as it amounts to redefining the original problem.

### 3.5 Solution by Fourier analysis

Equation (30) contains multiple $`x`$-derivatives but involves only one power of $`x`$ explicitly. Like Airy’s equation, it is amenable to Fourier analysis in $`x`$. Let

$$F(k,\theta )=_{\mathrm{}}^{\mathrm{}}f(x,\theta )e^{ikx}𝑑x$$
(32)
be the Fourier transform of $`f`$, which exists provided that $`f`$ decays sufficiently rapidly as $`|x|\mathrm{}`$. (In fact, by allowing for $`F`$ to be a generalized function, we may permit $`f`$ to grow at most algebraically in $`x`$.) The equation transforms to

$$\left[i_k_\theta +ika\mathrm{sin}\theta +\stackrel{~}{\nu }k^2\left(\frac{b+k^2}{1+k^2}\right)\right]G=2i\mathrm{sin}\theta \delta (k),$$
(33)
where

$$G(k,\theta )=\frac{1}{i\pi a}(1+k^2)F(k,\theta ).$$
(34)
Equation (33) could be solved numerically as a partial differential equation, but is amenable to further Fourier analysis in $`\theta `$. Let

$$G(k,\theta )=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}G_n(k)e^{in\theta },$$
(35)
then we obtain the system of ordinary differential equations

$$n\frac{dG_n}{dk}+\stackrel{~}{\nu }k^2\left(\frac{b+k^2}{1+k^2}\right)G_n\frac{ka}{2}\left(G_{n+1}G_{n1}\right)=\left(\delta _{n,1}\delta _{n,1}\right)\delta (k).$$
(36)
Note that the symmetry property

$$G_n(k)=\left[G_n(k)\right]^{}$$
(37)
follows from the fact that the enthalpy perturbation $`h^{}`$ is real. As a boundary condition, $`G_n(k)`$ is required to tend to zero as $`|k|\mathrm{}`$, in order that it be the Fourier transform of a continuous function.

The total tidal torque exerted on the disk is

$$T=_{\mathrm{}}^{\mathrm{}}_0^{2\pi }_0^{\mathrm{}}(\rho +\rho ^{})_\varphi \mathrm{\Phi }^{}rdrd\varphi dz.$$
(38)
When expressed in terms of the scaled variables introduced in Section 3.4, the torque exerted in the corotation region is

$$T_\mathrm{c}=\frac{mrc^2\mathrm{\Psi }}{4\mathrm{\Omega }}\frac{d}{dr}\left(\frac{\mathrm{\Sigma }}{B}\right)_0^{2\pi }_{\mathrm{}}^{\mathrm{}}f(x,\theta )\mathrm{sin}\theta dxd\theta ,$$
(39)
where the prefactor is understood to be evaluated at $`r=r_\mathrm{c}`$. In terms of the Fourier representation, this further simplifies to

$$T_\mathrm{c}=t_\mathrm{c}T_{\mathrm{GT}},$$
(40)
where

$$t_\mathrm{c}=G_1(0)G_1(0)$$
(41)
is dimensionless, and

$$T_{\mathrm{GT}}=\frac{m\pi ^2\mathrm{\Psi }^2}{2(d\mathrm{\Omega }/dr)}\frac{d}{dr}\left(\frac{\mathrm{\Sigma }}{B}\right)$$
(42)
is the torque formula given by GT79. Although equation (36) requires $`G_{\pm 1}(k)`$ to have unit discontinuities at $`k=0`$, the dimensionless torque $`t_\mathrm{c}`$ is nevertheless well defined as the quantity $`G_1(k)G_1(k)`$ is continuous at $`k=0`$.

### 3.6 A simplifying assumption

A considerable simplification in equation (36) is obtained when $`b=1`$. This occurs in a Keplerian disk when conditions in the corotation region are such that $`\nu \mathrm{\Sigma }\mathrm{\Sigma }^{1/3}`$ locally. In the absence of any detailed knowledge of the effective viscosity, we adopt this convenient assumption, anticipating that our results will not depend sensitively on it.

It is then natural to rescale the wavenumber to

$$\stackrel{~}{k}=\stackrel{~}{\nu }^{1/3}k,$$
(43)
so that

$$n\frac{dG_n}{d\stackrel{~}{k}}+\stackrel{~}{k}^2G_np\stackrel{~}{k}\left(G_{n+1}G_{n1}\right)=\left(\delta _{n,1}\delta _{n,1}\right)\delta (\stackrel{~}{k}),$$
(44)
where

$$p=\frac{1}{2}a\stackrel{~}{\nu }^{2/3}$$
(45)
is the only remaining parameter of the problem. In physical terms,

$$p=\left(\frac{\mathrm{\Psi }}{\kappa ^2d\mathrm{ln}\mathrm{\Omega }/d\mathrm{ln}r}\right)\left(\frac{md\mathrm{\Omega }/dr}{\nu }\right)^{2/3}$$
(46)
measures the strength of the potential relative to the viscosity (i.e., the nonlinearity relative to the dissipation). In terms of the three characteristic length-scales defined in equation (11), $`p=(\delta _\mathrm{\Psi }/\delta _\nu )^2`$.

We may assume without loss of generality that $`p>0`$. Our aim now is to evaluate the function $`t_\mathrm{c}(p)`$ that determines how the formula of GT79 is modified.

### 3.7 Linear theory

In linear theory the coupling parameter $`p`$ is set to zero and only modes $`n=\pm 1`$ are excited. The solution is

$$G_{\pm 1}(\stackrel{~}{k})=\pm H(\pm \stackrel{~}{k})\mathrm{exp}(\frac{1}{3}\stackrel{~}{k}^3),$$
(47)
where $`H`$ is the Heaviside function. The dimensionless torque (eq. ) is then $`t_\mathrm{c}=1`$. This shows that the torque formula of GT79, which was derived for a two-dimensional, inviscid disk, also applies to a three-dimensional, barotropic and viscous disk, provided that the linearization is valid (i.e., $`p1`$).

### 3.8 Feedback

The axisymmetric component $`G_0`$ of the solution has a special role. The case $`n=0`$ in equation (44) is not a differential equation but states that

$$\stackrel{~}{k}^2G_0=p\stackrel{~}{k}(G_1G_1).$$
(48)
Now $`G_0`$ also appears in the equations for $`G_{\pm 1}`$ in the form $`p\stackrel{~}{k}G_0`$. We solve equation (48) within the space of generalized functions to obtain

$$\stackrel{~}{k}G_0=p(G_1G_1)+C_1\delta (\stackrel{~}{k}),$$
(49)
where $`C_1`$ is an arbitrary constant. Equation (44) with $`n=\pm 1`$ then implies

$$\frac{dG_1}{d\stackrel{~}{k}}+\stackrel{~}{k}^2G_1p\stackrel{~}{k}G_2+p^2(G_1G_1)=(1C_1p)\delta (\stackrel{~}{k}),$$
(50)
$$\frac{dG_1}{d\stackrel{~}{k}}\stackrel{~}{k}^2G_1p\stackrel{~}{k}G_2+p^2(G_1G_1)=(1C_1p)\delta (\stackrel{~}{k}).$$
(51)
The solution of equation (49) is of the form

$$G_0=G_0^{\mathrm{reg}}(\stackrel{~}{k})+\frac{pt_\mathrm{c}}{\stackrel{~}{k}}C_1\delta ^{}(\stackrel{~}{k})+iC_2\delta (\stackrel{~}{k}),$$
(52)
where $`G_0^{\mathrm{reg}}(\stackrel{~}{k})`$ is regular at $`\stackrel{~}{k}=0`$, and $`C_2`$ is a second arbitrary constant. Consider the inverse Fourier transforms

$$_{\mathrm{}}^{\mathrm{}}k^1e^{ikx}𝑑k=i\pi \mathrm{sgn}(x),_{\mathrm{}}^{\mathrm{}}\delta ^{}(k)e^{ikx}𝑑k=ix,_{\mathrm{}}^{\mathrm{}}i\delta (k)e^{ikx}𝑑k=i,$$
(53)
the first of which is a principal value. By allowing $`C_10`$ we would admit a perturbation that redefined the basic vortensity gradient of the disk across the whole corotation region. This should be excluded, as discussed in Section 3.4. We may also set $`C_2=0`$ without loss of generality, as it simply adds a constant to the basic surface density.

The solution generated by the $`\stackrel{~}{k}^1`$ singularity in $`G_0`$ corresponds to a gradual, axisymmetric step in the enthalpy perturbation across the corotation region, in the sense that $`f(x,\theta )`$ tends to well defined limits as $`x\pm \mathrm{}`$, but these limiting values are unequal. This occurs because the tidal torque causes a redistribution of the surface density of the disk. Consider the standard evolutionary equations for an accretion disk in which the corotation torque is added as an infinitely localized source,

$$\frac{\mathrm{\Sigma }}{t}+\frac{1}{r}\frac{}{r}(r\mathrm{\Sigma }v_r)=0,$$
(54)
$$\mathrm{\Sigma }v_r\frac{d}{dr}(r^2\mathrm{\Omega })=\frac{1}{r}\frac{}{r}\left(\nu \mathrm{\Sigma }r^3\frac{d\mathrm{\Omega }}{dr}\right)+\frac{T_\mathrm{c}}{2\pi r_\mathrm{c}}\delta (rr_\mathrm{c}),$$
(55)
where $`v_r`$ is the mean radial velocity. In a steady state these integrate to

$$2\pi r\mathrm{\Sigma }v_r=\dot{M}=\mathrm{constant},$$
(56)
$$\dot{M}r^2\mathrm{\Omega }+2\pi \nu \mathrm{\Sigma }r^3\frac{d\mathrm{\Omega }}{dr}+T_\mathrm{c}H(rr_\mathrm{c})=\mathrm{constant}.$$
(57)
The step in $`\nu \mathrm{\Sigma }`$ across the resonant region associated with the $`\stackrel{~}{k}^1`$ singularity in $`G_0`$ can be shown exactly to balance the corotation torque $`T_\mathrm{c}`$ in this equation. Note that this simplified description does not treat the internal structure of the corotation region, but is intended to show how the torque affects the disk on a larger scale.

### 3.9 Approximate solutions for small and large $`p`$

Approximate solutions can be obtained by neglecting $`G_2`$ and $`G_2`$ in equations (50) and (51). This amounts to a severe Fourier truncation of the problem, in which we allow for the solution to feed back on the axisymmetric part of the disk but neglect the excitation of higher harmonics. Nevertheless, this method is found to give results in close agreement with the more accurate numerical solutions described in the next section.

Defining $`G_\pm =G_1\pm G_1`$ and taking the sum and difference of equations (50) and (51), we obtain

$$\frac{dG_+}{d\stackrel{~}{k}}+(\stackrel{~}{k}^2+2p^2)G_{}=2\delta (\stackrel{~}{k}),$$
(58)
$$\frac{dG_{}}{d\stackrel{~}{k}}+\stackrel{~}{k}^2G_+=0,$$
(59)
and so

$$\frac{d}{d\stackrel{~}{k}}\left(\frac{1}{\stackrel{~}{k}^2}\frac{dG_{}}{d\stackrel{~}{k}}\right)+(\stackrel{~}{k}^2+2p^2)G_{}=2\delta (\stackrel{~}{k}).$$
(60)
For $`p^21`$ the solution may be expanded in powers of $`p^2`$,

$$G_{}=G_{}^{(0)}+p^2G_{}^{(1)}+O(p^4).$$
(61)
Using the fact that the bounded solution of the equation

$$\frac{d^2y}{dx^2}+y=2f(x)$$
(62)
is

$$y(x)=_{\mathrm{}}^{\mathrm{}}f(x^{})e^{|xx^{}|}𝑑x^{},$$
(63)
we find, at successive orders in $`p^2`$,

$$G_{}^{(0)}(\stackrel{~}{k})=\mathrm{exp}(\frac{1}{3}|\stackrel{~}{k}^3|),$$
(64)
$$G_{}^{(1)}(\stackrel{~}{k})=_{\mathrm{}}^{\mathrm{}}\mathrm{exp}(\frac{1}{3}|\stackrel{~}{k}^3|\frac{1}{3}|\stackrel{~}{k}^3\stackrel{~}{k}^3|)𝑑\stackrel{~}{k}^{}.$$
(65)
The dimensionless torque is then

$$t_\mathrm{c}=G_{}(0)=1\left(\frac{2}{3}\right)^{2/3}\mathrm{\Gamma }\left(\frac{1}{3}\right)p^2+O(p^4)12.044p^2.$$
(66)
This asymptotic form can be shown to be correct even when the coupling to higher harmonics is taken into account. In this limit of small $`p`$, only wavenumbers $`|\stackrel{~}{k}|1`$ contribute significantly to the solution.

For $`p^21`$ we may neglect $`\stackrel{~}{k}^2`$ relative to $`2p^2`$ in equation (60). The solution is then

$$G_{}(\stackrel{~}{k})\frac{2^{1/8}}{\mathrm{\Gamma }\left(\frac{1}{4}\right)}p^{3/4}|\stackrel{~}{k}|^{3/2}K_{3/4}\left(2^{1/2}p\stackrel{~}{k}^2\right),$$
(67)
where $`K`$ is the modified Bessel function of the second kind, and the dimensionless torque is

$$t_\mathrm{c}=G_{}(0)\frac{\mathrm{\Gamma }\left(\frac{3}{4}\right)}{\mathrm{\Gamma }\left(\frac{1}{4}\right)}2^{1/4}p^{3/2}0.4019p^{3/2}.$$
(68)
Since $`p\nu ^{2/3}`$, this implies that the torque is proportional to the viscosity when the viscosity is small. This is to be expected on general grounds, but to neglect the coupling to higher harmonics is difficult to justify formally when $`p`$ is large. The above form of $`G_{}`$ indicates that only wavenumbers $`|\stackrel{~}{k}|p^{1/2}`$ contribute significantly.

### 3.10 Numerical solution

For a numerical solution, we truncate the system at some order $`N`$ and set $`G_n`$ to zero for $`|n|>N`$. We then solve equations (50) and (51) together with equation (44) for $`2|n|N`$. Now the symmetry property (37) allows us to consider $`\stackrel{~}{k}>0`$ only, and the jump conditions at $`\stackrel{~}{k}=0`$ imply

$$G_n(0+)+G_n(0+)=\delta _{n,1},n>0.$$
(69)
We integrate the equations from $`\stackrel{~}{k}=0+`$ to a finite value $`\stackrel{~}{k}=\stackrel{~}{k}_{\mathrm{max}}`$. A ‘particular solution’ $`G_n^{(\mathrm{p})}`$ satisfying the inhomogeneous boundary condition (69) is generated by starting from the initial condition $`G_1=1`$ ($`G_n=0`$ otherwise), and $`N`$ ‘complementary functions’ $`G_n^{(\mathrm{c},q)}`$ by starting from initial conditions $`G_q=G_q=1`$ ($`G_n=0`$ otherwise) for each $`q=1,\mathrm{},N`$. The amplitudes of the complementary functions appearing in the desired solution are determined by requiring that $`G_n(\stackrel{~}{k}_{\mathrm{max}})=0`$ for $`Nn1`$. This simulates the requirement that $`G_n`$ should tend to zero as $`|\stackrel{~}{k}|\mathrm{}`$. For $`n<0`$, $`G_n`$ has a tendency to grow as $`\mathrm{exp}(\frac{1}{3}\stackrel{~}{k}^3/|n|)`$ as $`\stackrel{~}{k}\mathrm{}`$, and our boundary condition is designed to eliminate such growing solutions. For $`n>0`$, $`G_n`$ tends to decay naturally as $`\stackrel{~}{k}\mathrm{}`$.

In Figure 1 we plot the dimensionless torque $`t_\mathrm{c}`$ determined from the numerical solution as a function of the coupling parameter $`p`$. The solution converges rapidly as $`N`$ and $`\stackrel{~}{k}_{\mathrm{max}}`$ are increased, and good agreement is found with the approximate solutions (66) and (68) for small and large $`p`$.

### 3.11 Interpretation

We now attempt to give a real-space interpretation of these results. Multiplying our governing equation (26) by $`\mathrm{\Sigma }r/(2\mathrm{\Omega })`$ and integrating over $`\varphi `$, we obtain

$$_\xi ^2\left\{\mathrm{\Sigma }\left[\frac{1}{\kappa ^2}(_\varphi \mathrm{\Phi }^{})_\xi h^{}+\frac{\nu r}{2\mathrm{\Omega }}\left(_\xi ^2+2r\mathrm{\Omega }\frac{d\mathrm{\Omega }}{dr}\frac{D_\mu }{c^2}\right)h^{}\right]\right\}d\varphi =\frac{\mathrm{\Sigma }}{c^2}(_\varphi \mathrm{\Phi }^{})_\xi h^{}d\varphi ,$$
(70)
where the very first term has been taken over to the right-hand side. The first integral with respect to $`\xi `$ is

$$_\xi \left\{\mathrm{\Sigma }\left[\frac{1}{\kappa ^2}(_\varphi \mathrm{\Phi }^{})_\xi h^{}+\frac{\nu r}{2\mathrm{\Omega }}\left(_\xi ^2+2r\mathrm{\Omega }\frac{d\mathrm{\Omega }}{dr}\frac{D_\mu }{c^2}\right)h^{}\right]\right\}d\varphi =\frac{\mathrm{\Sigma }}{c^2}h^{}_\varphi \mathrm{\Phi }^{}d\varphi ,$$
(71)
plus a possible constant of integration.

Equation (71) expresses the conservation of angular momentum in a steady state, although the terms that balance in the unperturbed disk do not appear. The left-hand side is the divergence of the angular momentum flux, and the right-hand side is the tidal torque per unit radius. To see this, the first term in the flux can be identified as the integrated Reynolds stress

$$r\rho u^{}v^{}𝑑\varphi 𝑑z,$$
(72)
using equations (18) and (19). The remaining terms in the flux are viscous stresses. The right-hand side of equation (71) can also be written as

$$\rho ^{}_\varphi \mathrm{\Phi }^{}d\varphi dz,$$
(73)
which is the tidal torque per unit radius.

In the Fourier representation, when suitably normalized, equation (71) reads

$$\frac{k^2a}{2}(F_1F_1)\stackrel{~}{\nu }k(k^2+b)F_0=\frac{a}{2}(F_1F_1),$$
(74)
plus a possible multiple of $`\delta (k)`$. The same equation could be written for $`G`$ instead of $`F`$, and is then equivalent to our equation (36) in the case $`n=0`$. (The possible additional term proportional to $`\delta (k)`$ is equivalent to the $`C_1`$ term that arises in passing from equation (48) to (49), and can be excluded on the same grounds.) This shows how equation (36) partitions into terms associated with Reynolds stress, viscous torque and tidal torque.

In Section 3.9 we worked with the quantity $`G_{}(1+k^2)(F_1F_1)`$. In these terms, the Fourier transforms of the Reynolds, viscous and tidal terms are proportional to

$$\left(\frac{k^2}{1+k^2}\right)G_{}(k),G_{}(k),\left(\frac{1}{1+k^2}\right)G_{}(k),$$
(75)
respectively. When $`p1`$ we found that $`G_{}(k)`$ is cut off for wavenumbers $`|\stackrel{~}{k}|1`$, i.e. physical wavenumbers $`|k/\delta _c|1/\delta _\nu `$. It follows that the characteristic scale of the viscous torque is the viscous scale $`\delta _\nu `$ defined in equation (11). If the viscosity is small, so that $`\delta _\nu \delta _c`$, then the scale of the Reynolds stress is also $`\delta _\nu `$. However, the scale of the tidal torque is $`\delta _c`$, because its Fourier transform is cut off for physical wavenumbers $`|k/\delta _c|1/\delta _c`$ by the factor $`(1+k^2)^1`$.

When $`p1`$ we found that $`G_{}(k)`$ is cut off for wavenumbers $`|\stackrel{~}{k}|p^{1/2}`$, i.e. physical wavenumbers $`|k/\delta _c|1/\delta _\mathrm{\Psi }`$. If the forcing amplitude is small, so that $`\delta _\mathrm{\Psi }\delta _c`$, then the scale of the Reynolds stress is also $`\delta _\mathrm{\Psi }`$, but the scale of the tidal torque is again $`\delta _c`$.

To summarize, if $`\delta _c`$ is the largest of the three scales, then the tidal torque is spread over a region of scale $`\delta _c`$, while the Reynolds and viscous stresses have a smaller scale $`\mathrm{max}(\delta _\nu ,\delta _\mathrm{\Psi })`$. This is consistent with the impression obtained from the original theory (GT79), that the tidal torque creates a disturbance in the disk on the scale $`\delta _c`$, which subsequently transfers its angular momentum to the background disk, through dissipation or nonlinearity, on a finer scale.

## 4 Application to eccentric resonances

The reduction of the torque exchanged between the perturber and the disk has consequences for the eccentricity evolution of young planets orbiting in a protoplanetary disk. We consider the first-order eccentric corotation resonances associated with a planet of mass ratio $`q=M_\mathrm{p}/M_{}`$ executing a orbit of eccentricity $`e`$ within a Keplerian disk.

GT80 provide expressions for the locations of these resonances and the forcing potentials. Inner and outer resonances occur at radii

$$r=\left(\frac{m}{m\pm 1}\right)^{2/3}a_\mathrm{p},$$
(76)
where $`a_\mathrm{p}`$ is the semi-major axis of the planet’s orbit. The amplitude of the forcing potential may be written in the form

$$\frac{\mathrm{\Psi }}{r^2\mathrm{\Omega }^2}=0.8020C_m^\pm meq,$$
(77)
where $`C_m^\pm `$ is a correction factor that tends to unity for large $`m`$. Accordingly, the saturation parameter of our analysis is

$$p=0.7006C_m^\pm m^{5/3}eq\alpha ^{2/3}\left(\frac{H}{r}\right)^{4/3},$$
(78)
where we write $`\nu =\alpha cH`$ and $`H=c/\mathrm{\Omega }`$. When $`p=1`$, i.e., when

$$e=e_{0.63}=\frac{1.427}{C_m^\pm }m^{5/3}q^1\alpha ^{2/3}\left(\frac{H}{r}\right)^{4/3},$$
(79)
$`63\%`$ saturation occurs because $`t_\mathrm{c}0.37`$. For $`5\%`$ saturation, only $`p0.16`$ is required.

Table 1 contains results for the lowest-order outer eccentric resonances ($`l=m1`$ in the notation of GT80) for a disk with $`H/r=0.05`$ and a planet of mass ratio $`q=0.001`$. The first column labels the azimuthal wavenumber of the potential component being considered. The second column gives the location of the resonance in units of $`a_\mathrm{p}`$. The third column contains the correction factor defined in equation (77). The fourth column gives a measure of the relative importance of the resonances: neglecting the torque cut-off and any differences in the disk parameters between different resonant radii, the unsaturated torque scales $`m\mathrm{\Psi }^2C_m^2m^3`$. The fifth column, labeled $`e_{0.63}`$, is the eccentricity required to reduce the eccentric corotation torque $`63\%`$ below its unsaturated value for a turbulent viscosity parameter $`\alpha =0.004`$. The sixth column, $`e_{0.05}`$, is similar to the fifth column, but gives the eccentricity needed to reduce the torque $`5\%`$ below its unsaturated level. The final column, $`\alpha _{0.63}`$, is the value of $`\alpha `$ corresponding to $`63\%`$ saturation when the eccentricity is $`e=0.1`$. Table 2 gives the same data for the inner resonances ($`l=m+1`$).

Resonances of high order occur very close to the planet’s orbit, and are present only if the planet does not clear a substantial gap in the disk. An analysis of such resonances requires a consideration of a number of effects that we have neglected. When $`m`$ becomes comparable to $`r/H`$, additional terms in the dynamical equations must be included, and the torque is reduced from the standard value (GT79). The proximity of the resonance to the planet also means that the perturbing potential varies significantly with $`z`$ within the disk, and the relevant quantity is the density-weighted vertical average of the potential, rather than the mid-plane value (Ward 1986, 1989). These two effects apply equally to Lindblad and corotation resonances, and so do not alter the $`5\%`$ balance between torques that excite and damp eccentricity. A further effect is that neighboring resonances of sufficiently high order overlap one another, which may alter the resonant torques.

## 5 Summary and discussion

We have determined the torque exerted in a steady state by an external potential on a three-dimensional gaseous disk at a non-coorbital corotation resonance. Our model accounts for the feedback of the torque on the surface density and vorticity in the corotation region, and assumes that the disk has a barotropic equation of state and a nonzero effective viscosity. The torque formula of Goldreich & Tremaine (1979) (eq. ) must be modified by a reduction factor $`t_\mathrm{c}(p)`$, plotted in Figure 1, which quantifies the extent to which the resonance is saturated. This factor depends essentially on a single dimensionless parameter $`p`$, defined in equation (46), which measures the strength of the potential $`\mathrm{\Psi }`$ relative to the viscosity $`\nu `$ (i.e., the nonlinearity relative to the dissipation).

In Section 3.10, we determined that the characteristic radial width of the resonance in the limit of large $`p`$ (low viscosity) is $`\delta \sqrt{\mathrm{\Psi }}/\mathrm{\Omega }`$. (This refers to the scale of Reynolds and viscous stresses; the tidal torque is spread over a region of characteristic width $`c/\kappa `$.) We found that the torque is reduced by a factor $`p^{3/2}`$, which can be regarded as the ratio of the viscous diffusion rate $`\nu /\delta ^2`$ over the libration region to the libration rate $`m|d\mathrm{\Omega }/dr|\delta `$. In this regime, our results are broadly consistent with the saturation model of Ward (1992), which involves the same ratio of rates, but applied to the coorbital region. Earlier, Goldreich & Tremaine (1981) had argued that, in a disk of collisional particles, saturation would occur when (in our notation) $`p1`$.

Goldreich & Sari (2002) have recently developed an evolutionary model for planetary eccentricity in which eccentricity growth occurs through a finite-amplitude instability. For infinitesimal eccentricity, the damping caused by corotation resonances just overcomes the eccentricity growth due to Lindblad resonances. Above a critical level of eccentricity, however, the corotation resonances become sufficiently saturated that growth occurs. Using our evaluation of the saturation function $`t_\mathrm{c}(p)`$, Goldreich & Sari (2002) were able to determine the critical value of eccentricity required within the context of a specific disk model. The results in Section 4 and Table 1 show that for a Jupiter-mass planet orbiting within a typical protoplanetary disk, eccentricities of a few percent are adequate to saturate all first-order eccentric corotation resonances except those of the lowest $`m`$-values ($`m4`$). To achieve a $`5\%`$ reduction in torque at such resonances, as may be sufficient to change the balance from eccentricity damping to growth (GT80), requires only eccentricities of $`1\%`$ or less.

Exactly which resonances contribute most to eccentricity evolution depends on the extent of the gap cleared by the planet, which in turn depends on the mass ratio and the properties of the disk. As seen in Table 1, the largest corotation torques are those associated with the highest $`m`$-values, and they are the easiest to saturate. A lower-mass planet that opens a smaller gap may excite many resonances so that the dominant torque comes from approximately the cut-off $`m`$-value, of order $`r/H`$ (GT80). In that case, $`63\%`$ saturation is achieved optimistically at the torque cut-off when $`e1.4q^1\alpha ^{2/3}(H/r)^3`$. Note that, although the large value of $`m`$ is advantageous for saturation, the low mass of the planet is unfavorable (see eq. ). On the other hand, if the planet mass is small enough so that there is no gap in the disk, then the coorbital resonances cause eccentricity decay (Ward 1988; Artymowicz 1993).

We have neglected a number of potential complications. Numerical simulations of a Jupiter-mass planet in a circular orbit with the disk parameters adopted in Section 4 suggest that the disk edge is not sharp (i.e., its radial extent is much greater than $`H`$; see Figure 1 of Lubow, Seibert, & Artymowicz 1999). It is evident, however, that within the radial extent of the planet’s Roche lobe, material is captured by the planet and the resonant effects considered here do not play a role. It is not clear which eccentric corotation resonances are excited in the presence of an eccentric planet. The complications are due at least in part to nonlinear effects other than those considered here. For example, shocks associated with the planet’s wake cause a non-closure of streamlines in the vicinity of the disk edge, leading to a drift of material through the corotation regions. In addition, material within the Roche lobe of the planet exerts torques that may contribute to the eccentricity balance.

As usual in accretion disk theory, our knowledge of the effective viscosity of the disk is limited. In Section 3.6 we adopted a convenient assumption in order to make analytical progress. Although we have succeeded in analyzing a three-dimensional disk, thereby generalizing existing theories, we assumed that the disk was barotropic. The effects of buoyancy or baroclinicity on the corotation region remain uninvestigated.

In recent work we were able to verify analytical theories of the torques exerted at Lindblad and vertical resonances through direct numerical simulations (Bate et al. 2002). It would be valuable to conduct simulations of a non-coorbital corotation resonance to test the findings of the present paper.

After eccentric corotation resonances saturate at small eccentricity, the eccentricity growth may be limited at intermediate values. The limitation may be due to the overlap of resonances or alternatively through the excitation of higher order Lindblad resonances, some of which cause eccentricity damping. A contribution to eccentricity damping occurs at an eccentric inner (outer) Lindblad resonance that lies outside (inside) the planet’s orbit. SPH simulations of eccentric-orbit binary stars suggest that little eccentricity growth via resonances occurs for eccentricities in the range of $`0.50.7`$ or higher (Lubow & Artymowicz 1993), and similar limits may occur for planets.

In conclusion, we have found a simple quantitative measure of the saturation of a corotation resonance in a gaseous disk. This analysis suggests that planets may plausibly experience a net growth of eccentricity through their interaction with the disk in a variety of circumstances, provided that the eccentricity is not extremely small to begin with.

We are grateful to Peter Goldreich and Re’em Sari for allowing us to coordinate our presentation of results with theirs prior to publication, and also for correcting some inaccuracies and omissions in an earlier version of this paper. We also thank the anonymous referee for a very thorough report that led to an improvement of the paper. This work was initiated at the Aspen Center for Physics, and we are grateful for their hospitality. GIO acknowledges the support of the Royal Society through a University Research Fellowship. SHL acknowledges support from NASA grant NAG5-10732.


format: text

text:
# Creating superpositions that correspond to efficiently integrable probability distributions

## Abstract

We give a simple and efficient process for generating a quantum superposition of states which form a discrete approximation of any efficiently integrable (such as log concave) probability density functions.

It is well known that probabilistic algorithms are sometimes more powerful (and often more practical!) than their deterministic counterparts in classical computing. Such algorithms generically make use of some source of random bits. It has been shown aharonov  however, that giving a quantum computer access to random bits does not, in fact, further empower it.

There has been some attention recently to the problem of using quantum computers to generate certain probability distributions, in particular by looking at the properties of quantum random walks aharonov2 .

In this short note, we will look at the question of whether, given a certain probability distribution $`\{p_i\}`$, we can efficiently create a quantum superposition of the form

$$|\psi (\{p_i\})=\underset{i}{}\sqrt{p_i}|i,$$
(1)
where the $`|i`$ are an orthonormal set of states. We are, of course, particularly interested in the case where the index $`i`$ ranges over some exponentially large set of possibilities $`N`$.

We have not solved the general problem of when a state of the form (1) can be efficiently created. However we will show here that if an efficient classical algorithm exists to *integrate* a certain probability distribution $`p(x)`$, then we can efficiently create a state of the form (1) for the discretized version $`\{p_i\}`$ of $`p(x)`$. A well known set of probability density functions which are efficiently integrable by monte carlo methods are *log-concave* distributions. (A log-concave distribution is one for which $`\frac{^2\mathrm{log}p(x)}{x^2}<0`$.) Most important statistical distributions (including the exponential and normal families) are log-concave. The task of sampling and/or integrating such distributions frequently arises for problems of determining the (multidimensional) volumes of solids applegate .

Clearly if an efficient classical algorithm exists, then we can certainly create a quantum computer in the mixed state $`\rho =_ip_i|ii|`$. Our interest in creating a state of the form (1) is thus not simply to be able to measure the state $`|i`$ with probability $`p_i`$. Rather we hope that generating the distribution in a coherent fashion will allow *further* (uniquely quantum) processing of this distribution. At the end we discuss some possibilities along these lines.

We will illustrate the procedure considering only a probability distribution over a single random variable, although the method can be easily extended to multivariable distributions.

We let $`n=\mathrm{log}N`$, where $`N`$ is the total number of points over which we wish to discretize this distribution. Imagine the distribution is divided into some number $`2^m`$ of regions, and that we already have a $`m`$ qubit state

$$|\psi _m=\underset{i=0}{\overset{2^m1}{}}\sqrt{p_i^{(m)}}|i,$$
(2)
where, as usual, the integer $`i`$, is interpreted in binary as the tensor product sequence of qubit states. Here $`p_i^{(m)}`$ is the probability for the random variable $`x`$ to lie region $`i`$; thus $`p_0^{(m)}`$ is the probability for $`x`$ to lie in the far leftmost region, $`p_1^{(m)}`$ is the probability for it to lie in the region adjacent to this and so on. Our goal is to show that we can now efficiently subdivide these $`2^m`$ regions to yield a $`2^{m+1}`$ region discretization of $`p(x)`$. That is we wish to add one qubit to the state (2), such that we achieve the evolution

$$\sqrt{p_i^{(m)}}|i\sqrt{\alpha _i}|i|0+\sqrt{\beta _i}|i|1,$$
(3)
where $`\alpha `$ (resp. $`\beta `$) is the probability for $`x`$ to lie in the *left* (resp. *right*) half of region $`i`$. If we can achieve such an evolution, then the new state of our $`m+1`$ qubits is

$$|\psi _{m+1}=\underset{i=0}{\overset{2^{m+1}1}{}}\sqrt{p_i^{(m+1)}}|i.$$
(4)
This process is then repeated until $`m=n`$, i.e. until we have created the desired superposition (1) over all $`N=2^n`$ states. If such an evolution is possible, then clearly it is efficient.

To show how an evolution of the form (3) can be achieved, we first define $`x_L^i`$ and $`x_R^i`$ to be the left and right boundaries of region $`i`$. We then define the function

$$f(i)=\frac{_{x_L^i}^{\frac{x_R^ix_L^i}{2}}p(x)𝑑x}{_{x_L^i}^{x_R^i}p(x)𝑑x}.$$

$`f(i)`$ is simply the probability that, given $`x`$ lies in region $`i`$, it also lies in the left half of this region. Since $`_a^bp(x)𝑑x`$ is efficiently computable classically, we can take an ancilla register initially in the state $`|0\mathrm{}0`$, and construct a circuit which efficiently performs the computation

$$\sqrt{p_i^{(m)}}|i|0\mathrm{}0\sqrt{p_i^{(m)}}|i|\theta _i,$$
(5)
where $`\theta _i\mathrm{arccos}\sqrt{f(i)}`$. We now perform a controlled rotation of angle $`\theta _i`$ on the $`m+1`$’th qubit:

$$\sqrt{p_i^{(m)}}|i|\theta _i|0\sqrt{p_i^{(m)}}|i|\theta _i(\mathrm{cos}\theta _i|0+\mathrm{sin}\theta _i|1)$$
(6)
and we then uncompute the register containing $`|\theta _i`$ to leave us in a state of the form (3) as desired. Note that the efficient classical algorithm to perform the integration may well be probabilistic, however it can always be implemented on a quantum computer in the required coherent manner by using an ancilla register initially prepared containing the string of random bits to be used in the computation.

Having seen that we can create a state of the form (1), we now give some examples of why creating such a state might be interesting.

(i) *Non-uniform priors for quantum searching* The quantum search algorithm lov  is typically envisioned as a search over a large number of possible solutions to a problem, where every potential solution is a-priori equally likely. In many instances however, we might believe that the correct solution is more likely to be found in a certain region of state space than another. Thus we might wish to input a non-uniform prior distribution into the search algorithm. If this prior distribution is log-concave (e.g. gaussian, poissonian etc) then the procedure outlined above can be used to generate it.

(ii) *Sampling of non-log-concave distributions* A state of the form (1) is clearly of little value if we were to just measure the qubits and obtain outcome $`|i`$ with probability $`p_i`$, since such sampling is possible classically. The fundamental difference between a quantum superposition and a classical sample however, is *quantum interference*. In particular, we may perform unitary transformations on this state to take it to some other state of the form

$$|\psi ^{}=\underset{j}{}\sqrt{q_j(\{p_i\})}|j.$$
(7)
Each of the new probabilities $`q_j`$ can in some sense be functions of the full set of $`p_i`$’s, and we may obtain “interference” between the probabilities that has no classical analogue. For example, a simple Walsh-Hadamard transformation on each of the the qubits produces the new distribution

$$q_j=\frac{1}{N}(\underset{i}{}(1)^{ij}\sqrt{p_i})^2$$
(8)
and this distribution is *not* log-concave.

(iii) *Estimating the magnitude of a fourier component.* A state of the form (1) can be fourier transformed efficiently shor . If we are interested in the magnitude of a partciular fourier component, then techniques exist for the process of *amplitude estimation*, and these proceed quadratically faster than their classical counterparts.

###### Acknowledgements.

This research was supported by the NSA & ARO under contract No. DAAG55-98-C-0040.