TYPE I PLANET MIGRATION IN A MAGNETIZED DISK. II. EFFECT OF VERTICAL ANGULAR MOMENTUM TRANSPORT

The dynamics of a magnetized protoplanetary disk is governed by the momentum conservation, mass conservation, and induction equations, given respectively by

$$\begin{eqnarray}&&\rho \left[ \frac{\partial {\boldsymbol{v}} }{\partial t}+({\boldsymbol{v}} \cdot \nabla ){\boldsymbol{v}} \right]=-\nabla P+{\boldsymbol{F}} -\rho \nabla \psi ,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho {\boldsymbol{v}} )=0,\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{B}} }{\partial t}+\nabla \times ({\boldsymbol{v}} \times {\boldsymbol{B}} )=0,\end{eqnarray} \tag{ 3 }$$

where ρ is the mass density, $\psi =G{{M}_{{\rm p}}}/|r-{{r}_{{\rm p}}}|$ the gravitational potential (with ${{M}_{{\rm p}}}$ and ${{r}_{{\rm p}}}$ being the planet's mass and orbital radius, respectively, and G the gravitational constant), ${\boldsymbol{v}}$ the velocity field, ${\boldsymbol{B}}$ the magnetic field, and ${\boldsymbol{F}}$ the Lorentz force density

$$\begin{eqnarray}&&{\boldsymbol{F}} =(\nabla \times {\boldsymbol{B}} )\times {\boldsymbol{B}} /\mu \end{eqnarray} \tag{ 4 }$$

(with $\mu \equiv 4\pi$). The field is required to satisfy the solenoidal condition $\nabla \cdot {\boldsymbol{B}} =0$.

The induction equation was written in its ideal-MHD form, which deserves a clarification. As noted in Section 1, the giant-planet formation zone in protoplanetary disks is expected to be weakly ionized and hence magnetically diffusive. This observation applies in particular to the disk midplane, where the planets reside, because this region is most strongly shielded from the main ionization sources of the disk (external cosmic rays, X-rays, and UV radiation). As a result, the ideal-MHD limit can be employed only under certain constraints. Since the focus of the current study is the effect of the magnetic field on the torque that the disk exerts on the planet (which arises from the nonaxisymmetric perturbations that the planet induces in the disk), it seems reasonable to require that the effective coupling time between the neutral disk component (which dominates the gravitational interaction with the planet) and the local magnetic field be shorter than the Keplerian rotation period (because any nonaxisymmetric perturbations would be washed out over timescales much larger than the inverse of the Keplerian angular velocity ${{{\Omega }}_{{\rm K}}}$). The ratio of these two timescales is quantified by the Elsasser number ${\Lambda }\equiv v_{{\rm A}}^{2}/{{{\Omega }}_{{\rm K}}}{{\eta }_{\bot }}$, where $v_{{\rm A}}^{2}={{B}^{2}}/\mu \rho$ is the square of the Alfvén speed and ${{\eta }_{\bot }}$ is the magnetic diffusivity perpendicular to the field (e.g., Königl et al. 2010). It can be expected that a magnetized disk would influence planet migration in an ideal-MHD–like way if ${\Lambda }\gg 1$ at the midplane. Perhaps not surprisingly, this is also the minimum coupling condition for sustaining angular momentum transport through either MRI-induced turbulence or a centrifugally driven wind (e.g., Königl & Salmeron 2011, pp. 283–352). It is conceivable that the diffusivity of a real disk could give rise to interesting effects that cannot be captured with the current formulation, but it nevertheless seems worthwhile to study the simpler ideal-MHD limit first.

Throughout this work we use an inertial, nonrotating, cylindrical coordinate system $(r,\varphi ,z)$ centered on the star. The adopted field geometry and equilibrium disk conditions are described below. We also adopt an isothermal equation of state, $P={{c}^{2}}\rho$, where c is the isothermal speed of sound. Although our focus is on the disk midplane, we allow for vertical wave propagation, which entails taking account of the 3D structure of the disk (e.g., Tanaka et al. 2002; Muto et al. 2008). We simplify the treatment by averaging over the effective thickness $2h$ of the disk (assumed to be $\ll r$ at the location of the planet) while imposing an even field symmetry. In this approach, the density ρ is replaced by the surface density ${\Sigma }$, although we will also be referring to the vertically averaged equilibrium density ${{\rho }_{{\rm av}}}={{{\Sigma }}_{0}}/2h$.

Following Muto et al. (2008), we further simplify our treatment by neglecting the effects of the vertical component of the stellar gravitational field on the disk structure.

In the remainder of this paper we normalize the spatial variables r, z, and h by the orbital radius of the planet ${{r}_{{\rm p}}}$, frequencies by the planet's orbital frequency ${{{\Omega }}_{{\rm p}}}$, time by ${\Omega }_{{\rm p}}^{-1}$, and velocities by ${{r}_{{\rm p}}}{{{\Omega }}_{{\rm p}}}$.

We adopt a simplified equilibrium magnetic field configuration intended to capture the basic geometry of an even symmetry open magnetic field that threads a geometrically thin disk:

$$\begin{eqnarray}&&{{{\boldsymbol{B}} }_{0}}=\{{{B}_{0r}}(z),{{B}_{0\varphi }}(z),{{B}_{0z}}(r)\}\ .\end{eqnarray} \tag{ 5 }$$

The geometrical thinness $h\ll 1$ implies, by the solenoidal condition, that B_0z inside the disk is very nearly constant with height. Under the assumed even symmetry, the radial and azimuthal field components vanish at z = 0 and increase in magnitude upon going away from the midplane (with opposite signs in the two hemispheres). Based on a simplified treatment of a wind-driving disk (e.g., Wardle & Königl 1993), we approximate ${{B}_{0\varphi }}(z)={{B}_{\varphi 0}}z$, where ${{B}_{\varphi 0}}$ is a constant, and similarly for ${{B}_{0r}}(z)$. Since the dominant contribution to ${{B}_{\varphi }}$ in such disks is the shearing of B_r by the differentially rotating gas, ${{B}_{\varphi 0}}$ generally has the opposite sign of ${{B}_{r0}}$. We henceforth concentrate on the case where B_0z and B_0r are $\gt 0$ and ${{B}_{0\varphi }}$ is <0.⁴ We adopt a power-law form ${{B}_{0z}}(r)={{B}_{z0}}{{r}^{{{q}_{z}}}}$ (with B_z0 constant) for the radial dependence of the vertical field component, which is motivated by self-similar models of wind-driving disks (e.g., Blandford & Payne 1982; Contopoulos & Lovelace 1994). The other equilibrium magnetic field components would in general also vary with $r$, but we neglect this dependence in this paper. (It is, however, incorporated into the functional form we adopt for ${{B}_{0\varphi }}$ in Paper I.)

Using the above field morphology in Equation (4), we obtain for the vertically averaged equilibrium Lorentz force

$$\begin{eqnarray}\begin{array}{rcl} \mu \langle {\boldsymbol{F}} \rangle & = & \left( {{B}_{r0}}{{B}_{z0}}{{r}^{{{q}_{z}}}}-{{q}_{z}}B_{z0}^{2}{{r}^{2{{q}_{z}}-1}}-B_{\varphi 0}^{2}\frac{{{h}^{2}}}{3r}, \right. \\ {} & {} & \left. {{B}_{\varphi 0}}{{B}_{z0}}{{r}^{{{q}_{z}}}}+{{B}_{\varphi 0}}{{B}_{r0}}\frac{{{h}^{2}}}{3r},0 \right), \\ \end{array}\end{eqnarray} \tag{ 6 }$$

where the triangular brackets $\langle \rangle$ denote a vertically averaged (between $-h$ and h) quantity. The first two terms in the radial component of this expression represent the magnetic tension and pressure forces, respectively, which act to reduce the inward pull of gravity. (The third radial term is generally much smaller due to the h² factor.) As we noted in Section 1, the strong shearing of a radial field component cannot be self-consistently incorporated into an equilibrium ideal-MHD model, so, in the applications given in this paper and in Paper I, we set ${{B}_{0r}}=0$. This eliminates the radial tension, which would typically dominate the radial Lorentz force in a real wind-driving disk. The first (and dominant) term in the azimuthal component of Equation (6) represents the braking torque acting on the disk. When substituted into the φ component of Equation (1) (which describes the conservation of angular momentum), we obtain an expression for the induced accretion velocity v_0r,

$$\begin{eqnarray}&&{{v}_{0r}}=\frac{2{\Omega }\langle {{F}_{\varphi }}\rangle }{{{\rho }_{{\rm av}}}{{\kappa }^{2}}},\end{eqnarray} \tag{ 7 }$$

where ${\Omega }$ is the angular velocity and κ is the epicyclic frequency (${{\kappa }^{2}}=\frac{2{\Omega }}{r}\frac{\partial ({{r}^{2}}{\Omega })}{\partial r}$).

Equation (7) highlights the fact that a disk threaded by a vertical magnetic field and possessing a nonzero vertical gradient of the azimuthal magnetic field would experience vertical angular momentum transport. This is a key property of disks that drive centrifugal outflows. Because of the strong vertical gradients that characterize geometrically thin disks, this transport could be very efficient. In fact, when such a wind is launched, the resulting accretion speed (assuming adequate magnetic coupling, ${\Lambda }\gg 1$, in the disk) is of the order of c (e.g., Salmeron et al. 2007), which is much higher than the speeds that arise from radial angular momentum transport by a small-scale, disordered field. As noted in Section 1, the field-line bending induced by the inflowing gas would give rise to a radial field component that, in turn, would lead to a rapidly growing azimuthal field component and a departure from equilibrium unless magnetic diffusivity were taken into account. Since we intend to include nonideal-MHD effects in a future work, we retain the v_0r term in the basic formulation presented in this paper, and assume an equilibrium velocity field of the form ${\boldsymbol{v}} =({{v}_{0r}},r{\Omega },0)$ in the ensuing perturbation analysis. However, the semianalytic results that we obtain are derived in the limit ${{v}_{0r}}=0$, which should be an adequate approximation given that $c/r{\Omega }$ (which is of the order of h) is $\ll 1$.⁵

To explore how our equilibrium disk model responds to the perturbing potential of the planet, $\psi ^{\prime}$ (normalized by $r_{{\rm p}}^{2}{\Omega }_{{\rm p}}^{2}$), we follow the method of Terquem (2003) and earlier analyses by considering small perturbations and linearizing the basic equations. We allow the planet to perturb the disk in all three directions, and Fourier transform in time as well as in space along the $\hat{\varphi }$ and $\hat{z}$ directions by considering all Fourier quantities to be of the form $X_{m,{{k}_{z}}}^{\prime }(r){{e}^{i(m\varphi +{{k}_{z}}z-\omega t)}}$, with the frequency taken to be $\omega =m\;{{{\Omega }}_{{\rm p}}}$. As in Terquem (2003), we denote all perturbed quantities with a prime. The planet's potential can be expanded in a Fourier series, ${{\psi }^{\prime }}=\sum _{m=0}^{\infty }\psi _{m}^{\prime }{\rm cos} m\phi$, where $\phi \equiv \varphi -{{{\Omega }}_{{\rm p}}}t$, and expressed in terms of the generalized Laplace coefficients $b_{1/2}^{m}(\alpha )$ (e.g., Ward 1989; Korycansky & Pollack 1993; Terquem 2003):

$$\begin{eqnarray}\psi _{m}^{\prime }=\left\{ \begin{array}{ccccccccccccccc} -\frac{{{M}_{{\rm p}}}}{{{M}_{*}}r}{{\delta }_{m}}b_{1/2}^{m}(\alpha ) & \quad \ \ \ \ r\gt 1 \\ -\frac{{{M}_{{\rm p}}}}{{{M}_{*}}}\ {{\delta }_{m}}b_{1/2}^{m}(\alpha ) & \quad \ \ \ \ r\lt 1, \\ \end{array} \right.\end{eqnarray} \tag{ 8 }$$

where M_* is the stellar mass (assumed to be $\gg {{M}_{{\rm p}}}$), the coefficient ${{\delta }_{m}}$ is equal to 0.5 for m = 0 and to 1.0 for all the other (positive integer) values of m,

$$\begin{eqnarray}&&b_{1/2}^{m}(r)=\frac{2}{\pi }\int _{0}^{\pi }\frac{{\rm cos} m\phi \ {\rm d}\phi }{{{({{q}^{2}}+{{p}^{2}}{{\alpha }^{2}}-2\alpha {\rm cos} \phi )}^{1/2}}},\end{eqnarray} \tag{ 9 }$$

and where, for $r\gt 1$, $\alpha =1/r$, q = 1, and ${{p}^{2}}=1+{{\epsilon }^{2}}$; whereas for $r\lt 1$, $\alpha =r$, p = 1, and ${{q}^{2}}=1+{{\epsilon }^{2}}$ (with , a real number of magnitude $\ll 1$, being the potential softening parameter; see Equation (I.20)).

Linearizing Equations (1) and (2) yields

$$\begin{eqnarray}&&\begin{array}{l} im\sigma v_{r}^{\prime }-2v_{\varphi }^{\prime }{\Omega }+\left( {{v}_{0r}}\frac{\partial v_{r}^{\prime }}{\partial r}+v_{r}^{\prime }\frac{\partial {{v}_{0r}}}{\partial r}+v_{z}^{\prime }\frac{\partial {{v}_{or}}}{\partial z} \right) \\ \quad =-\frac{\partial }{\partial r}({{\psi }^{\prime }}+{{W}^{\prime }})+\frac{2h}{{{{\Sigma }}_{0}}}\left( -\frac{{{W}^{\prime }}}{{{c}^{2}}}\langle {{F}_{r}}\rangle +\langle F_{r}^{\prime }\rangle \right), \\ \end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\begin{array}{l} im\sigma v_{\varphi }^{\prime }+\frac{{{\kappa }^{2}}}{2{\Omega }}v_{r}^{\prime }+{{v}_{0r}}\left( \frac{\partial v_{\varphi }^{\prime }}{\partial r}+\frac{v_{\varphi }^{\prime }}{r} \right) \\ \quad =-\frac{im}{r}\left( {{\psi }^{\prime }}+{{W}^{\prime }} \right)+\frac{2h}{{{{\Sigma }}_{0}}}\left( -\frac{{{W}^{\prime }}}{{{c}^{2}}}\langle {{F}_{\varphi }}\rangle +\langle F_{\varphi }^{\prime }\rangle \right), \\ \end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&im{{\sigma }_{0}}v_{z}^{\prime }+{{v}_{0r}}\frac{\partial v_{z}^{\prime }}{\partial r}=-i{{k}_{z}}{{W}^{\prime }}-i{{k}_{z}}{{\psi }^{\prime }}+2h\frac{\langle F_{z}^{\prime }\rangle }{{{{\Sigma }}_{0}}},\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}\begin{array}{rcl} \frac{im\sigma {{W}^{\prime }}}{{{c}^{2}}} & = & \frac{-1}{r{{{\Sigma }}_{0}}}\frac{\partial }{\partial r}\left( r{{{\Sigma }}_{0}}v_{r}^{\prime } \right) \\ {} & {} & -\frac{im}{r}v_{\varphi }^{\prime }-i{{k}_{z}}v_{z}^{\prime } \\ {} & {} & -\frac{1}{r{{{\Sigma }}_{0}}}\frac{\partial }{\partial r}\left( r{{{\Sigma }}^{\prime }}{{v}_{0r}} \right), \\ \end{array}\end{eqnarray} \tag{ 13 }$$

where ${{W}^{\prime }}={{{\Sigma }}^{\prime }}{{c}^{2}}/{{{\Sigma }}_{0}}$ is the enthalpy perturbation, $\sigma \equiv {\Omega }-{{{\Omega }}_{p}}$ ($=\;{\Omega }-1$ in normalized units), and where we suppressed the Fourier labels (m and k_z) on the primed quantities to simplify the presentation. To calculate the perturbed Lorentz force components in Equations (10)–(12), we use

$$\begin{eqnarray}&&\mu {{{\boldsymbol{F}} }^{\prime }}=(\nabla \times {{{\boldsymbol{B}} }^{\prime }})\times {{{\boldsymbol{B}} }_{0}}+\left( \nabla \times {{{\boldsymbol{B}} }_{0}} \right)\times {{{\boldsymbol{B}} }^{\prime }}\end{eqnarray} \tag{ 14 }$$

and the following expression from Chandrasekhar & Fermi (1953),

$$\begin{eqnarray}&&{{{\boldsymbol{B}} }^{\prime }}=\nabla \times ({\boldsymbol{ \xi }} \times {{{\boldsymbol{B}} }_{0}}),\end{eqnarray} \tag{ 15 }$$

which relates the perturbed magnetic field to the Lagrangian displacement ${\boldsymbol{ \xi }}$. This relation expresses the ideal-MHD condition of the magnetic field lines being "frozen" into the matter, and can be derived from the induction equation in the limit of a negligible equilibrium velocity. Equation (15) does not strictly hold when v_0r is finite, but its use in the present analysis is justified by the fact that the results we derive are obtained in the limit ${{v}_{0r}}=0$.

Note that the vertical averaging in the perturbation equations is done after all the spatial derivatives are evaluated. Thus, for example, the perturbed vertical momentum equation contains a term $\propto \partial {{B}_{0\varphi }}/\partial z$ (if it is nonzero) even though the equilibrium vertical magnetic force vanishes in our model (see Equation (6)). It is also worth noting that, in general, the perturbation equations do not contain terms involving the spatial derivative of the equilibrium density if the temperature is assumed to be a spatial constant (as is done in this paper, following, e.g., Korycansky & Pollack 1993; Terquem 2003). Consequently, the linearization results are not sensitive to our neglect of a vertical density gradient (induced by either tidal gravity or a vertical magnetic pressure gradient) in the equilibrium momentum equation.

The Lagrangian displacement is related to the Eulerian velocity perturbations through

$$\begin{eqnarray}&&{{{\boldsymbol{v}} }^{\prime }}=\frac{\partial {\boldsymbol{ \xi }} }{\partial t}+({{{\boldsymbol{v}} }_{0}}\cdot \nabla ){\boldsymbol{ \xi }} -(\xi \cdot \nabla ){{{\boldsymbol{v}} }_{0}}\end{eqnarray} \tag{ 16 }$$

(e.g., Zhang & Lovelace 2005). Thus, continuing from now on to suppress the Fourier labels on the perturbed quantities,

$$\begin{eqnarray}&&v_{r}^{\prime }=im\sigma {{\xi }_{r}}+{{v}_{0r}}\frac{\partial {{\xi }_{r}}}{\partial r}-\frac{\partial {{v}_{or}}}{\partial r}{{\xi }_{r}},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&v_{\varphi }^{\prime }=im\sigma {{\xi }_{\varphi }}-r{{\xi }_{r}}\frac{\partial {\Omega }}{\partial r}+{{v}_{0r}}\frac{\partial {{\xi }_{\varphi }}}{\partial r}-\frac{{{v}_{0r}}{{\xi }_{\varphi }}}{r},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&v_{z}^{\prime }=im\sigma {{\xi }_{z}}+{{v}_{0r}}\frac{\partial {{\xi }_{z}}}{\partial r}.\end{eqnarray} \tag{ 19 }$$

Substituting these expressions into Equations (10)–(13) and carrying out the necessary algebra, one obtains a second-order differential equation for ${{\xi }_{r}}$ of the form

$$\begin{eqnarray}&&\begin{array}{l} {{A}_{2}}(r)\frac{{{\partial }^{2}}{{\xi }_{r}}}{\partial {{r}^{2}}}+{{A}_{1}}(r)\frac{\partial {{\xi }_{r}}}{\partial r}+{{A}_{0}}(r){{\xi }_{r}} \\ \quad ={{S}_{1}}(r)\frac{\partial {{\psi }^{\prime }}}{\partial r}+{{S}_{0}}(r){{\psi }^{\prime }}, \\ \end{array}\end{eqnarray} \tag{ 20 }$$

which can be integrated numerically (see Section 3).⁶ The derived solution for ${{\xi }_{r}}$ can then be plugged back into the linearized equations of motion and continuity to obtain ${{W}^{\prime }}$ and to evaluate the torque exerted by the planet on the disk.

The solution of Equation (20) becomes singular at the locations where the leading-order coefficient, ${{A}_{2}}(r)$, goes to zero. These locations correspond to the resonances of the differential equation and are manifested by a divergence of the density perturbation. Their presence also signals the appearance of distinct regions of wave propagation in the disk. The torque exerted in the regions surrounding the resonances can potentially dominate the response of the disk and thus determine the rate and direction of planet migration. It is therefore important to calculate the resonant locations and examine their contribution to the integrated torque.

In general, the differential Equation (20) describes wave propagation in the disk, and we can rewrite a homogeneous version of it in terms of just the second- and zeroth-order terms. We do this, following Terquem (2003), by defining ${{\lambda }_{r}}\equiv {{\xi }_{r}}\ {\rm exp} \left( \frac{1}{2}\int \frac{{{A}_{1}}}{{{A}_{2}}}\partial r \right)$ as a new dependent variable, which yields

$$\begin{eqnarray}&&\frac{{{\partial }^{2}}{{\lambda }_{r}}}{\partial {{r}^{2}}}+K{{\lambda }_{r}}=0,\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&K=\frac{{{A}_{0}}}{{{A}_{2}}}-\frac{1}{4}{{\left( \frac{{{A}_{1}}}{{{A}_{2}}} \right)}^{2}}-\frac{1}{2}\frac{\partial }{\partial r}\left( \frac{{{A}_{1}}}{{{A}_{2}}} \right).\end{eqnarray} \tag{ 22 }$$

When K is real, Equation (21) has wave-like solutions if $K\gt 0$ and evanescent solutions if $K\lt 0$. More generally, when K is complex (with real and imaginary parts given by $\Re \{K\}$ and $\Im \{K\}$, respectively), this equation always has wave-like solutions, but whether or not the waves propagate is determined by the magnitudes and signs of $\Re \{K\}$ and $\Im \{K\}$ (see Appendix A for details). In our model, K is complex for disks that contain a vertically variable equilibrium field component (nonzero $\partial {{B}_{\varphi }}/\partial z$ or $\partial {{B}_{r}}/\partial z$), and real for vertically constant field configurations. The locations where the solution changes its character from wave propagation to evanescence (i.e., where $K=0$) are known as the turning points of the differential equation. Since it is the density wave propagation from the turning points that is responsible for the exchange of angular momentum with the planet (e.g., Goldreich & Tremaine 1979), finding the locations of the turning points makes it possible to identify the regions of the disk that exert a back torque on the planet. In the next subsection we evaluate the locations of the resonances and turning points of a simplified version of Equation (20) and show how their properties depend on the magnetic field configuration.

With a nonzero v_0r, Equations (10)–(13) all contain both zeroth-order and derivative terms of the perturbed variables, which means that it is impossible to eliminate all variables besides ${{\xi }_{r}}$ from the equations. As discussed in Section 2.2, v_0r is an effectively small velocity, so in this initial study we drop all terms in which it appears in the equations. Equation (15) then manifestly satisfies the linearized induction equation, and the (vertically averaged) perturbed field is given by

$$\begin{eqnarray}\begin{array}{rcl} \langle {\boldsymbol{B}} ^{\prime} \rangle & = & \left( -{{B}_{r0}}{{\xi }_{z}}+i{{k}_{z}}{{B}_{z0}}{{r}^{{{q}_{z}}}}{{\xi }_{r}},{{B}_{\varphi 0}}{{\xi }_{z}}+i{{k}_{z}}{{B}_{z0}}{{r}^{{{q}_{z}}}}{{\xi }_{\varphi }}, \right. \\ {} & {} & -\frac{im{{B}_{z0}}{{r}^{{{q}_{z}}}}{{\xi }_{\varphi }}}{r}-\frac{({{q}_{z}}+1){{B}_{z0}}{{r}^{{{q}_{z}}}}{{\xi }_{r}}}{r} \\ {} & {} & -\left. {{B}_{z0}}{{r}^{{{q}_{z}}}}\frac{\partial {{\xi }_{r}}}{\partial r} \right). \\ \end{array}\end{eqnarray} \tag{ 23 }$$

It can be seen from Equation (14) that ${{F}_{\varphi }}^{\prime}$ includes a term of the form ${{B}_{r0}}\partial {{\xi }_{\varphi }}/\partial r$, which would result in the inclusion of a derivative of ${{\xi }_{\varphi }}$ in Equation (11) and thereby make the system (10)–(13) difficult to solve. The leading coefficient in this case does not have a simple closed form that can be used to extract the resonances. Therefore, in the interest of maximizing our insight into the relevant resonances, we further simplify the treatment by letting B_0r be zero. Since the appearance of a radial field component is tied physically to the development of a vertically sheared radial velocity field in the disk (in particular, when B_z and $\partial {{v}_{or}}/\partial z$ are nonzero; see Equation (3)), this approximation is consistent with our neglect of the v_0r terms.

With these two simplifications, the real part of the leading coefficient ${{A}_{2}}(r)$ of the resulting second-order differential equation has two obvious roots, which yield two resonance conditions. The first is

$$\begin{eqnarray}&&\begin{array}{l} {{m}^{2}}{{\sigma }^{2}}\;= \\ \quad \frac{k_{z}^{2}{{c}^{2}}v_{{\rm A}z}^{2}+\frac{{{h}^{2}}}{3}\left( v_{{\rm A}z}^{2}\partial v_{{\rm A}\varphi }^{2}(k_{z}^{2}+\frac{{{m}^{2}}}{{{r}^{2}}})+\partial v_{{\rm A}\varphi }^{2}{{c}^{2}}\frac{{{m}^{2}}}{{{r}^{2}}} \right)}{{{c}^{2}}+v_{{\rm A}z}^{2}+\frac{{{h}^{2}}}{3}\partial v_{{\rm A}\varphi }^{2}}. \\ \end{array}\end{eqnarray} \tag{ 24 }$$

Here ${{v}_{{\rm A}z}}$ is the (rescaled) Alfvén velocity associated with the $z$ component of the field, $v_{{\rm A}z}^{2}=\frac{B_{z0}^{2}}{\mu {\Sigma }}$, and $\partial v_{{\rm A}\varphi }^{2}\equiv \frac{{{(\partial {{B}_{\varphi }}/\partial z)}^{2}}}{\mu {\Sigma }}=\frac{B_{\varphi 0}^{2}}{\mu {\Sigma }}$.⁷ Ignoring order-h² terms, the above equation indicates that one type of resonance occurs at the locations (one interior and one exterior to the planet's radius) where the corotating perturbation frequency $m\sigma$ is equal to the frequency of an SMS wave propagating along B_z. As in Paper I, we refer to these as magnetic resonances (MRs). The second resonance condition is

$$\begin{eqnarray}&&{{m}^{2}}{{\sigma }^{2}}=k_{z}^{2}v_{{\rm A}z}^{2}-\partial v_{{\rm A}\varphi }^{2}\left( 1-\frac{{{h}^{2}}{{m}^{2}}}{3{{r}^{2}}} \right).\end{eqnarray} \tag{ 25 }$$

If not for the vertical gradient in the $\varphi$ component, Equation (25) would indicate that there is another type of resonance at the locations where the perturbation frequency in the corotating frame matches that of a shear Alfvén wave propagating in the $\hat{z}$ direction. As in Paper I, we refer to these as Alfvén resonances (ARs). Equations (24) and (25) indicate that the MRs and ARs for the field configuration considered in this paper differ from those derived for the multicomponent field explored in Paper I (which is distinguished from the one discussed here by having a nonzero azimuthal component at the disk midplane). We find that the MRs depend strongly on the vertical field and very little on the azimuthal field gradient (which only appears in conjunction with factors of order h²), whereas for the configuration considered in Paper I, the MR locations are sensitive also to the azimuthal field component (see Equation (I.13)). Furthermore, in the limit of ${{k}_{z}}\ne 0$ and no vertical field, the ARs are defined by ${{m}^{2}}{{\sigma }^{2}}=-\partial v_{{\rm A}\varphi }^{2}(1-\frac{{{h}^{2}}\;{{m}^{2}}}{3{{r}^{2}}})$. Thus, the ARs in this case only exist for values of m that satisfy $\frac{{{h}^{2}}\;{{m}^{2}}}{3{{r}^{2}}}\gt 1$, in stark contrast to the case of an azimuthal field with ${{k}_{z}}\ne 0$ considered in Paper I (see Equations (I.11) and (I.12)), in which the ARs always exist. However, the ARs disappear, regardless of the field geometry, in the strictly 2D limit (see Section I.2.1). The 2D cases of $\partial {{B}_{\varphi }}/\partial z\ne 0$ do exhibit weak MRs (see Equation (24)), which is consistent with the results for a uniform azimuthal field (e.g., Terquem 2003), but they turn out to have a minimal effect on the net torque.

The dependence of the resonance locations on the field configuration and on the wavenumber k_z is illustrated in Figure 1. The MRs are shown by the black lines, and the ARs by the gray ones. For a pure vertical field, the ARs are always located outside the MRs (i.e., farther away from the corotation radius, where the disk's angular velocity is equal to that of the planet). Both resonances have locations that depend on the azimuthal wavenumber m, getting closer to the planet with increasing m. For a pure-B_z field, both of these resonances move further away from the planet as k_z increases, diminishing the contributions of the resonances to the torque. For the ${{B}_{\varphi 0}}\ne 0,\ {{B}_{z0}}=0$ case and ${{k}_{z}}\ne 0$, the MRs are located symmetrically with respect to the corotation radius, and their radial position does not vary with m. This is similar to the behavior of the MR resonances for the 2D, pure-${{B}_{\varphi }}$ case (e.g., Terquem 2003), except that in the case considered here the MRs only appear away from the midplane.⁸ In this case, an additional AR appears interior to the MR. As noted above, this AR only exists for wave numbers m above a critical value that depends on the disk thickness. For a thin disk with $h\ll 1$, this critical value is large, which makes the torque contribution of these ARs negligible.

Zoom In Zoom Out Reset image size

Figure 2 shows the outer ($r\gt 1$) wave propagation regions and their relation to the locations of the turning points and resonances for a pure vertical field and finite value of k_z. The magnetic field strength is parameterized by the value of ${{\beta }_{z}}\equiv {{c}^{2}}/v_{{\rm A}z}^{2}$, where β is half the thermal-to-magnetic pressure ratio. (This ratio is generally $\gtrsim 1$ near the midplane of a magnetically well-coupled, wind-driving disk; e.g., Königl & Salmeron 2011, pp. 283–352.) Each of the resonances (MRs and ARs) is associated with a pair of turning points that straddle it. We label the turning points according to the resonances they surround. From the planet's location outward, there are four such turning points—${{R}_{M-}}$, ${{R}_{M+}}$, ${{R}_{A-}}$, and ${{R}_{A+}}$—and they are followed by ${{R}_{L+}}$, the modified effective outer Lindblad resonance. Two regions of wave propagation, from which torque is exerted on the planet, surround each resonance and are bounded by the turning points. Waves also propagate away from the outer Lindblad turning point. The integrated torque exerted on the planet from the region $r\gt 1$ depends on the properties of these three regions, as well as on the point-like contributions from each of the resonances. The net torque, in turn, is determined by the balance between the opposing effects of the outer and inner disk regions, and by the contribution of the co-orbital region. Note that, as the value of the azimuthal mode number m increases, the resonances move closer to the planet but their associated regions of wave propagation become narrower. It is thus not obvious a priori which mode numbers dominate the torque. We address this issue on the basis of our numerical results in Section 4.2.

Zoom In Zoom Out Reset image size

Muto et al. (2008) studied the pure-B_z, ${{k}_{z}}\ne 0$ case using the shearing-sheet model and the WKB approximation. Similarly to what we find, they identified three regions of wave propagation, corresponding (in order of increasing distance from the planet) to the SMS, Alfvén, and FMS waves. Their Figures 1 and 2(c) (for which ${{\beta }_{z}}=0.9$) are evidently most closely related to the results in our Figure 2. Their inner evanescence region, whose outer boundary (which they term LR-) corresponds to our ${{R}_{M-}}$, exists only for $k_{z}^{2}v_{{\rm A}z}^{2}\gt 3\;{\Omega }_{{\rm p}}^{2}$.⁹ In our work we also find that this region is present only when $k_{z}^{2}v_{{\rm A}z}^{2}$ is sufficiently large, although the lower bound that we obtain with our more general treatment ($3.04\;{\Omega }_{{\rm p}}^{2}$) is slightly different from the analytic limit. Muto et al. (2008) identified the inner boundary of the intermediate evanescence region with the locus of the MRs rather than with that of the turning points ${{R}_{M+}}$ (which their derivation misses); we note in this connection that, even though we clearly distinguish between these two loci, in practice they are very close to each other in our example for all values of m. The inner boundary of the intermediate evanescence region, given by the turning points ${{R}_{A-}}$, appears to correspond to what these authors call vertical resonances.¹⁰ Their derivation, however, also misses the turning points ${{R}_{A+}}$, and thus they identify the inner boundary of the outer evanescence region with the locus of the Alfvén resonances. In our example, the Alfvén resonances are spatially well separated from the turning points ${{R}_{A+}}$ except at low values of m, where the WKB approximation (which requires $m\ll {{k}_{z}}r$) applies. This comparison illustrates the limitations of the WKB approximation in providing an accurate description of the disk–planet interaction. The outer boundary of the outer evanescence region is, however, correctly identified by this approximation with the locus of the modified effective Lindblad resonances.

In the limit ${{k}_{z}}=0$ of the pure-B_z case, both resonances vanish and there is only one set of turning points—the effective Lindblad resonances modified by the presence of the field, as described by Equation (37) of Muto et al. (2008). As noted above, our formalism does not involve the shearing-sheet approximation that was employed by these authors, and we are therefore able to evaluate the difference between the torques exerted by the inner and outer disk regions, and hence the net torque acting on the planet. Figure 3 shows the distances of the turning points on the two sides of the planet and their dependence on the field strength. The basic structure is similar to that in a purely HD disk, in that the outer effective Lindblad resonances, which exert a negative torque on the planet, lie closer to the planet than the inner turning points and can therefore be expected to dominate the interaction. This is true for all azimuthal mode numbers, although the asymmetry is stronger for low values of m. Under these circumstances, the planet is predicted to migrate inward. As the field strength increases, the turning points move farther away from the planet, which suggests that the magnitude of the net torque, and hence the rate of inward migration, are reduced. This heuristic inference is borne out by our numerical results (see Figure 5) as well as by the 2D numerical simulations presented in Paper I (see Figure I.3). Using numerical simulations, we further verified that the decrease in the net torque with increasing B_z also applies in 3D.

Zoom In Zoom Out Reset image size

When both ${{B}_{\varphi 0}}$ and B_z0 are nonzero, the coefficients A₀, A₁, and A₂, and hence K (Equation (22)) are complex. In this case, we expect wave propagation when $\Re \{K\}\gt 0$ and $\Re \{K\}\gg |\Im \{K\}|$, as well as when $|\Im \{K\}|\gg |\Re \{K\}|$ (see Appendix A). Figure 4 shows the predicted regions of wave propagation (at a given m) for this field configuration for both 2D and 3D modes. The ${{k}_{z}}\ne 0$ modes behave similarly to the 3D pure-B_z case, which can be understood from the fact that, even at the disk surface (subscript "s"), the amplitude of ${{B}_{\varphi }}$ remains a small fraction of B_z ($|{{B}_{\varphi {\rm s}}}/{{B}_{z0}}|=0.13$ in this example). However, when ${{k}_{z}}=0$, there is an interior wave propagation region that starts very close to the planet (corresponding to a region of large $|\Im \{K\}|$), which differs from the behavior of a pure-B_z configuration in 2D. Although Figure 4 only shows the wave propagation regions for a given value of m, we have verified that the area of the interior propagation region in both the 2D and 3D cases decreases with increasing m (much as it does in the 3D pure-B_z model shown in Figure 2). The presence of an inner wave-propagation region in 2D for this field configuration suggests that, in contrast with the pure-B_z disk, the torque would be enhanced (rather than reduced) in this case in comparison with the HD limit (see Section 4.4).

Zoom In Zoom Out Reset image size

As we already pointed out, a vertical gradient of the radial field component cannot be self-consistently incorporated into our ideal-MHD disk model. However, in view of the fact that the magnetic tension force ($\propto {{B}_{z}}\partial {{B}_{r}}/\partial z$) could play an important dynamical role in real wind-driving disks (e.g., Wardle & Königl 1993), it is of interest to at least look for clues to its possible effect on planet migration. As noted in Section 2.5, when ${{B}_{r0}}\ne 0$, one cannot obtain a differential equation in the form of Equation (20). We therefore attempt to apply an iterative approach to the procedure that yielded this equation. First, we ignore the ${{B}_{r0}}\partial {{\xi }_{\varphi }}/\partial r$ term in Equation (11) and solve for ${{\xi }_{\varphi }}$ as a function of ${{\xi }_{r}}$ and its derivatives only. We then differentiate the resulting expression for ${{\xi }_{\varphi }}$ and use it to evaluate the term we originally ignored in order to obtain an updated value for ${{\xi }_{\varphi }}$. We solve the rest of the system of equations using this updated value. In general, the leading-order coefficient of Equation (20) that is obtained in this way is rather complicated, and there is no simple form for its roots. However, one can obtain approximate resonance conditions by keeping the dominant terms in the limit ${{B}_{\varphi 0}}\to 0$ and also dropping the smaller-order $\propto h$ terms. In this way we find the following approximate conditions:

$$\begin{eqnarray}&&{{m}^{2}}{{\sigma }^{2}}-v_{{\rm A}z}^{2}k_{z}^{2}-2{{v}_{{\rm A}z}}\partial {{v}_{{\rm A}r}}\frac{\partial {\Omega }/\partial r}{\sigma }=0,\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\begin{array}{l} {{m}^{2}}{{\sigma }^{2}}(v_{{\rm A}z}^{2}+{{c}^{2}})-v_{{\rm A}z}^{2}{{c}^{2}}k_{z}^{2}+\partial v_{{\rm A}r}^{2}(v_{{\rm A}z}^{2}+{{c}^{2}}) \\ \quad -\partial {{v}_{{\rm A}r}}\frac{\partial {{v}_{{\rm A}z}}}{\partial r}(v_{{\rm A}z}^{2}+{{c}^{2}})=0, \\ \end{array}\end{eqnarray} \tag{ 27 }$$

where $\partial {{v}_{{\rm A}r}}$ is defined analogously to $\partial {{v}_{{\rm A}\varphi }}$ in Section 2.5.

In the absence of a radial field, Equation (26) reduces to the Alfvén resonance condition and Equation (27) to the MR condition for a pure-B_z configuration. However, Equation (26) indicates that, when $\partial {{v}_{{\rm A}r}}\ne 0$, there is an inherent asymmetry between the inner and outer Alfvén-like resonances. This equation is formally cubic in σ, and since $\partial {\Omega }/\partial r$ is negative, it yields an additional resonance at $r\gt 1$. This extra resonance persists even when ${{k}_{z}}=0$. This suggests that the presence of a vertical gradient of B_r could enhance the torque exerted by waves exterior to the planet's orbit, which induces inward migration. Although we cannot check on the validity of this prediction in the present study, we hope to pursue this case further in a future investigation.

We ran our scheme with $B=0$ to check against the HD results of Korycansky & Pollack (1993). In this limit, waves propagate away from the effective Lindblad resonances and there is a genuine resonance at the corotation radius. The total torque exerted by the planet is positive, representing the standard inward-migration result. We found that our results reproduce those of Korycansky & Pollack (1993) well. In particular, we obtained a good match to their plot (in Figure 2) of the real and imaginary parts of ${{W}^{\prime }}(r)$ for m = 10, and our calculated cumulative dimensionless torque of ∼1300 is in good agreement with the value reported in their Figure 14. Our results similarly match the HD calculations presented in Terquem (2003), and we further validated our numerical scheme by reproducing the azimuthal field results given in that paper (in particular, the plot of T_m versus m for the uniform field case, shown in the top panel of Figure 6 of that paper).¹¹

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also considered the 3D torque for this case. Tanaka et al. (2002) originally found that the 3D torque is weaker than the 2D torque by a factor of $\gtrsim 2$. For the ${{k}_{z}}h$ values used in this paper, we found the ratio to be closer to ∼1. We note, however, that an exact correspondence is not expected because Tanaka et al. (2002) used different disk parameters and included vertical stratification. Furthermore, the value of the 3D torque in our analysis also depends on the range of ${{k}_{z}}h$ that we use (which in practice corresponds to the choice of h, the effective disk half-thickness). Our chosen minimum value of ${{k}_{z}}h$ ($=1.56$) is based on applying our adopted value of ${{\beta }_{z}}$ ($=0.8$) to the condition for the existence of the inner evanescence region (or, equivalently, the condition for the existence of the innermost turning point) in the pure-B_z case (see Section 2.5). Regarding the upper limit of the adopted range, it is in practice limited by the computational difficulties encountered for large values of ${{k}_{z}}h$ (see Section 3). We were, however, able to explore the 3D torque in this case for sufficiently large values of this product (represented by ever thinner disks) to confirm that it is indeed weaker than the 2D HD torque, which is in qualitative agreement with the Tanaka et al. (2002) result.

Although Muto et al. (2008) investigated the case of a pure-B_z field, they formulated the problem in the shearing-sheet approximation and considered the 3D regime only for disks that are strongly magnetically dominated (${{\beta }_{z}}\leqslant {{10}^{-2}}$), and thus do not correspond to wind-driving systems. We therefore take a fresh look at this case.

Figure 5 shows the integrated torque T_m as a function of m for a uniform disk in the 2D limit and for two finite values of k_z. We also plot the HD result for comparison. The 2D integrated torque is invariably lower than in the unmagnetized case, resulting in a cumulative dimensionless torque for the chosen disk parameters that is a factor ∼2 smaller than the HD value. As already pointed out by Muto et al. (2008; see also Section 2.5 and Paper I), this reduction can be attributed to the outward shift of the effective Lindblad resonances in a magnetized disk, which increases with the field strength (see Figure 3). The behavior of the 3D modes is more complicated. For k_zh = 1.56, T_m is positive and generally larger than the corresponding value for k_z = 0. This is due to the appearance of the magnetic and Alfvén resonances and of their associated regions of wave propagation. However, for the larger vertical wavenumber considered in this figure (k_zh = 3.12), the net torque for each value of m is negative. This is a consequence of the fact—revealed by a close examination of Figure 1—that the inner MR lies slightly closer to the planet than the outer MR. This asymmetry becomes more pronounced, with a corresponding increase in its effect, as k_z increases. We thus expect the integrated torque to remain negative as k_zh grows further; however, its contribution to the total torque would be smaller than that of the k_zh = 3.12 mode because the resonances and associated wave-propagation regions move outward with increasing k_z. We estimate the total torque in this case by adding up the net torques for the three values of k_z shown in Figure 5, and find a value that is still smaller (∼0.7) than the 2D HD case shown in the same plot. This ballpark estimate is in good agreement with the total torque obtained in the numerical simulations presented in Paper I (see Figure I.10). In Section 3.4.1 of that work we point out that, for the adopted model parameters, our semianalytic estimates yield comparable values for the 3D and the 2D HD torques. Hence, our conclusion that the torque in a disk with a pure-B_z field is roughly half that of an unmagnetized disk applies in both 2D and 3D.

The breakdown of the contributions to the integrated torque from different regions in the disk for k_zh = 3.12, ${{\beta }_{z}}=0.8$, and two values of m (20 and 100) is shown in Table 1. For m = 20, we also show results for a field that is twice as strong (${{\beta }_{z}}=0.4$). The regions are chosen as follows (see Figure 2 for the meaning of the different labels): from the inner/outer edge of the disk to the inner/outer ${{R}_{L+}}$ (T_Li and T_Lo, respectively); from the inner/outer ${{R}_{L+}}$ to the inner/outer ${{R}_{A-}}$ (T_ARi and T_ARo, respectively), and from the inner/outer ${{R}_{A-}}$ to the planet's location (T_MRi and T_MRo, respectively). These sectors encompass the propagation regions of FMS, Alfvén, and SMS waves, respectively. Our ability to distinguish between the contributions of the regions on the inner and outer sides of the planet helps us gain insights that could not be obtained with a shearing-sheet model.

Table 1.  Net Torques for Different Field Configurations in 3D (${{k}_{z}}h=3.12$)

${\boldsymbol{B}}$	${{\beta }_{z}}$	${{\beta }_{\varphi }}$	m		TLi	TARi	TMRi	TMRo	TARo	TLo	${{T}_{{\rm integrated}}}$
Bz	0.8	⋯	20		−0.5	−16.5	−40.7	32.1	11.4	0.2	−14.0
Bz	0.8	⋯	100		−0.1	−2.3	−65.2	63.6	2.2	0.1	−1.7
Bz	0.4	⋯	20		>−0.01	−2.2	−54.5	43.0	0.1	<0.01	−13.6
${{B}_{z}}+\partial {{B}_{\varphi }}/\partial z$	0.8	0.04	20		−0.5	−20.4	15.3	15.3	−7.0	0.1	2.8
${{B}_{z}}+\partial {{B}_{\varphi }}/\partial z$	0.8	0.04	100		−0.2	0.2	60.9	59.7	0.8	0.1	121.5

Download table as:  ASCIITypeset image

A comparison between the m = 20 and m = 100 results for ${{\beta }_{z}}=0.8$ reveals that the contribution of the MR region dominates, and its magnitude increases with m, on each of the two (inner and outer) sides of the planet. However, the contributions from the two sides are closer in magnitude, resulting in a smaller net torque, for the higher-m case. On the other hand, the contribution from the Alfvén resonance region (from either side of the planet as well as the net one) decreases with m. These two trends combine to make the integrated torque much higher (with comparable contributions from the magnetic and Alfvén resonance regions) for m = 20. This behavior also explains why T_m peaks at a comparatively low value of m (see Figure 5).

The torque arising from the regions around a resonance is a combination of the point-like contribution from the resonance location and that from the associated wave-propagation zone. In the case of the MR, we find, similar to Terquem (2003) for the pure-${{B}_{\varphi }}$ field, that the point-like contribution to T_m typically has a smaller magnitude than the contribution from the surrounding region, and that these two contributions can have opposite signs. (For the most part, the direction of the point-like torque oscillates at low values of m and converges toward a set sign at high m.) The same behavior holds true for the AR, but in general the point-like contribution from the AR is much smaller in magnitude in comparison with both the contribution of the surrounding AR region and the point-like MR torque.

Table 1 also shows the breakdown of the integrated torque for a field that is twice as strong (i.e., ${{\beta }_{z}}=0.4$) at m = 20. Equations (24) and (25) indicate that, when the field becomes stronger, the magnetic and Alfvén resonances move further away from the planet, which has the effect of decreasing the magnitude of the planetary potential at their locations. We find, however, that while the strength of the disk response (as measured by the magnitude of $\Im \{{{W}^{\prime }}\}$) indeed decreases in the AR region (as well as in the FMS wave propagation region beyond the outer Lindblad turning point) when B_z0 is increased, other factors conspire to increase the enthalpy amplitude in the MR region. We verified this trend with numerical simulations, which show a more defined wake from SMS waves—but a weaker Lindblad wake—when ${{\beta }_{z}}$ is decreased. The decrease in torque from the Alfvén and FMS wave propagation regions even as the contribution from the MR region goes up results in a slight reduction in the total torque when ${{\beta }_{z}}$ is decreased from 0.8 to 0.4. Since the torque at a given value of m is not always representative of the cumulative torque (typically because of an ancillary effect such as the variation in the sign of the point-like contribution described in the preceding paragraph), we calculated the torque for the ${{\beta }_{z}}=0.4$ case over the same range of m and ${{k}_{z}}h$ that was employed in Figure 5. We found that the total torque calculated in this way was slightly smaller for the stronger field. We also carried out numerical simulations for the two values of ${{\beta }_{z}}$ considered in Table 1, and confirmed that the torque is reduced for the higher value of B_z.

The numerical estimates given in the caption to Figure 5 indicate that, for ${{\beta }_{z}}=0.8$, the total torque on the planet is roughly the same in 2D and 3D. If this inference is correct, the implied behavior would contrast with that of HD disks, in which the 2D torque dominates (e.g., Tanaka et al. 2002). Muto et al. (2008) inferred that the 3D torque exerted on a given side of the planet exceeds the corresponding 2D torque in the limit ${{\beta }_{z}}\ll 1$ (for which the MR contribution dominates), but this regime is not relevant to the formation region of giant planets. Furthermore, their analysis was carried out in the context of the shearing-sheet approximation, and therefore could not anticipate our findings (for ${{\beta }_{z}}\lesssim 1$) that the 3D torque decreases only slightly with increasing field strength for given values of k_zh and m because of the opposing effects of the two sides of the planet, and that, in addition, the ${{k}_{z}}\ne 0$ contributions to the total torque largely cancel each other out.¹²

We consider this field configuration here and in Paper I even though it is not realistic in an attempt to isolate the effect of the $\partial {{B}_{\varphi }}/\partial z$ term on planet migration in wind-driving disks. Figure 6 shows the net torque T_m versus m for this configuration in the 2D limit and for the two vertical wavenumbers considered in Figure 5. As in the pure-B_z case shown in the latter figure, we were unable to derive results for a few low-m modes in the k_zh = 3.12 case, but their effect is likely not crucial in this case either. We find that the contributions from the two ${{k}_{z}}\ne 0$ integrals we present roughly cancel each other out, as they did in the pure-B_z case, and we thus approximate the total torque by the contribution of the ${{k}_{z}}=0$ mode. In this way, we infer a value for the 2D torque that is nearly equal to that in the HD limit, which agrees with the result of the numerical simulations presented in Paper I (see Figure I.10).

As we discussed in Section 2.5, MRs appear above and below the midplane in the 3D solutions for this case, and the evidence for them—seen in the behavior of the perturbed enthalpy—persists into the 2D limit. In this limit, the MRs can be regarded as the "smeared out" (over the vertical extent of the disk) version of the 2D, pure-${{B}_{\varphi }}$ resonances found by Terquem (2003). These diluted MRs are, however, quite weak, and their contribution to the torque is very small. In addition, there are no ARs for this field configuration in the 2D limit. This is because the terms in the coefficients of Equation (20) that would have given rise to these resonances cancel out in the limit ${{k}_{z}},{{B}_{z}}\to 0$.¹³ These facts explain our finding that the magnitude of the 2D torque for this case is close to its HD value.

We also already pointed out (in connection with Equation (25)) that, in 3D, ARs only exist in this case for ${{m}^{2}}\gt 3{{r}^{2}}/{{h}^{2}}$. When this condition is satisfied, the ARs appear very close to the planet, even more so than the MRs (see Figure 1). The larger the value of h, the smaller the critical value of m, and since for low values of m the planet's potential couples well to the disk (see Equations (8) and (9)), the stronger could be the influence of the planet on the disk. In the example presented in Figure 6, we chose the disk half-thickness to be of the order of what would have been the density scale height ($c/{\Omega }$) had we explicitly included the vertical stellar gravity. This was sufficient to guarantee that, despite the fact that we assumed a uniform-density disk, the effect of the ARs on the integrated torque did not become artificially large.

Having considered the effects of the ${{B}_{z}}\ne 0$ and $\partial {{B}_{\varphi }}/\partial z\ne 0$ field configurations separately in the previous two subsections, we now investigate the outcome of combining them together. As we discussed in connection with Equation (7), the simultaneous presence of these two terms in the angular momentum conservation equation gives rise to vertical angular momentum transport in a magnetically threaded disk, and they are thus a key ingredient of a wind-driving disk model. We parametrize this configuration through the values of $\partial {{B}_{\varphi }}/\partial z$ (which equals ${{B}_{\varphi 0}}$; see footnote 3) and B_z0 (or, equivalently, of ${{\beta }_{z}}$ and ${{\beta }_{\varphi }}$). The ratio of the azimuthal field amplitude at the disk surface to the vertical field component can be expressed as $|{{B}_{\varphi {\rm s}}}/{{B}_{z0}}|=h\sqrt{{{\beta }_{z}}/{{\beta }_{\varphi }}}$, and is equal to 0.13 for our chosen parameters. This value is consistent with models of wind-driving disks that are magnetically active at the midplane, for which this ratio is typically inferred to be $\ll 1$ (e.g., Wardle & Königl 1993).

Not surprisingly, the behavior of the "combined" field configuration differs from that of its separate components. We noted in Section 2.5 that one distinctive feature of adding the ${{B}_{\varphi }}$ gradient term to the uniform B_z term in the 2D limit is the emergence of a wave propagation region in the immediate vicinity of the planet (see top panel of Figure 4). We illustrate the detailed structure of this region in Figure 7, which plots the radial dependence of the perturbed enthalpy ${{W}^{\prime }}$ for two different values of ${{\beta }_{\varphi }}$ (while keeping the magnitude of ${{\beta }_{z}}$ as well as the value of m fixed). It is apparent from the shape of the imaginary component of ${{W}^{\prime }}$ (which enters into the torque expression, Equation (30)) that these solutions exhibit wave-like behavior near the planet, and that this behavior becomes more pronounced as $|{{B}_{\varphi 0}}/{{B}_{z0}}|$ increases. This region can be expected to make a strong contribution to the torque, and we indeed find (see Figure 8) that the integrated torque in this case continues to grow with decreasing m, and that, correspondingly, the cumulative torque is significantly larger than the value for either of the 2D cases shown in Figures 5 and 6. Another noteworthy finding is that the contributions to the torque from this region have the same (positive) sign on both sides of the planet. This behavior, which was not previously encountered in planet–disk interactions, is also exhibited by the ${{k}_{z}}\ne 0$ modes for this field configuration.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the bottom panel of Figure 4 illustrates, the distribution of the wave propagation regions for the ${{k}_{z}}\ne 0$ modes of this field configuration is qualitatively similar to that of a pure-B_z field. In fact, given the smallness of the factor ${{h}^{2}}\partial v_{{\rm A}\varphi }^{2}/v_{{\rm A}z}^{2}=B_{\varphi {\rm s}}^{2}/B_{z0}^{2}\approx 0.02$ in Equations (24) and (25), the locations of the MRs and of the ARs are essentially the same in these two cases. However, the disk response is very different. This is demonstrated by the last two lines in Table 1, which provide the breakdown of the contributions to the integrated torque in this case for the same two values of m as in the first two lines of this table, which list the corresponding contributions for the pure-B_z field configuration (see Section 4.2). It is seen that, whereas the contribution of the MR regions increases with m and that of the AR regions decreases with m in both cases, other major aspects of the disk behavior have changed. In particular, while the regions interior and exterior to the planet contribute with opposite signs in the pure-B_z case (as they also do in an HD disk), in this case the contributions from the MR and AR regions have the same sign on both sides of the planet for sufficiently high values of m ($\gtrsim 20$ in this case): those from the MRs are positive everywhere, whereas those from the ARs change from negative to positive as m increases (and we verified that in each case the point-like contribution has the same sign as that from the surroundings). Thus, in contrast with the behavior of the pure-B_z field, where the contributions from the two sides of the planet offset each other and lead to comparatively small integrated torques, in this case the growth of the MR torque with increasing m, and the fact that the contributions from the two sides of the planet add up, can result in very large integrated torques.

Calculating the cumulative torque for the ${{k}_{z}}\ne 0$ modes in this case is complicated by the fact that the values of T_m(m) are still growing even for m near 100 (see Figure 8). Integrating $\partial {{T}_{m}}/\partial A$ at high values of m is computationally challenging because of the high oscillation frequency, which necessitates taking ever smaller step sizes. We have overcome this difficulty by reducing the integration-domain boundaries for high values of m, making sure that, with each such reduction, the torque already calculated for a given value of m did not change. The inset in Figure 8 shows the result of this calculation up to m = 1500 for the lowest vertical mode (k_zh = 1.56). Since, for high values of m, the contribution from the regions farthest away from the planet (the "Lindblad" and AR regions; ${{T}_{{\rm Li}/{\rm o}}}$ and ${{T}_{{\rm ARi}/{\rm o}}}$ in Table 1) are negligible, the value of the integral is not meaningfully affected by shrinking the integration domain. At these very high values of m, the contribution from the MRs—which is added with the same sign from both sides of the planet—dominates. However, as the mode number continues to increase, the physical width of the resonance becomes progressively smaller (corresponding to a decreasing width of the wave propagation region, akin to the high-m behavior of the pure-B_z field seen in Figure 4), and eventually (for $m\gtrsim 700$) the point-like contribution becomes negligible and T_m converges to zero. The behavior of the other vertical modes is similar and they evidently also contribute a net positive torque (the case k_zh = 3.12 is shown in Figure 8 for $m\leqslant 100$); however, since the contribution of the k_zh = 1.56 mode is by far the dominant one, we have not pursued the higher k_z modes to full convergence.