Type I Planet Migration in a Magnetized Disk. II. Effect of Vertical Angular Momentum Transport
We study the effects of a large-scale, ordered magnetic field in protoplanetary disks on Type I planet migration using a linear perturbation analysis in the ideal-MHD limit. We focus on wind-driving disks, in which a magnetic torque (where and are the equilibrium vertical and azimuthal field components) induces vertical angular momentum transport. We derive the governing differential equation for the disk response and identify its resonances and turning points. For a disk containing a slightly subthermal, pure- field, the total 3D torque is close to its value in the 2D limit but remains lower than the hydrodynamic torque. In contrast with the 2D pure- field model considered by Terquem (2003), inward migration is not reduced in this case when the field amplitude decreases with radius. The presence of a subdominant component whose amplitude increases from zero at has little effect on the torque when acting alone, but in conjunction with a component it gives rise to a strong torque that speeds up the inward migration by a factor . This factor could, however, be reduced in a real disk by dissipation and magnetic diffusivity effects. Unlike all previously studied disk migration models, in the case the dominant contributions to the torque add with the same sign from the two sides of the planet. We attribute this behavior to a new mode of interaction wherein a planet moves inward by plugging into the disk’s underlying angular momentum transport mechanism.
Subject headings:accretion, accretion disks – MHD – planet–disk interactions – protoplanetary disks
Type I migration is the radial drift induced in a planet that is embedded in a protoplanetary disk by its linear gravitational interaction with the surrounding gas. This process, which is applicable to low-mass planets (typically up to Neptune’s mass), is of great relevance to theories of planet formation and evolution (see, e.g., Lubow & Ida 2010, Kley & Nelson 2012, Baruteau & Masset 2013, and Baruteau 2013 for recent reviews). Early treatments of this interaction focused on isothermal disks with a smooth density distribution and inferred rapid inward migration. The subsequent discovery of “hot Jupiters,” giant planets orbiting within stellar radii of their host stars, has exacerbated the conundrum elicited by this finding. These large, mostly gaseous, planets must have formed much farther out in their natal protoplanetary disks, where, according to the currently favored core-accretion model, they needed to assemble a large-enough ( a few times ) rocky core before a fast gas accretion phase (which created their envelopes) was triggered. However, the rapid inward drift predicted by the early models of Type I migration would prevent the postulated embryonic cores from lingering long enough at their formation sites for consistency with the observed distribution of hot Jupiters (e.g., Ida & Lin, 2008).
In the companion paper (Uribe et al. 2014, hereafter Paper I) we briefly review some of the more recent attempts to study Type I migration in the context of more realistic protoplanetary disk models, and we note a few of the promising results that have already been obtained. These include the effects of a disk magnetic field, which could be either a small-scale, disordered field associated with magnetohydrodynamic (MHD) turbulence (in particular, the turbulence induced by the magnetorotational instability) or a large-scale, ordered field that permeates the disk. These two field configurations likely play a central role in, respectively, the radial and vertical angular momentum transport in protoplanetary disks (e.g., Balbus, 2011; Königl & Salmeron, 2011). The effects of a small-scale, disordered field on planet migration include the stochastic motions induced in the smallest planets by MHD turbulence, the effective viscosity that prevents saturation of the corotation torque arising from the co-orbital region (which can be the dominant component of the torque on a planet undergoing Type I migration in a nonisothermal disk), and the magnitude of the mass threshold for opening a gap in the disk (i.e., for the onset of Type II migration), which depends on the value of the effective viscosity. A novel feature of the presence of a large-scale, ordered field is the appearance of new types of waves that propagate in the disk (specifically, slow- and fast-magnetosonic [SMS and FMS, respectively], as well as Alfvén), which supplant the sound waves that characterize a hydrodynamic (HD) disk. The MHD waves modify the way in which the disk responds to the gravitational perturbation induced by the planet, and, in turn, the way in which the disk acts back on the planet to cause it to migrate. The studies conducted so far have demonstrated that a purely azimuthal field in 2D can slow down, and even reverse, inward migration if the magnetic field pressure is not much smaller than the thermal pressure and if the field amplitude decreases with radius fast enough (Terquem, 2003; Fromang et al., 2005). It was also shown that, while a strong azimuthal field suppresses the corotation torque, a weaker one modifies it in a way that can potentially lead to outward migration (e.g., Guilet et al., 2013). In the case of a purely vertical field it was found by Muto et al. (2008) that the only effect in 2D is the replacement of the HD sound waves by the MHD FMS waves, with a resulting reduction in the magnitude of the torque, but that, in 3D, SMS and Alfvén waves also play a role (although these authors could not determine the net torque since their calculations were carried out in the shearing-sheet approximation).
As we observed in Paper I, a field configuration that is either purely azimuthal (and uniform with height) or purely vertical, as assumed in the aforementioned papers, is not representative of real protoplanetary disks. Such disks are likely permeated by open field lines that correspond to the interstellar field that originally threaded the parent molecular cloud cores and was dragged inward when the cores underwent gravitational collapse and formed the central stars and their associated circumstellar disks. The radially advected field lines assume a pinched (or hourglass) morphology and, in addition, are twisted by the differential rotation in the disk. In general, the field would thus have vertical, radial, and azimuthal components, which, near the disk surfaces, could have comparable magnitudes. This morphology is conducive to the formation of centrifugally driven outflows that can efficiently transport angular momentum away from the disk. This picture is supported by observations of protostellar disks as well as by theoretical considerations (see, e.g., Königl & Salmeron 2011 and references therein, as well as Bans & Königl 2012). Despite the strong azimuthal shear in the (typically Keplerian) protoplanetary disk, the magnetic field morphology in the giant-planet formation region can be expected to maintain an equilibrium structure on a timescale that greatly exceeds the local rotation period because the bulk of the gas at that distance (a few AU) is generally weakly ionized and therefore magnetically diffusive. In the simplest representation of such an equilibrium open-field configuration, the field has an even symmetry about the midplane ( in cylindrical coordinates ), with the equilibrium (subscript ‘0’) radial () and azimuthal () field components changing sign there (but with the vertical component remaining nonzero).
The aim of Paper I and this paper is to conduct a systematic investigation of the effects of a multi-component, ordered, large-scale magnetic field on Type I planet migration. The ultimate goal of this study is to model a planet that is embedded in a wind-driving disk, which, as the discussion above indicates, involves all three spatial components of the field and must, for self-consistency, be treated within the framework of nonideal MHD. Given the complexity of this problem, its full treatment is deferred to a future publication, and the present study concentrates on simpler field configurations that can be investigated using ideal MHD. We employ the complementary approaches of a linear perturbation analysis (the focus of this paper) and of numerical simulations (the main focus of Paper I). The linearization procedure is formulated in full generality in Section 2 of this paper, the numerical procedure employed in calculating the net torque on the planet is described in Section 3, and semianalytic solutions in 2 and 3 dimensions for three representative field configurations are presented in Section 4. The most general field configuration that can be treated self-consistently within the framework of ideal MHD involves a field that has both vertical and azimuthal components that can vary with but not with . This case is discussed extensively in Paper I (which, in Section I.2, also presents analytic results obtained from our perturbation analysis for this case).
2. Linear Analysis
2.1. Basic Equations
The dynamics of a magnetized protoplanetary disk is governed by the momentum conservation, mass conservation, and induction equations, given respectively by
where is the mass density, the gravitational potential (with and being the planet’s mass and orbital radius, respectively, and the gravitational constant), v the velocity field, B the magnetic field, and F the Lorentz force density
(with ). The field is required to satisfy the solenoidal condition .
The induction equation was written in its ideal-MHD form, which deserves a clarification. As already 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, since 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 (since any nonaxisymmetric perturbations would be washed out over timescales much larger than the inverse of the Keplerian angular velocity ). The ratio of these two timescales is quantified by the Elsasser number , where is the square of the Alfvén speed and 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 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). 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 centered on the star. The adopted field geometry and equilibrium disk conditions are described below. We also adopt an isothermal equation of state, , where 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 of the disk (assumed to be at the location of the planet) while imposing an even field symmetry. In this approach, the density is replaced by the surface density , although we will also be referring to the vertically averaged equilibrium density .
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 , , and by the orbital radius of the planet , frequencies by the planet’s orbital frequency , time by , and velocities by .
2.2. Equilibrium Conditions
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:
The geometrical thinness implies, by the solenoidal condition, that inside the disk is very nearly constant with height. Under the assumed even symmetry, the radial and azimuthal field components vanish at and increase in magnitude on 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 , where is a constant, and similarly for . Since the dominant contribution to in such disks is the shearing of by the differentially rotating gas, generally has the opposite sign of . We henceforth concentrate on the case where and are and
Using the above field morphology in Equation (4), we obtain for the vertically averaged equilibrium Lorentz force
where the triangular brackets denote a vertically averaged (between and ) 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 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 . 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 ,
where is the angular velocity and is the epicyclic frequency ().
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, , in the disk) is of the order of (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 was 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 term in the basic formulation
presented in this paper and assume an equilibrium velocity field of the form in the ensuing perturbation analysis. However, the semianalytic results that we obtain are derived in the limit , which should be an adequate approximation given that (which is of the order of ) is .
2.3. The Disk Response
To explore how our equilibrium disk model responds to the perturbing potential of the planet, (normalized by ), 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 3 directions, and Fourier transform in time as well as in space along the and directions by considering all Fourier quantities to be of the form , with the frequency taken to be . As in Terquem (2003), we denote all perturbed quantities with a prime. The planet’s potential itself can be expanded in a Fourier series, , where , and expressed in terms of the generalized Laplace coefficients (e.g., Ward, 1989; Korycansky & Pollack, 1993; Terquem, 2003):
where is the stellar mass (assumed to be ), the coefficient is equal to 0.5 for and to 1.0 for all the other (positive integer) values of ,
and where, for , , , and ; whereas for , , , and (with , a real number of magnitude , being the potential softening parameter; see Equation (I.20)).
where is the enthalpy perturbation, ( in normalized units), and where we suppressed the Fourier labels ( and ) on the primed quantities to simplify the presentation. To calculate the perturbed Lorentz force components in Equations (10) – (12), we use
and the following expression from Chandrasekhar & Fermi (1953),
which relates the perturbed magnetic field to the Lagrangian displacement . 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 is finite, but its use in the present analysis is justified by the fact that the results we derive are obtained in the limit .
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 (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 and 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
(e.g., Zhang & Lovelace, 2005). Thus, continuing from now on to suppress the Fourier labels on the perturbed quantities,
which can be integrated numerically (see Section 3).
2.4. Resonances and Turning Points
The solution of Equation (20) becomes singular at the locations where the leading-order coefficient, , 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 as a new dependent variable, which yields
When is real, Equation (21) has wave-like solutions if and evanescent solutions if . More generally, when is complex (with real and imaginary parts given by and , respectively), this equation always has wave-like solutions, but whether or not the waves propagate is determined by the magnitude and signs of and (see Appendix A for details). In our model, is complex for disks that contain a vertically variable equilibrium field component (nonzero or ), and real for vertically constant field configurations. The locations where the solution changes its character from wave propagation to evanescence (i.e., where ) 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.
2.5. Two- and Three-Dimensional Modes
With a nonzero , 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 from the equations. As discussed in Section 2.2, is an effectively small velocity, so in this initial study we drop all the 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
It can be seen from Equation (14) that includes a term of the form , which would result in the inclusion of a derivative of 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 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 and are nonzero; see Equation (3)), this approximation is consistent with our neglect of the terms.
With these two simplifications, the real part of the leading coefficient of the resulting second-order differential equation has two obvious roots, which yield two resonance conditions. The first one is
Here is the (rescaled) Alfvén velocity associated with the component of the field, , and .
If not for the vertical gradient in the 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 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 ), 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 and no vertical field, the ARs are defined by . Thus, the ARs in this case only exist for values of that satisfy , in stark contrast to the case of an azimuthal field with considered in Paper I (see Equations (I.11) and (I.12)), in which the ARs always exist. However, the ARs disappear, irrespective of the field geometry, in the strictly 2D limit (see Section I.2.1). The 2D cases of do, however, 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 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 , getting closer to the planet with increasing . For a pure- field, both of these resonances move further away from the planet as increases, diminishing the contributions of the resonances to the torque. For the case and , the MRs are located symmetrically with respect to the corotation radius, and their radial position does not vary with . This is similar to the behavior of the MR resonances for the 2D, pure- case (e.g., Terquem, 2003), except that in the case considered here the MRs only appear away from the midplane.
Figure 2 shows the outer () wave propagation regions and their relation to the locations of the turning points and resonances for the case of a pure vertical field and a finite value of . The magnetic field strength is parameterized by the value of , where is half the thermal-to-magnetic pressure ratio. (This ratio is generally near the midplane of a magnetically well-coupled, wind-driving disk; e.g., Königl & Salmeron 2011.) 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, , , , and , and then an additional one, , 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 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 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.
Muto et al. (2008) studied the pure-, 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 ) are evidently most closely related to the results we show in Figure 2. Their inner evanescence region, whose outer boundary (which they term LR-) corresponds to our ,
exists only for .
In the limit of the pure- 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. It is seen that 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 . 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 below) as well as by the 2D numerical simulations presented in Paper I (see Figure I.3). We further verified by numerical simulations that the decrease in the net torque with increasing applies also in 3D.
When both and are nonzero, the coefficients , , and , and hence (Equation (22)) are complex. In this case, we expect wave propagation when and as well as when (see Appendix A). Figure 4 shows the predicted regions of wave propagation (at a given ) for this field configuration for both 2D and 3D modes. The modes behave similarly to the 3D pure- case, which can be understood from the fact that, even at the disk surface (subscript ‘s’), the amplitude of remains a small fraction of ( in this example). However, when , there is an interior wave propagation region that starts very close to the planet (corresponding to a region of large ), which differs from the behavior of a pure- configuration in 2D. Although Figure 4 only shows the wave propagation regions for a given value of , we have verified that the area of the interior propagation region in both the 2D and 3D cases decreases with increasing (much as it does in the 3D pure- 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- disk, the torque would be enhanced (rather than reduced) in this case in comparison with the HD limit (see Section 4.4).
2.6. Effect of Vertical Gradient of
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 () 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 we noted in Section 2.5, when , 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 term in Equation (11) and solve for as a function of and its derivatives only. We then differentiate the resulting expression for and use it to evaluate the term we originally ignored so as to obtain an updated value for . 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 and also dropping the smaller-order terms. In this way we find the following approximate conditions:
where is defined analogously to 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 magnetic resonance condition for a pure- configuration. However, Equation (26) indicates that, when , there is an inherent asymmetry between the inner and outer Alfvén-like resonances. This equation is formally cubic in , and since is negative, it yields an additional resonance at . This extra resonance persists even when . This suggests that the presence of a vertical gradient of 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.
3. Numerical Scheme
To solve Equation (20), we employ the same methodology used by Korycansky & Pollack (1993) and Terquem (2003). To save computational time, we employ known recurrence formulas to evaluate the Laplace coefficients (Equation (9)) for a given . The relations we use are outlined in the appendix of Korycansky & Pollack (1993) and also in Brouwer & Clemence (1961) and Hagihara (1972).
The integration of Equation (20) is generally tricky due to the steepness of the planet potential and the presence of the various resonances. To achieve the required accuracy while using a standard initial-value integration routine, it is necessary to adopt small steps and require high precision; here we apply a fifth-order Runge-Kutta method with adaptive stepsize control (Press et al., 1992) and use 16-bit precision. To deal with any poles at , we follow Korycansky & Pollack (1993) and Terquem (2003) and displace the pole from the real axis by , where is a real number with , i.e., . This effectively also displaces the resonances from the real axis. We start the integration from the planet and integrate toward the inner and outer domain boundaries (specified below). We find that this methodology produces a smoother run than a “monotonic” integration from the inner boundary to the outer one.
At the boundaries of our integration domain, where , we neglect in comparison to , which is justified by the fact that, to leading order, . In this limit, the homogeneous version of Equation (20) becomes
where we retain only the leading-order term in each coefficient. Applying the WKB approximation by assuming that the solution of Equation (28) takes the form yields . For the case of a field with nonzero and , is given by
Although the terms that include are small in comparison with the other terms on the right-hand side of Equation (29), they can influence the results at small values of . Note also that, in the limit and , reduces to the dispersion relation for acoustic density waves in a 2D HD disk.
We use this WKB result as a boundary condition, , at the inner and outer edges of the integration domain. We still need initial values for and to start the integration, so, similarly to the standard shooting method and following Korycansky & Pollack (1993) and Terquem (2003), we start with random values and use an iterative approach to capture the true solution. For the sake of completeness, we provide a summary of this procedure below; additional details can be found in Korycansky & Pollack (1993).
Starting with random values implies that the solution obtained by the numerical scheme will in general contain a linear combination of solutions of the homogeneous equation (subscript ‘hs’) and of the particular solution of the inhomogeneous equation (subscript ‘ps’), i.e., . To attain convergence toward the true solution, we simultaneously integrate the two linearly independent homogeneous solutions, and . After one full integration, we can use the boundary conditions to constrain the constants and (specifically, our integrated solutions for and their derivatives yield a 2x2 system for and when plugged into at the disk boundaries). After solving for these constants, we perform one more integration, this time initializing as — this allows to converge toward . As in Terquem (2003), we find that good convergence is attained already after two iterations.
In contrast with the behavior of the 2D, purely azimuthal configuration considered by Terquem (2003), in which, for example, the location of the magnetic resonance was independent of , for the field configurations that we consider the locations of the resonances vary with the azimuthal mode number as well as with the vertical wavenumber . In our case, large vertical wavenumbers and small azimuthal mode numbers correspond to resonances (and, more generally, wave propagation regions) that lie far from the planet. Calculating the torque can then be computationally expensive, as it requires integration over a large region to capture the full effect of the resonances. To overcome this challenge, we limit ourselves to low values of (although we also calculate the torque for a few higher values of this product) and interpolate between some of the larger values of . For very high values of , we find it necessary to reduce the size of the integration domain (see Section 4.4). As was also the case in Korycansky & Pollack (1993), we find that for certain disk parameters we cannot calculate the contribution of some low azimuthal mode numbers (the differential equation becomes stiff and the needed step size to integrate is excessively small). In these cases we extrapolate to estimate the torque contribution from these low . This, in turn, introduces an error into the summed (over ) torque value, but we find that the qualitative behavior of the solution is not affected by this procedure. Note that, as in Korycansky & Pollack (1993) and Terquem (2003), we only attempt to calculate modes with . For , there is no effect on the torque since (see Equation (30)), whereas for , the results involve a contribution from the disk’s outer edge (see Shu et al., 1990), which we eliminate by enforcing an outgoing wave condition at the radial boundaries.
Solutions to Equation (20) are highly oscillatory, leading to strong cancellations of the torque contributions from adjacent zones, particularly in the outer regions of the disk. For this reason, the numerical results can depend on the size of the integration domain. In our treatment, we use and as the integration boundaries for . This domain is slightly larger than the one used by Korycansky & Pollack (1993) and Terquem (2003) since we need to ensure that the integration region is sufficiently large also for the cases, in which resonances and turning points lie farther away from the planet. For we adopt an even larger domain, and . We assume that the torque value for a given has effectively converged if the value does not differ by more than from the torque obtained using a slightly larger range. Following the practice in Korycansky & Pollack (1993) and Terquem (2003), we set the values of the parameters and to be and , respectively; however, our results are not sensitive to the specific choices of these values provided that they remain small (either one can be a factor of larger before any changes are detected in the results).
(see Korycansky & Pollack, 1993). To get the cumulative torque, we sum over (which in practice we do through an integral, assuming that a given value of characterizes a small interval of values). To obtain the total torque in the 3D cases, we approximate the additional sum over the vertical wave numbers by adding up the contributions of the specific values that we evaluate.
In this paper we use the convention, commonly adopted in linear analyses of the type we carry out (e.g., Korycansky & Pollack, 1993; Terquem, 2003), wherein one calculates the torque exerted by the planet on the disk. In this case, a positive torque implies inward migration whereas a negative torque indicates outward migration. This sign convention is the opposite of that commonly used in numerical studies, which concentrate on the effect that the disk has on the planet. The latter convention is also adopted in Paper I, where the focus is on numerical simulations. We caution the reader to be mindful of this difference when comparing the results of the two papers.
4.1. Hydrodynamic Limit
We ran our scheme with 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 for , and our calculated cumulative dimensionless torque of 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 vs. for the uniform-field case, shown in the top panel of Figure 6 of that paper).
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 . For the values used in this paper, we found the ratio to be closer to . We note, however, that an exact correspondence is not expected since 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 that we use (which in practice corresponds to the choice of , the effective disk half-thickness). Our chosen minimum value of () is based on applying our adopted value of () 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- 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 (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, in qualitative agreement with the Tanaka et al. (2002) result.
4.2. , Field Configuration
Although Muto et al. (2008) investigated the case of a pure- field, they formulated the problem in the shearing-sheet approximation and considered the 3D regime only for disks that are strongly magnetically dominated () and thus do not correspond to wind-driving systems. We therefore take a fresh look at this case.
Figure 5 shows the integrated torque as a function of for a uniform disk in the 2D limit and for two finite values of . 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 smaller than the HD value. As was 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 , is positive and generally larger than the corresponding value for . This is due to the appearance in this case 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 (), the net torque for each value of 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 increases. We thus expect the integrated torque to remain negative as grows further; however, its contribution to the total torque would be smaller than that of the mode because the resonances and associated wave-propagation regions move outward with increasing . We estimate the total torque in this case by adding up the net torques for the three values of shown in Figure 5; we find a value that is still smaller () 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- 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 , , and two values of (20 and 100) is shown in Table LABEL:tab1. For , we also show results for a field that is twice as strong (). 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 ( and , respectively); from the inner/outer to the inner/outer ( and , respectively), and from the inner/outer to the planet’s location ( and , 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.
The comparison between the and results for reveals that the contribution of the magnetic resonance region dominates, and its magnitude increases with , 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- 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 . These two trends combine to make the integrated torque much higher (with comparable contributions from the magnetic and Alfvén resonance regions) for . This behavior also explains why peaks at a comparatively low value of ( 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, similarly to Terquem (2003) for the pure- field, that the point-like contribution to 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 and converges toward a set sign at high .) 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 LABEL:tab1 also shows the breakdown of the integrated torque for a field that is twice as strong, i.e. , at . 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 ) indeed decreases in the AR region (as well as in the FMS wave propagation region beyond the outer Lindblad turning point) when is increased, other factors conspire in this case to increase the enthalpy amplitude in the MR region. We have verified this trend with numerical simulations, which show a more defined wake from SMS waves but a weaker Lindblad wake when 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 is decreased from 0.8 to 0.4. Since the torque at a given value of is not always representative of the cumulative torque (typically on account 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 case over the same range of and that was employed in Figure 5. We found that the total torque calculated in this way was again slightly smaller for the stronger field. We also carried out numerical simulations for the two values of considered in Table LABEL:tab1 and confirmed that the torque is reduced for the higher value of .
The numerical estimates given in the caption to Figure 5 indicate that, for , 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 (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 ) that the 3D torque decreases only slightly with increasing field strength for given values of and because of the opposing effects of the two sides of the planet, and that, in addition, the contributions to the total torque
largely cancel each other out.
4.3. , Field Configuration
We consider this field configuration both here and in Paper I even though it is not realistic in an attempt to isolate the effect of the term on planet migration in wind-driving disks. Figure 6 shows the net torque vs. for this configuration in the 2D limit and for the two vertical wavenumbers considered in Figure 5. As in the pure- case shown in the latter figure, we were unable to derive results for a few low- modes in the case, but their effect is likely not crucial in this case either. We find that the contributions from the two integrals that we present roughly cancel each other out as they did in the pure- case, and we thus approximate the total torque by the contribution of the 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- 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 .
We also already pointed out (in connection with Equation (25)) that, in 3D, ARs only exist in this case for . 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 , the smaller the critical value of , and since for low values of 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 () 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.
4.4. , Field Configuration
Having considered the effects of the and 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 (which equals ; see footnote 5) and (or, equivalently, of and ). The ratio of the azimuthal field amplitude at the disk surface to the vertical field component can be expressed as , 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 (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 gradient term to the uniform 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 for two different values of (while keeping the magnitude of as well as the value of fixed). It is apparent from the shape of the imaginary component of (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 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 , 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 modes for this field configuration (see next paragraph).
As the bottom panel of Figure 4 illustrates, the distribution of the wave propagation regions for the modes of this field configuration is qualitatively similar to that of a pure- field. In fact, given the smallness of the factor 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 LABEL:tab1, which provide the breakdown of the contributions to the integrated torque in this case for the same two values of as in the first two lines of this table, which list the corresponding contributions for the pure- field configuration (see Section 4.2). It is seen that, whereas the contribution of the MR regions increases with and that of the AR regions decreases with 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- 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 ( in this case): Those from the MRs are positive everywhere, whereas those from the ARs change from negative to positive as 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- 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 , 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 modes in this case is complicated by the fact that the values of are still growing even as approaches (see Figure 8). Integrating at high values of 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 , making however sure that, with each such reduction, the torque already calculated for a given value of did not change. The inset in Figure 8 shows the result of this calculation up to for the lowest vertical mode (). Since, for high values of , the contribution from the regions farthest away from the planet (the “Lindblad” and AR regions; and in Table LABEL:tab1) are negligible, the value of the integral is not meaningfully affected by shrinking the integration domain. At these very high values of , 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- behavior of the pure- field seen in Figure 4), and eventually (for ) the point-like contribution becomes negligible and converges to zero. The behavior of the other vertical modes is similar and they evidently also contribute a net positive torque (the case is shown in Figure 8 for ); however, since the contribution of the mode is by far the dominant one, we have not pursued the higher modes to full convergence.
The main goal of the present analysis has been to gain insights into Type I planet migration in wind-driving, magnetized disks. A key property of such disks is vertical angular momentum transport, which is brought about by the magnetic torque term . Given that the equilibrium structures of such systems cannot be studied self-consistently in the context of ideal MHD, we considered separately the behavior of pure- and pure- disks, and then attempted a linear perturbation analysis also for the combined field configuration.
In the case of a uniform vertical field component, we determined that the torque is reduced in comparison with the pure-HD torque in both 2D and 3D. This implies that planet migration would still be inward, but at a lower (by a factor of ) rate than in a disk not threaded by a large-scale field. For the case of a uniform azimuthal field in 2D, Terquem (2003) and Fromang et al. (2005) found that the torque is actually enhanced over the HD torque (resulting in faster inward migration). However, these authors also demonstrated that if decreases sufficiently rapidly with radius then inward migration can be strongly suppressed for () and even reversed (for ). The physical reason for this behavior is that an outward decrease in the field amplitude reduces the contribution of the outer MR region (where the torque exerted by the planet is positive) relative to that of the inner MR region (which has the opposite effect). An outward-decreasing field amplitude (corresponding to ; see Section 2.2) can be expected to arise naturally in the case of an open magnetic field that is dragged inward by the disk accretion flow, and it is thus of interest to inquire whether a similar reversal in the direction of planet migration is possible also for the astrophysically more relevant case of a vertical midplane field. We now show that the answer to this question is no.
In the 2D limit, a stronger vertical field leads to a decreased torque (see Section 4.2), in contrast with the pure