STABILITY OF MAGNETIZED DISKS AND IMPLICATIONS FOR PLANET FORMATION

During the gravitational collapse that forms star/disk systems, magnetic fields are dragged in from the interstellar medium (e.g., Galli et al. 2006; Shu et al. 2006). Additional fields can be generated by the central star. Consideration of mean field magnetohydrodynamics (MHD) in these disks shows that magnetic effects produce substantial departures from Keplerian rotation curves through both magnetic pressure and magnetic tension (Shu et al. 2007, hereafter S07). This sub-Keplerian rotation has important consequences for the ejection of disk winds (Shu et al. 2008) and the migration of extrasolar planets in protoplanetary disks (Adams et al. 2009). On the other hand, conservation of angular momentum implies that most of the material that eventually accretes onto the forming star initially lands on the disk (Cassen & Moosman 1981). As a result, disk surface densities can be high enough to support gravitational instability. In the limit of axisymmetric perturbations, the criterion for gravitational instability is determined by the value of the parameter Q_T

$$\begin{equation} Q_T \equiv \frac{a \kappa }{\pi G \Sigma }, \end{equation} \tag{ 1.1 }$$

where κ = ϖ⁻¹[∂(ϖ²Ω)²/∂ϖ]^1/2 is the epicyclic frequency, Ω is the angular rotation rate, a is the sound speed, and Σ is the surface density (Toomre 1964). In the presence of magnetic fields, however, the conditions required for gravitational instability are modified. The principal goal of this paper is to generalize the criterion of Equation (1.1) to include the effects of magnetic fields. More specifically, we derive a generalized stability parameter Q_M that characterizes magnetized disks.

We note that gravitational instability can play two important roles in circumstellar disks during the star formation process. If the instabilities grow into the nonlinear regime, they can produce secondary bodies within the disk, such as brown dwarfs and giant planets. If the growing perturbations saturate, the gravitational torques can lead to redistribution of angular momentum and disk accretion. Both processes require the onset of gravitational instability, which is determined by the parameter Q_M derived in this paper.

The properties and evolution of magnetized disks also depend on the dimensionless mass-to-flux ratio λ, defined by

$$\begin{equation} \lambda \equiv \frac{2\pi G^{1/2} \Sigma _0}{B_{z0}}. \end{equation} \tag{ 1.2 }$$

For example, gravitational collapse requires λ>1. As found by S07 and discussed herein, for realistic magnetic field strengths, this constraint inhibits the formation of giant planets by gravitational instability in circumstellar disks. In addition, the generalized stability parameter Q_M derived in this paper must be a function of λ.

As we show in this paper, the inclusion of magnetic fields leads to competing effects, some of which inhibit and some of which enhance gravitational instability and planet formation. However, as outlined above, magnetic fields will be present within circumstellar disks. As a result, in order to understand disk physics, one must include the effects of magnetic fields, and the goal of this paper is to provide an assessment of these effects.

This paper is organized as follows. We specify the equations of motion for magnetized disks in Section 2 and find their linearized counterparts in Section 3. This procedure leads to the dispersion relation for spiral density waves and the generalized stability parameter Q_M. In Section 4, we present numerical examples and apply the results to the observed protostellar source Ceph A HW2. The condition Q_M>1 is necessary for stability and implies a corresponding maximum disk mass, as shown in Section 5. The onset of instability and the derivation of Q_M can be determined by setting the resistivity η = 0; however, realistic disks have η ≠ 0 and Section 6 outlines the corresponding effects of magnetic diffusion. We then consider giant planet formation in Section 7. In addition to deriving modified constraints on planet formation via gravitational instability due to magnetic effects, we find that magnetic disks require an additional constraint: the need to remove magnetic flux places a lower bound on the electrical resistivity η. Finally, we conclude in Section 8 with a summary and discussion of our results.

This section specifies the equations of motion for this problem. We include the effects of a poloidal magnetic field dragged into the disk during the gravitational collapse of the natal cloud that produces a newly born star/disk system. This field threads vertically through the circumstellar disk and is pinched radially inward by viscous disk accretion. The accretion in these disks is believed to occur via the magneto-rotational instability (MRI; see, e.g., the review of Balbus & Hawley 1998). In fact, an empirical formulation of the MRI viscosity in thin disks has been obtained by S07 using mixing length arguments.

Consider the evolution of gas and magnetic field in a thin axisymmetric, viscously accreting disk of half-thickness z₀, surrounding a young star with mass M_⋆ at the origin of a cylindrical coordinate system (ϖ, z). We denote the surface density of the disk by Σ, the radial velocity of accretion in the plane by u, the azimuthal velocity about the z-axis by v, the component of the magnetic field threading vertically through the disk by B_z, and the radial component of the magnetic field just above the disk that responds to the radial accretion flow by B⁺_ϖ. The component of the Lorentz force per unit area in the plane of the disk can be written as

$$\begin{equation} {\bf f}_{\parallel }=\frac{{\bf B}^+_\parallel B_z}{2 \pi } -\nabla _\parallel \int \frac{B_z^2}{8 \pi } \, dz, \end{equation} \tag{ 2.1 }$$

where ${\bf B}^+_\parallel = B^+_\varpi \hat{\bf e}_\varpi +B^+_\varphi \hat{\bf e}_\varphi$ and $\nabla _\parallel = \hat{\bf e}_\varpi \partial /\partial \varpi +\hat{\bf e}_\varphi (1/\varpi)\partial /\partial \varphi$. The two terms in Equation (2.1) are associated with the effects of the magnetic tension and magnetic pressure (see Shu & Li 1997, who, however, adopt a different definition of magnetic pressure). In this paper, we adopt the approximation

$$\begin{equation} \int \frac{B_z^2}{8\pi }\, dz \approx \frac{z_0 B_z^2}{4 \pi }, \end{equation} \tag{ 2.2 }$$

where B_z on the right-hand side is evaluated in the midplane of the disk. In what follows, we neglect the toroidal component of the magnetic field threading the disk. In contrast, the stability of a disk with a purely azimuthal magnetic field was studied by Lynden-Bell (1966). With these specifications, the governing MHD equations in cylindrical coordinates include the equation of continuity,

$$\begin{equation} \frac{\partial \Sigma }{\partial t} + \frac{1}{\varpi } \frac{\partial }{\partial \varpi } (\varpi \Sigma u) + \frac{1}{\varpi ^2} \frac{\partial }{\partial \varphi } (\Sigma \varpi v) = 0; \end{equation} \tag{ 2.3 }$$

and the radial and azimuthal components of the equation of momentum,

$$\begin{eqnarray} &&\Sigma {\left(\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial \varpi } + \frac{v}{\varpi } \frac{\partial u}{\partial \varphi } - \frac{v^2}{\varpi }\right)} = -\frac{\partial \Pi }{\partial \varpi } - \Sigma {\left(\frac{\partial V}{\partial \varpi }\right)}_{z=0}\nonumber\\ &&\quad + \frac{B_z B_\varpi ^+}{2\pi } - \frac{z_0}{4\pi } \frac{\partial B_z^2}{\partial \varpi }, \end{eqnarray} \tag{ 2.4 }$$

$$\begin{eqnarray} &&\Sigma {\left(\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial \varpi } + \frac{v}{\varpi } \frac{\partial v}{\partial \varphi }\right)} + \frac{u v}{\varpi } = - \frac{1}{\varpi }\frac{\partial \Pi }{\partial \varphi }\nonumber\\ &&\quad - \frac{1}{\varpi }\Sigma {\left(\frac{\partial V}{\partial \varphi }\right)}_{z=0} - \frac{z_0}{4\pi \varpi } \frac{\partial }{\partial \varphi } \big(B_{\varpi }^{+2} + B_z^2\big), \end{eqnarray} \tag{ 2.5 }$$

where V is the gravitational potential and Π is the gas pressure integrated over the disk thickness. The vertical component of the induction equation takes the form

$$\begin{eqnarray} &&\frac{\partial B_z}{\partial t} + \frac{1}{\varpi } {\left[\frac{\partial }{\partial \varpi } (\varpi B_z u) + \frac{\partial }{\partial \varphi } (B_z v)\right]} = -\frac{1}{\varpi }\frac{\partial }{\partial \varpi }\nonumber\\ &&\quad\times \left[\eta \varpi \left(\frac{B_\varpi ^+}{z_0}-\frac{\partial B_z}{\partial \varpi }\right)\right] +\frac{1}{\varpi }\frac{\partial }{\partial \varphi } \left(\frac{\eta }{\varpi }\frac{\partial B_z}{\partial \varphi }\right) \,, \end{eqnarray} \tag{ 2.6 }$$

where η is the electric resistivity. Note that we ignore the viscous torque in the azimuthal component of the momentum Equation (2.5) because the viscous timescale is much longer than the gravitational instability timescale. Also, the azimuthal component of the magnetic field in the disk that arises from the stretching of the poloidal field by differential rotation averages to zero in the vertical integration.

The vacuum fields above and below the disk are treated using Green's function technique (see, e.g., Shu & Li 1997). Since B is current-free outside the disk, the magnetic field can be derived from a scalar potential

$$\begin{equation} {\bf B} = \boldsymbol{\nabla }\Psi \,. \end{equation} \tag{ 2.7 }$$

The condition ∇ ·  B = 0 then implies that Ψ satisfies Laplace's equation

$$\begin{equation} \nabla ^2 \Psi = 0\quad \mbox{ for $z \ne 0$} \,, \end{equation} \tag{ 2.8 }$$

subject to the boundary condition

$$\begin{equation} \left(\frac{\partial \Psi }{\partial z}\right)_{z=0} = B_z \,. \end{equation} \tag{ 2.9 }$$

The gravitational potential is given by the stellar and disk contributions, V = −GM_⋆/ϖ +  V_d, where V_d satisfies Poisson's equation for a thin disk,

$$\begin{equation} \nabla ^2 V_d = 4 \pi G \Sigma \,\delta (z) \,. \end{equation} \tag{ 2.10 }$$

For simplicity, in the following we assume an isothermal equation of state, Π = a²Σ, where a is the local sound speed.

This section presents a perturbation analysis of the equations of motion and derives a dispersion relation for spiral density waves in magnetized disks. First we expand to first order, and use the subscript 0 to denote zeroth order variables and 1 to denote first-order perturbations. It is understood that every variable is evaluated in the midplane (or just above the midplane). We thus look for solutions that have the Fourier decomposition

$$\begin{equation} F(\varpi, \varphi, t) = F_0 (\varpi) + F_1 (\varpi) e^{i(\omega t - m \varphi)} \,, \end{equation} \tag{ 3.1 }$$

where ω is a complex frequency and m is a positive integer. To zeroth order, we assume the disk to be in a state of axisymmetric radial equilibrium:

$$\begin{equation} \Omega ^2\varpi -\frac{a^2}{\Sigma _0}\frac{\partial \Sigma _0}{\partial \varpi } - \left(\frac{\partial V_0}{\partial \varpi }\right)_{z=0} + \frac{B_{z0}B_{\varpi 0}^+}{2\pi \Sigma _0} - \frac{z_0}{4\pi \Sigma _0}\frac{\partial B_{z0}^2}{\partial \varpi } = 0 \,, \end{equation} \tag{ 3.2 }$$

where Ω = v₀/ϖ.

The first-order equations become

$$\begin{equation} i(\omega - m \Omega)\frac{\Sigma _1}{\Sigma _0} + \frac{1}{\Sigma _0 \varpi } \frac{\partial }{\partial \varpi } (\varpi \Sigma _0 u_1) - i m \frac{v_1}{\varpi } = 0, \end{equation} \tag{ 3.3 }$$

$$\begin{eqnarray} &&\hspace{-34pt}i(\omega - m\Omega) u_1 - 2 \Omega v_1\nonumber\\ &&\hspace{-34pt}\quad = \frac{\Sigma _1}{\Sigma _0} {\left(- \frac{B_{z0} B_{\varpi 0}^+}{2\pi \Sigma _0} + \frac{a^2}{\Sigma _0}\frac{\partial \Sigma _0}{\partial \varpi } + \frac{z_0}{4\pi \Sigma _0} \frac{\partial B_{z0}^2}{\partial \varpi }\right)}\nonumber\\ &&\hspace{-34pt}\qquad -\frac{a^2}{\Sigma _0} \frac{\partial \Sigma _1}{\partial \varpi } - \frac{\partial V_1}{\partial \varpi } + \frac{B_{z0} B_{\varpi 1}^+}{2\pi \Sigma _0} + \frac{B_{z1} B_{\varpi 0}^+}{2\pi \Sigma _0}\nonumber\\ &&\hspace{-34pt}\qquad - \frac{z_0 B_{z1}}{2\pi \Sigma _0} \frac{\partial B_{z0}}{\partial \varpi } - \frac{z_0 B_{z0}}{2\pi \Sigma _0} \frac{\partial B_{z1}}{\partial \varpi }, \end{eqnarray} \tag{ 3.4 }$$

$$\begin{eqnarray} &&\hspace{-42pt} i(\omega - m \Omega) v_1 + u_1 \frac{\kappa ^2}{2\Omega } = im \frac{a^2}{\varpi \Sigma _0}\Sigma _1 + im \frac{V_1}{\varpi }\nonumber\\ &&\hspace{-42pt}\quad + im \frac{z_0}{\pi \varpi \Sigma _0} \big(B_{\varpi 0}^+ B_{\varpi 1}^+ + B_{z0} B_{z1}\big), \end{eqnarray} \tag{ 3.5 }$$

$$\begin{eqnarray} &&\hspace{-26pt} i(\omega - m\Omega) B_{z1} + \frac{1}{\varpi } {\left[\frac{\partial }{\partial \varpi } (\varpi B_{z0} u_1) -i m B_{z0} v_1\right]}\nonumber\\ &&\hspace{-26pt}\quad = \frac{i}{\varpi }\frac{\partial }{\partial \varpi }\left(\eta \frac{\varpi B_{z1}}{z_0}\right) -\eta m^2\frac{B_{z1}}{\varpi ^2}. \end{eqnarray} \tag{ 3.6 }$$

The last two terms in the parenthesis in Equation (3.4) correspond to the thermal and magnetic pressure of the zeroth order solution, and both of these terms are much smaller than the magnetic tension term. As a result, we can drop them in the following analysis. To proceed further, we invoke the WKB approximation and write the perturbed quantities as

$$\begin{equation} F_1 (\varpi) = {\hat{F}}_1(\varpi) e^{i k \varpi } \,, \end{equation} \tag{ 3.7 }$$

where k is the radial wavenumber. We make the further assumption that |k|ϖ ≫ 1, i.e., the spiral perturbations are tightly wrapped. As a result, we can ignore all derivatives of the amplitude, or division by ϖ, compared to derivatives of the phase. In other words, we may replace the radial derivative of a variable with multiplication by ik times its amplitude (for simplicity, we will omit the circumflex accent on the slowly varying amplitudes). As an additional simplification, we assume 1/|k|ϖ is of the same order as a/ϖΩ. This specification implies that the parameter |k|z₀ remains order unity. Next we note that the first-order perturbation in surface density and gravitational potential are related by

$$\begin{equation} \Sigma _1 = - \frac{|k| V_1}{2\pi G} \,, \end{equation} \tag{ 3.8 }$$

which follows from the asymptotic solution of the Poisson equation for V_d to the leading WKB order (see Shu 1992). Following a similar argument for Poisson's equation for the magnetic potential Ψ, we also obtain

$$\begin{equation} B_{z1} = - |k| \Psi _1^+ \,, \end{equation} \tag{ 3.9 }$$

where Ψ⁺₁ is the value of Ψ₁ just above the disk (Ψ₁ is an odd function of z: if B_z1>0, then Ψ₁ < 0 for z>0 and Ψ₁>0 for z < 0). From the definition (2.7) we have

$$\begin{equation} B_{\varpi 1}^+ \equiv \frac{\partial \Psi _1^+}{\partial \varpi } = -\frac{1}{|k|}\frac{\partial B_{z1}}{\partial \varpi } = -i \frac{k}{|k|}B_{z1} \,. \end{equation} \tag{ 3.10 }$$

Using the WKB approximation in conjunction with the above relationships, the first-order equations become

$$\begin{eqnarray} &i(\omega - m \Omega)\frac{|k|V_1}{2\pi G \Sigma _0} - ik u_1 = 0, \end{eqnarray} \tag{ 3.11 }$$

$$\begin{eqnarray} &&i(\omega - m\Omega) u_1 - 2 \Omega v_1 = \frac{|k|V_1}{2\pi G \Sigma _0} \left(\frac{B_{z0} B_{\varpi 0}^+}{2\pi \Sigma _0} + ik a^2 \right) - ikV_1\nonumber\\ &&\quad - i\frac{k}{|k|}(1+|k|z_0)\frac{B_{z0} B_{z 1}}{2\pi \Sigma _0} + \frac{B_{z1} B_{\varpi 0}^+}{2\pi \Sigma _0} - \frac{z_0 B_{z1}}{2\pi \Sigma _0} \frac{\partial B_{z0}}{\partial \varpi }, \end{eqnarray} \tag{ 3.12 }$$

$$\begin{eqnarray} && i(\omega - m \Omega) v_1 + u_1\frac{\kappa ^2}{2\Omega } = - im \frac{a^2|k|V_1}{2\pi G \Sigma _0\varpi } + im \frac{V_1}{\varpi }\nonumber\\ &&\quad +\, im \frac{z_0 B_{z1}}{\pi \Sigma _0\varpi } \left(B_{z0}-i\frac{k}{|k|}B_{\varpi 0}\right), \end{eqnarray} \tag{ 3.13 }$$

$$\begin{eqnarray} &i(\omega - m\Omega) B_{z1} + ik B_{z0} u_1 = -\eta |k| \frac{B_{z1}}{z_0}. \end{eqnarray} \tag{ 3.14 }$$

We analyze first the case of ideal MHD (we consider the effects of a non-zero resistivity in Section 6). Setting η = 0 and solving the first and last equations for u₁ and B_z1, we obtain

$$\begin{eqnarray} &u_1 = (\omega - m \Omega)\frac{|k|V_1}{2\pi G \Sigma _0k}, \end{eqnarray} \tag{ 3.15 }$$

$$\begin{eqnarray} &B_{z1} = \frac{B_{z0}}{\Sigma _0}\Sigma _1 = - \frac{|k|V_1}{\lambda G^{1/2}}\,, \end{eqnarray} \tag{ 3.16 }$$

where the mass-to-flux ratio is defined through Equation (1.2). After substituting these results, the remaining equations become

$$\begin{eqnarray} && i(\omega - m\Omega)^2\frac{|k|}{2\pi G \Sigma _0k} - 2 \Omega \frac{v_1}{V_1}\nonumber\\ &&\quad = ik a^2\frac{|k|}{2\pi G \Sigma _0} - ik + ik \frac{1}{\lambda ^2} (1 + |k|z_0), \end{eqnarray} \tag{ 3.17 }$$

$$\begin{eqnarray} &&\hspace{-1pc} i(\omega - m \Omega) \frac{v_1}{V_1} + (\omega - m \Omega)\frac{|k|}{2\pi G \Sigma _0k} \frac{\kappa ^2}{2\Omega }\nonumber\\ &&\quad\hspace{-1pc} = \frac{im}{\varpi } {\left[1-\frac{a^2|k|}{2\pi G \Sigma _0} - \frac{2|k|z_0}{\lambda ^2}\left(1-i\frac{kB_{\varpi 0}^+}{|k|B_{z 0}}\right)\right]} \,. \end{eqnarray} \tag{ 3.18 }$$

After further algebraic manipulation, the leading order dispersion relation takes the form

$$\begin{equation} (\omega - m \Omega)^2 = \kappa ^2 - 2\pi G \Sigma _0 |k|\epsilon + k^2 \Theta a^2 \,, \end{equation} \tag{ 3.19 }$$

where we define

$$\begin{equation} \Theta \equiv 1+\frac{B_{z0}^2z_0}{2\pi \Sigma _0 a^2}\quad \mbox{ and } \epsilon \equiv 1 - \frac{1}{\lambda ^2} \,. \end{equation} \tag{ 3.20 }$$

In the limit of vanishing magnetic field, this dispersion relation reduces to the familiar form for spiral density waves in a gaseous disk (e.g., Shu 1992). The magnetic field threading the disk modifies the standard dispersion relation for an unmagnetized disk by (1) replacing the sound speed a by the magnetosonic speed Θ^1/2a = (a² + v²_A0)^1/2, where v_A0 = B²_z0z₀/2πΣ₀ is the Alfvén speed at the disk midplane, and (2) diluting the effects of gravity by a factor if λ>1. If λ < 1, the right-hand side of Equation (3.19) is always positive and no instability occurs.

Although Equation (3.19) is valid in general for a thin disk with any radial distribution of the mass-to-flux ratio λ, it has the same form as the dispersion relation obtained by Shu & Li (1997) for the special case of a disk with spatially uniform λ (an "isopedic" disk). In particular, the marginal stability of isopedic disks with λ = 1 was demonstrated explicitly by Zweibel & Lovelace (1997). The magnetically modified Toomre Q_M parameter, which provides the boundary of stability for axisymmetric (m = 0) perturbations, is thus given by

$$\begin{equation} Q_M = \frac{\Theta ^{1/2}a\kappa }{\pi \epsilon G \Sigma _0} \,. \end{equation} \tag{ 3.21 }$$

Note that the definition of Θ in Equation (3.20) differs from that of Shu & Li (1997). For Q_M < 1, perturbations with wavenumber between $k_\pm =k_{\rm max}(1\pm\sqrt{\vphantom{A^R}\smash{\hbox{${1-Q_M^2}$}}})$ are unstable, with k_max = (/Θ)k_J being the wavenumber of maximum growth, and where k_J = πGΣ₀/a² is the Jeans wavenumber. Since /Θ < 1, the effect of the magnetic field is to increase the length scale of the gravitational instability with respect to the Jeans length scale.

Another important factor that determines Q_M in Equation (3.21) is the epicyclic frequency. Disks around young stars that have dragged in magnetic fields from the interstellar medium by gravitational collapse do not rotate at Keplerian speeds because magnetic tension modifies the force balance equation (see, e.g., Equation (18) of S07 when magnetic tension dominates over magnetic and gas pressure). S07 showed that in magnetized disks that are viscously accreting by the MRI, the rotation curve is sub-Keplerian by a constant fraction f. In their models, the sub-Keplerian parameter f is determined by their Equation (73) that states that the magnetic flux brought in by star formation is conserved and is left behind in the disk. Thus, for a given mass-to-flux ratio λ, the factor f depends on the stellar mass, M_⋆ (necessary to recover the flux brought in by star formation), the mass accretion rate, $\dot{M}_d$, and the system age, t_age. For λ ∼ 4, S07 obtained values of f in the range 0.39–0.95 for disks around low-mass and massive young stars (see their Table 2). For sub-Keplerian disks, the epicyclic frequency is given by

$$\begin{equation} \kappa = f \Omega _K = f (G M_\ast /\varpi ^3)^{1/2} \,. \end{equation} \tag{ 3.22 }$$

Therefore, the inclusion of magnetic fields produces competing effects on the instability parameter Q_M: the strong fields enforce sub-Keplerian flow, which reduces Q_M and leads to greater instability. On the other hand, both magnetic pressure and magnetic tension act to increase Q_M and lead to enhanced stability.

To evaluate the numerical values of the quantities derived in the previous section, we write $\Sigma _0=\mu \rm m_H N_{\rm H}$, $a=(3kT/2\mu \rm m_H)^{1/2}$, $\kappa =\Omega ={\cal G}$, where μ is the molecular weight, N_H is the hydrogen column density, T is the gas temperature, and ${\cal G}$ is the velocity gradient. The mass-to-flux ratio thus becomes

$$\begin{equation} \lambda =2.71 \, \mu \left(\frac{N_{\rm H}}{\mbox{$10^{24}$ cm$^2$}}\right) \left(\frac{B_z}{\mbox{mG}}\right)^{-1}, \end{equation} \tag{ 4.1 }$$

$$\begin{equation} \Theta =1+1.15\times 10^{-2} \, \left(\frac{B_z}{\mbox{mG}}\right)^2 \left(\frac{z_0}{\mbox{AU}}\right) \left(\frac{N_{\rm H}}{\mbox{$10^{24}$ cm$^2$}}\right)^{-1} \left(\frac{T}{\mbox{K}}\right)^{-1}, \end{equation} \tag{ 4.2 }$$

$$\begin{equation} \epsilon =1-1.36\times 10^{-1} \, \frac{1}{\mu ^2} \left(\frac{B_z}{\mbox{mG}}\right)^2 \left(\frac{N_{\rm H}}{\mbox{$10^{24}$ cm$^2$}}\right)^{-2}, \end{equation} \tag{ 4.3 }$$

$$\begin{eqnarray} Q_M &=& 2.12 \, \frac{\Theta ^{1/2}}{\epsilon \mu ^{3/2}} \left(\frac{{\cal G}}{\mbox{$10^{-2}$ km s$^{-1}$ AU$^{-1}$}}\right) \left(\frac{T}{\mbox{K}}\right)^{1/2}\nonumber\\ &&\times \left(\frac{N_{\rm H}}{\mbox{$10^{24}$ cm$^2$}}\right)^{-1}. \end{eqnarray} \tag{ 4.4 }$$

It is useful to consider an observed star/disk system where these results can be applied. The disk around the massive protostar Cepheus A HW2 is threaded by a large-scale magnetic field of strength B_z0 ≈ 23 mG at a radius R₀ ≈ 650 AU, inferred from methanol masers polarization (Vlemmings et al. 2010). The disk, observed in the continuum and in several molecular tracers (Patel et al. 2005; Jiménez-Serra et al. 2007, 2009) is seen almost edge-on, with the magnetic field roughly perpendicular to the disk and inclined with respect to the line of sight by an angle i ≈ 73° (Vlemmings et al. 2010). This inclination is in agreement with the high aspect ratio of the CH₃CN and SO₂ integrated emission, corresponding to inclinations of ∼68° and ∼79°, respectively (Patel et al. 2005; Jiménez-Serra et al. 2007). The disk shows a velocity gradient ${\cal G} \sin i \approx 6$ km s⁻¹ over 05 (Patel et al. 2005), corresponding to a radial range of ∼360 AU at a distance of 725 pc. Assuming quasi-Keplerian rotation, the velocity gradient at the radius where the magnetic field has been measured implies that ${\cal G} \sin i \approx 7\times 10^{-3}$ km s⁻¹ AU⁻¹. The inferred mass of the central star then becomes M_⋆ ≈ 15sin⁻²i M_☉. The gas temperature in the outer regions of the molecular disk is T ≈ 250 K (Jiménez-Serra et al. 2009), corresponding to an isothermal sound speed of a ≈ 1.8 μ^−1/2 km s⁻¹. Assuming that the disk is thermally supported, the scale height z₀ ≈ a/Ω at R₀ is then z₀ ≈ 260 μ^−1/2sin i AU. If the density in the methanol maser region has the typical value n_H ≈ 10⁹ cm⁻³, the hydrogen column density is N_H ≈ 7.8 × 10²⁴ μ^−1/2sin i cm⁻².

Assuming i ≈ 73° for consistency with the magnetic field strength determination, and μ = 2.33, we obtain λ ≈ 1.3 (in agreement with the estimate of Vlemmings et al. 2010), Θ ≈ 1.8, ≈ 0.45. With these values, the combined effects of magnetic pressure and tension increase the value of the Toomre Q_T parameter by a factor of Θ^1/2/ ≈ 3, reinforcing the stability of the disk from a marginally stable Q_T ≈ 1.4, to a safer Q_M ≈ 4.3.

Figure 1 shows the theoretical dimensionless mass-to-flux ratio λ (Equation (1.2)) and the magnetic parameter Q_M (Equation (3.21)) as a function of radius R. The disk is unstable to fragmentation at large radii where λ>1 and Q_M < 1. The sound speed a, λ, and the column density Σ (required to compute Q_M) were calculated using the profiles of magnetized disks in S07, for a system with a central star, M_⋆ = 16 M_☉, and mass accretion rate, $\dot{M}=10^{-4}$ M_☉ yr⁻¹, which correspond to the expected values in the circumstellar disk around the massive protostar protostar Cepheus A HW2. The curves correspond to different values of the sub-Keplerian constant f. The small circles indicate the values of λ and Q_M derived in the text at R = 650 AU and the best fit corresponds to f = 0.92.

Zoom In Zoom Out Reset image size

The numerical values discussed here are certainly subject to considerable uncertainties. We estimate that the calculated Q_M could vary in the range from a minimum of 1.3 to a maximum of 7.3, depending on the uncertainty on the disk's inclination, column density, magnetic field, and velocity gradient. While it is conceivable that this rather large range will be reduced by future observations, this example serves as a "proof of concept" to illustrate the stabilizing effect of magnetic fields in circumstellar disks, and their importance for the process of planet formation, addressed in the following sections.

In this section, we estimate the maximum disk mass that can remain stable against gravitational collapse in the presence of magnetic effects. We first manipulate the critical stability equation into the form

$$\begin{equation} 1 + \left(\frac{\pi G \Sigma _{\rm max}}{a \kappa } \right) \left(\frac{2 \beta }{\lambda ^2}\right) = \left(\frac{\pi G \Sigma _{\rm max}}{a \kappa } \right)^2 \left(1-\frac{1}{\lambda ^2}\right)^2 Q_M^2 \,, \end{equation} \tag{ 5.1 }$$

where we defined the disk scale height parameter β = κz₀/a.

If we set Q_M = 1 and consider both the flux-to-mass ratio λ and β to be a known functions of radius, we can solve Equation (5.1) for the critical (maximum) surface density profile

$$\begin{equation} \Sigma _{\rm max}=\left(\frac{a\Omega _K}{\pi G}\right)f{\cal M}\,, \end{equation} \tag{ 5.2 }$$

where

$$\begin{equation} {\cal M}=\frac{\lambda ^2}{(\lambda ^2-1)^2} \lbrace \beta +[\beta ^2+(\lambda ^2-1)^2]^{1/2}\rbrace. \end{equation} \tag{ 5.3 }$$

The first factor in Equation (5.2) represents the critical surface density in the absence of magnetic effects, the second factor represents the reduction of Keplerian rotation by magnetic tension, and the third represents the increase in the critical mass due to magnetic support. Note that ${\cal M}\rightarrow 1$ in the limit of unmagnetized disk (λ → ∞). In general, the dimensionless quantities f, λ, and β will be functions of the radial coordinate ϖ in the disk. However, one can choose constant representative values to get an idea of how disk stability depends on these quantities. In particular, for a thermally supported Keplerian disk, the scale height z₀ ≈ a/κ, and thus β ≈ 1.

Figure 2 illustrates the manner in which the magnetic field affects the stability of a circumstellar disk against gravitational perturbations, as expressed by Equation (5.2). Each curve in the λ–Q_T plane, where Q_T = afΩ_K/πGΣ₀, represents the condition Q_M = 1 for different values of the disk scale height parameter β. The thick curve corresponds to β = 1 and the thin curves correspond to β = 0 and $\beta =\sqrt{2}$, which is the maximum value allowed for a thermally supported disk (S07). For each value of β, the portion of the plane below the corresponding curve represents the region of parameter space for which disks are unstable in the presence of magnetic effects (Q_M < 1). Above the curve, Q_M is larger than unity, and the disk is stable. The effects of magnetic pressure and magnetic tension make a disk more stable compared to its unmagnetized counterpart, whereas the effects of sub-Keplerian rotation destabilize the disk. For f = 1, the line Q_T = 1 defines the boundary for stability in the absence of magnetic effects. The same curves also show the value of the inverse of the function ${\cal M}$ defined in Equation (5.3).

Zoom In Zoom Out Reset image size

Next we define the benchmark disk mass M_max integrating the critical surface density when there are not magnetic effects,

$$\begin{equation} M_{\rm max} \equiv \int _{R_\ast }^{R_d} 2 \varpi d\varpi \frac{a \Omega _K}{\pi G} \,. \end{equation} \tag{ 5.4 }$$

As shown by many authors (starting with Adams et al. 1988), if one uses the observed spectral energy distributions of T Tauri star/disk systems to specify the radial distribution of temperature, and hence the sound speed profile a(ϖ), the mass scales resulting from Equation (5.4) lie in the range M_max = 0.3–1 M_☉, i.e., masses comparable to those of the central stars. Note that this benchmark mass scale is calculated by assuming that Q_T = 1 throughout the disk, and thus represents an upper limit on the disk mass that can be stable. In practice, the Toomre parameter depends on radius, Q_T = Q_T(ϖ), so that much of the disk will have larger Q_T and hence smaller surface density than used in this exercise. Considerations of global stability show that the maximum disk mass is lower, e.g., the maximum disk mass that is stable to a class of m = 1 modes is given by M_d/(M_* +  M_d) = 3/4π (Shu et al. 1990). Nonetheless, this mass scale of Equation (5.4) provides an interesting benchmark. For the case of magnetized disks, the maximum mass is increased, as shown below.

With the inclusion of magnetic effects, the maximum disk mass that is stable to gravitational perturbations takes the form

$$\begin{eqnarray} M_{{\rm max}, M} &\equiv& \int _{R_\ast }^{R_d} 2 \varpi d\varpi \left(\frac{a \Omega _K}{\pi G} \right) \, f \, \left(\frac{\lambda }{\lambda ^2 - 1}\right)^2\nonumber\\ &&\times\, \lbrace \beta + [ \beta ^2 + (\lambda ^2 - 1)^2 ]^{1/2} \rbrace, \end{eqnarray} \tag{ 5.5 }$$

where, in general, f, λ, and β are functions of ϖ. In the simplest case, however, we can take these parameters to be constant, with representative values. In this case, the maximum allowed disk mass takes the form

$$\begin{eqnarray} &&\hspace{-15pt} M_{{\rm max}, M} = f {\cal M} M_{\rm max}\nonumber\\ &&\hspace{-15pt}\quad = M_{\rm max} \, f \, \left(\frac{\lambda }{\lambda ^2 - 1}\right)^2 \lbrace \beta + [ \beta ^2 + (\lambda ^2 - 1)^2 ]^{1/2} \rbrace \,. \end{eqnarray} \tag{ 5.6 }$$

This result differs from the field-free case by a factor F = $f {\cal M}$. For example, for the choice of parameters (f, λ, β) = (0.9, 1.3, 1), roughly corresponding to the observed values for the protostellar source Cepheus A HW2 considered in the previous section, the factor is F≈ 8. For somewhat larger mass-to-flux ratios λ = 2, the factor F ≈ 1.8. As a result, magnetic fields can produce a significant increase in the allowed masses of stable circumstellar disks. In the limit λ → ∞, we expect f → 1 so that F → 1.

It is important to keep in mind that the condition Q_M>1 is necessary for the disk to remain stable to axisymmetric perturbations only. An important result from spiral density wave theory is that disks which are stable to m = 0 perturbations can still be unstable to spiral perturbations with m ≠ 0. As a result, Equation (5.5) represents an upper limit to the maximum stable disk mass in these systems. Similarly, Equation (5.4) represents an upper limit for the maximum stable disk mass in the absence of magnetic fields. The actual maximum disk mass is smaller than these benchmark scales by a factor of ∼2 (e.g., Adams et al. 1989; Shu et al. 1990).

If we include the effects of magnetic diffusion, η ≠ 0, then the dispersion relation derived above is modified to take the form

$$\begin{eqnarray} && (\omega - m \Omega)^2 \left(\omega - m \Omega - i \eta \frac{|k|}{z_0} \right) = (\omega - m \Omega)\nonumber\\ &&\quad\times \left[ \kappa ^2 + k^2 a^2 {\left(1 + \frac{B_{z0}^2z_0}{2\pi \Sigma _0 a^2}\right)} - 2\pi G \Sigma _0 |k| \left(1-\frac{1}{\lambda ^2}\right) \right]\nonumber\\ &&\quad -i\eta \frac{|k|}{z_0}\left(\kappa ^2+k^2 a^2-2\pi G\Sigma _0|k|\right) -\eta k^2\frac{B_{z 0}B_{\varpi 0}^+}{2 \pi \Sigma _0 z_0}. \end{eqnarray} \tag{ 6.1 }$$

A basic analysis of this expression shows that the dispersion relation has no solutions for which the frequency ω is purely real. As a result, ω must be complex and can be written in the form

$$\begin{equation} \omega = \omega _0 + i \gamma \,, \end{equation} \tag{ 6.2 }$$

where both ω₀ and γ are real. In addition, since the parameter η is expected to be small, the imaginary part of the frequency is expected to be much smaller than the real part, |γ| ≪ ω₀. It is useful to define the functions

$$\begin{equation} A(k) \equiv \kappa ^2 + k^2 a^2 - 2\pi G \Sigma _0 |k| \,, \end{equation} \tag{ 6.3 }$$

and

$$\begin{eqnarray} B(k) &\equiv & \kappa ^2 + k^2 a^2 {\left(1 + \frac{B_{z0}^2z_0}{2\pi \Sigma _0 a^2}\right)} - 2\pi G \Sigma _0 |k| \left(1-\frac{1}{\lambda ^2}\right)\nonumber\\ &=& A(k)+\frac{2\pi G\Sigma _0}{\lambda ^2}|k|(1+|k|z_0) \,. \end{eqnarray} \tag{ 6.4 }$$

We also specialize to the case of axial symmetry (m = 0) so that the dispersion relation has the form

$$\begin{equation} \omega ^2 \left(\omega - i \eta \frac{|k|}{z_0} \right) = \omega B(k) - i \eta \frac{|k|}{z_0}A(k) - \eta k^2 b^2 \,, \end{equation} \tag{ 6.5 }$$

where we have also defined b² ≡ B_z0B⁺_ϖ0/(2πΣz₀). Solving for the real and imaginary parts of the dispersion relation, and working to leading order, we find

$$\begin{equation} \omega _0 \left(\omega _0^2 + 2 \gamma \eta \frac{|k|}{z_0} \right) = \omega _0 B(k) - \eta b^2 k^2 \,, \end{equation} \tag{ 6.6 }$$

and

$$\begin{equation} \omega _0^2 \left(3 \gamma - \eta \frac{|k|}{z_0} \right) = \gamma B(k) - \eta \frac{k}{z_0} A(k) \,. \end{equation} \tag{ 6.7 }$$

To leading order, ω²₀ = B(k), so the imaginary part of the frequency (the growth rate) takes the simple form

$$\begin{equation} \gamma =\eta \frac{|k|}{2 z_0}\left[1-\frac{A(k)}{B(k)}\right] \,. \end{equation} \tag{ 6.8 }$$

Since B(k) ⩾ A(k), with equality only in the limit |k| → 0 or λ → ∞, the growth rate is always positive, so that the solutions decay like exp(−γt).

This decay of the perturbation solutions is expected, in general terms, because the disk must spread for the case η ≠ 0. However, the specific form γ ∼ η|k|/z₀ is less obvious. Nonetheless, this result can be derived by solving the magnetic induction equation in the limit of an infinitesimally thin disk (and this calculation is carried out in the Appendix).

Gravitational instabilities in circumstellar disks can, in principle, lead to the formation of giant planets, or somewhat larger secondary bodies such as brown dwarfs (e.g., Boss 2001). In disk systems with significant magnetic support, however, the formation of secondaries is highly suppressed, as outlined in this section. The formation of secondary bodies requires both the onset of gravitational instability and sufficiently short cooling timescales (e.g., Gammie 2001). The required compromise between these two constraints is modified by sufficiently strong magnetic fields and is calculated in this section (compare with Rafikov 2005). Even when these two constraints are met, so that gravitational instability could in principle produce secondaries, the magnetic flux problem remains; in other words, the forming protoplanet must reduce its magnetic flux in order to contract to planetary sizes. This latter issue is also addressed below.

We first consider the coupled constraints of gravitational instability and sufficient cooling. In order for gravitational perturbations to grow, the stability parameter Q_M must be sufficiently small. Although the growth of axisymmetric instabilities requires Q_M ⩽ 1, spiral gravitational instabilities can grow in many star/disk systems under the weaker condition Q_M < q_* ≈ 2 (e.g., Adams et al. 1989), where the parameter q_* defines the required threshold. This constraint, a necessary condition for gravitational instability, takes the form

$$\begin{equation} Q_M = \frac{\Theta ^{1/2} a \kappa }{\pi \epsilon G \Sigma } < { q_\ast }\,, \end{equation} \tag{ 7.1 }$$

where the dimensionless parameters Θ and are defined above (Equation (3.20)).

The survival of gravitational instabilities requires that the cooling time t_cool is sufficiently short, and this condition takes the form Ω_Kt_cool < α_cool, where the value of α_cool ∼ 3. Note that this constraint uses the Keplerian rotation rate Ω_K; any departures from Keplerian can be incorporated into the value of α_cool. Following previous authors (e.g., Gammie 2001), we specify the form of the cooling time according to

$$\begin{equation} { t_{\rm cool} }= \frac{\Sigma a^2}{\gamma - 1} \cdot \frac{\tau + 1/\tau }{2 \sigma T^4} \,, \end{equation} \tag{ 7.2 }$$

where τ is the optical depth of the disk, γ is the adiabatic index, and T is the temperature. The cooling time constraint then takes the form

$$\begin{equation} \frac{\Omega _K \Sigma a^2}{\sigma T^4} \, < { p_\ast }\,, \end{equation} \tag{ 7.3 }$$

where we have defined the dimensionless parameter p_* ≡ 2α_cool(γ − 1)/(τ + 1/τ). Note that for typical values α_cool = 3, γ = 5/3, and τ = 1, the parameter p_* = 2.

By combining the constraints of Equations (7.1) and (7.3), one can show that the disk temperature—at the radial location where secondary formation could occur—obeys the ordering

$$\begin{equation} \left(\frac{\mu }{k} \right)^3 \left(\frac{\pi G \Sigma }{{\cal F} \kappa } \right)^6 > T^3 > \frac{k}{\mu } \cdot \frac{\Omega _K \Sigma }{{ p_\ast }\sigma } \,, \end{equation} \tag{ 7.4 }$$

where we have used the relation T = μa²/k, and

$$\begin{equation} {\cal F} \equiv \frac{1}{{ q_\ast }^2} \left(\frac{\lambda }{\lambda ^2 - 1} \right)^2 \lbrace \beta + [ \beta ^2 + { q_\ast }^2 (\lambda ^2 - 1)^2 ]^{1/2} \rbrace \,. \end{equation} \tag{ 7.5 }$$

By eliminating the temperature using the end points of Equation (7.4), one finds the following constraint on the surface density

$$\begin{eqnarray} &&\left(\frac{\mu }{k} \right)^4 \left(\frac{\pi G \Sigma }{{\cal F} \kappa } \right)^6 > \frac{\Omega _K \Sigma }{{ p_\ast }\sigma }\nonumber\\ &&\quad \Rightarrow \Sigma \, > \, \left(\frac{f {\cal F}}{\pi G} \right)^{6/5} \left(\frac{k}{\mu } \right)^{4/5} \frac{\Omega _K^{7/5}}{({ p_\ast }\sigma)^{1/5}} \,. \end{eqnarray} \tag{ 7.6 }$$

If we specialize to the case of a solar-mass star, with hydrogen gas (μ = 2m_P), and p_* = 2, the surface density constraint takes the form

$$\begin{equation} \Sigma > (3 \times 10^5 \, {\rm g} \, \, {\rm cm}^{-2}) \left(f {\cal F} \right)^{6/5} \, \left(\frac{\varpi }{1 \, {\rm AU}} \right)^{-21/10} \,. \end{equation} \tag{ 7.7 }$$

For comparison, the Minimum Mass Solar Nebula (MMSN) is expected to have a much lower surface density at ϖ = 1 AU, where Σ₁ ∼ 4500 g cm⁻² or smaller (e.g., Kuchner 2004). If the surface density profile is a negative power law with index p < 2.1, then gravitational instability can operate in the outer disk. For the parameter values listed here, planet formation via gravitational instability can only take place at radii exceeding the bound

$$\begin{equation} {\varpi } > [ 67 (f {\cal F})^{6/5} ]^{10/(21-10p)} \, {\rm AU} \approx 1100 \left(f {\cal F} \right)^{2} \, {\rm AU}, \end{equation} \tag{ 7.8 }$$

where the final equality specializes to the index p = 3/2, often used for the MMSN. Recall that the parameter f determines the degree to which the disk is sub-Keplerian and ${\cal F}$ is defined through Equation (7.5). Although the result is a somewhat complicated function of disk parameters (especially those that characterize the magnetic field), the required radius is always large. Planet formation via gravitational instability can only take place in the outer regions of large disks. We note that most disks surrounding low-mass stars have outer radii of order 100 AU and hence will not generally extend out to 1100 AU where giant planets could form. In addition, any planets that form at such large distances would have trouble migrating inward to the locations where (most of) the current sample of extrasolar planets resides (a < 10 AU). This channel of planet formation, via gravitational instability, is thus expected to have a limited impact on current observations

The MMSN, with its benchmark mass scale M_d ≈ 0.05 M_☉, typically has an outer radius of only 30 AU. For disks that have the same form for their surface density as the MMSN, but extend out to the radial scales of Equation (7.8), the corresponding mass is larger, about M_d ≈ 0.2 M_☉ $(f {\cal F})^{1/2}$. Note that the constraints on radial location are weaker if the index p for the surface density profile is smaller. However, the constraint on the mass scale is nearly the same.

Next we consider the magnetic flux problem in the context of forming secondary bodies. Even if the disk is heavy enough to become unstable, and the cooling time is short enough to allow contraction, magnetic flux must be transferred out of the region where the planet forms. The length scale ΔR that contains enough mass to form a Jovian planet is given by the integral over a disk annulus that contains the planet mass, i.e.,

$$\begin{equation} 2 \pi \int _{R_1}^{R_2} \Sigma \varpi d \varpi = m_P, \end{equation} \tag{ 7.9 }$$

where m_P ≈ 1 m_J = 2 × 10³⁰ g. For an MMSN disk model, where Σ = Σ₁(ϖ/R₁)^−3/2, where R₁ = 1 AU, this length scale is given by

$$\begin{eqnarray} \Delta R &=& R_2 - R_1 \approx \frac{m_P}{2 \pi \Sigma _1 R_1} \left(\frac{\varpi }{1 \, {\rm AU}} \right)^{1/2}\nonumber\\ &\approx& 0.31 \, {\rm AU} \, \left(\frac{\varpi }{1 \, {\rm AU}} \right)^{1/2} \, \approx 10 \, {\rm AU}, \end{eqnarray} \tag{ 7.10 }$$

where ϖ is the radial location of the forming planet and where we have chosen the nominal location ϖ = 1100 AU as a reference radius (see Equation (7.8)). Gravitational instability can potentially form fragments with Jovian mass and this length scale (10 AU), which is much smaller than Hill's radius, R_H = ϖ(m_P/3 M_*)^1/3 ≈ 76 AU at this location. As a result, the fragment can remain gravitationally bound.

Let us now consider flux freezing. The initial fragment is threaded with a magnetic field B₀, which can be written in terms of the mass-to-flux ratio λ through the relation B₀≈ 2πG^1/2Σ/λ. If we use the required values of Σ from the above analysis, and for λ of order unity, the initial field strength is about B₀ ≈ 10⁻⁴ G at the large radii where gravitational perturbations can grow. To form a giant planet, the fragment must contract to planetary size scales R_P∼10¹⁰ cm. Flux freezing implies that

$$\begin{equation} B_0 R_0^2 \approx B_P R_P^2 \,, \end{equation} \tag{ 7.11 }$$

where B_P is the magnetic field strength expected on the planetary surface. For the values R₀ = 10 AU and R_P = few R_J, the surface field strength would be B_P ≈ 2 × 10⁴ G. As a result, planet formation requires flux freezing to be compromised.

In these disks, the resistivity η ≠ 0 plays the role of a diffusion constant for the magnetic field. The time required for the resistivity to remove magnetic field from a region of size ℓ is thus the diffusion time t_diff and is given approximately by

$$\begin{equation} t_{\rm diff} \sim \ell ^2 / \eta \,. \end{equation} \tag{ 7.12 }$$

In this case, ℓ ∼ (z₀ΔR)^1/2 ∼ 45 AU, where the disk half-thickness is given by z₀ = Aϖ, and the aspect ratio is taken to be A = 0.1(ϖ/100 AU)^1/4. The value of the resistivity is η ≈ 2 × 10²⁰ cm² s⁻¹, required to dissipate enough magnetic flux to meet the constraints posed by measurements of paleomagnetism in meteorites (Shu et al. 2006). With these values, the diffusion time is only t_diff≈ 74 yr. To leading order, the magnetic field strength decreases exponentially with time, with timescale t_diff. In order to decrease the field strength by the required factor of 10⁴, the system needs 4log 10 ≈ 9.2 timescales, or about 680 yr. The expected values of the resistivity are thus high enough to allow the magnetic field to diffuse out of the protoplanet on a short timescale (680 yr), much shorter than the orbital timescale at the large radii where giant planets could form. Nonetheless, the required flux loss places a lower bound on the resistivity, i.e., the diffusion timescale t_diff from Equation (7.12) must be shorter than the cooling time. This constraint can be written in the form

$$\begin{equation} \eta > { \alpha _{\rm diff} }\Omega \ell ^2 \approx { \alpha _{\rm diff} }\Omega \left(\frac{m_P A}{2 \pi \Sigma } \right) \, \sim 2.5 \times 10^{18}\, {\rm cm^2\, s^{-1} }, \end{equation} \tag{ 7.13 }$$

where α_diff is a dimensionless parameter of order unity and the numerical estimate was obtained for the density of the MMSN at 1100 AU.

To conclude this section, the loss of magnetic flux is crucial to the planet formation process, but this flux can be removed with reasonable values of the resistivity η, which must obey the constraint of Equation (7.13). The bottleneck in the planet formation process is given by the coupled constraints of gravitational instability and short cooling timescales, where these constraints are modified significantly by the presence of magnetic fields. As a result, we expect giant planet formation via gravitational instability to be highly suppressed in these systems.

This paper has generalized the dispersion relationship for spiral density waves in circumstellar disks to include the leading order effects of magnetic fields for any radial distribution of the mass-to-flux ratio λ (Equation (3.19)). This procedure results in a generalized version of the stability parameter (denoted here as Q_M) for gravitational instabilities in magnetized disks, where Q_M is given by Equation (3.21). Magnetic fields produce competing effects regarding the stability of disks to gravitational perturbations: the increased pressure support and magnetic tension lead to greater stability; however, these same forces lead to sub-Keplerian rotation curves, which in turn lead to greater instability. The supporting terms generally dominate (see Figure 2), so that magnetic effects lead to an overall suppression of gravitational instabilities. In particular, there exists a maximum disk mass that is stable to gravitational perturbations. This maximum disk mass is larger for magnetized disks compared to those with B = 0 (see Section 5 and Equation (5.6)). These ideas can be tested through observations, as illustrated by the case of the disk surrounding the high-mass protostar Ceph A HW2 (Section 4).

The inclusion of magnetic fields leads to significant modifications to the prospects for the formation of giant planets through gravitational instability. In order to form giant planets, the disk must be unstable so that Q_M is small, and the cooling time must be short enough. Even in field-free disks, these coupled constraints limit the formation of planets to take place at large radii. The constraints become tighter in the presence of magnetic fields (see Equations (7.6) and (7.7)). Even in the event that gravitational instabilities do occur, and cooling timescales are short enough, the gravitationally bound fragments must lose large amounts of magnetic flux in order to contract to planetary size scales. This requirement leads to another constraint on the planet formation process, namely, a lower bound on the electrical resistivity (see Equation (7.13)).

The results of this paper indicate that giant planet formation via gravitational instability is difficult and hence should occur rarely. Nonetheless, rare is not the same as never: circumstellar disks that are sufficiently large (in radius) and massive could meet the constraints developed in this paper and support secondary formation. Some of the planetary candidates that have discovered through direct imaging (e.g., Marois et al. 2008; Kalas et al. 2008; Lagrange et al. 2009) could represent examples where this process has taken place.

This paper represents an important step toward understanding the effects of magnetic support in circumstellar disks. However, the generalized stability parameter Q_M derived here applies only to the onset of gravitational instabilities, specifically, linear perturbations with azimuthal wavenumber m = 0. Spiral modes (with m ≠ 0) can grow when axisymmetric perturbations are stable, so that this work should be generalized to include higher wavenumbers. In addition, the long-term fate of gravitational perturbations depends on the nonlinear evolution of these magnetized disks, and hence fully time-dependent MHD calculations should be carried out (e.g., Inutsuka et al. 2010). Since disks are expected to contain magnetic fields, these studies are vital in understanding both disk accretion and the possible formation of secondary bodies in star/disk systems.

We thank the Centro de Radioastronomía y Astrofísica, UNAM (Morelia, Mexico), the Physics Department at U. C. San Diego (La Jolla, CA), and Academia Sinica (Taiwan) for hosting and helping to facilitate this collaboration. We thank Zhi-Yun Li and Frank Shu for useful discussions. We also thank Mohsen Shadmehri for sharing unpublished calculations on the stability of cooling disks. S.L. is supported by grant CONACyT 48901 and grant PAPIIT 101310; M.J.C. is supported by grant NSC-95-2112-M-001-044; and F.C.A. is supported by NASA grant NNX07AP17G and NSF grant DMS-0806756.

Equation (6.8) shows that the exponential decay rate of a magnetic perturbation with radial wavenumber k in an infinitesimally thin disk is given by

$$\begin{equation} \gamma \sim \eta \frac{k}{z_0} \,. \end{equation} \tag{ A1 }$$

Thus, the characteristic diffusion length scale ℓ in an infinitesimally thin disk is neither the radial wavelength of the perturbation ∼1/k nor the disk thickness ∼z₀, but rather the geometric mean of the two, ℓ ∼ (z₀/k)^1/2. Since this result is curious and at variance with the usual assumption that ℓ ∼ z₀ in disks (e.g., see Parker 1979, p. 808, and Zeldovich et al. 1983), in the following we consider the expression of the decay rate derived by Bräuer & Rädler (1988) and Krause (1990) for a magnetic perturbation in a disk of finite thickness (without the effects of self-gravity, pressure, and rotation), and then we take the limit for z₀/ϖ → 0.

Consider a disk of finite thickness 2z₀ and infinite radius, and assume that η is constant inside the disk. In this case, the magnetic induction equation becomes

$$\begin{equation} \frac{\partial \bf B}{\partial t}=-\eta \nabla \times (\nabla \times {\bf B}) \qquad \mbox{for $|z|<z_0$} \,. \end{equation} \tag{ A2 }$$

Next we assume a vacuum field outside the disk,

$$\begin{equation} \nabla \times {\bf B}=0 \qquad \mbox{for $|z|>z_0$}, \end{equation} \tag{ A3 }$$

and impose the condition of continuity of B at z = ±z₀.

Any solenoidal field can be decomposed into a poloidal and toroidal components defined by the scalar functions P and T, respectively. Note that for an axisymmetric system P = ∂A_φ/ ∂ϖ, T = B_φ, and the terms "poloidal" and "toroidal" assume their usual meaning. The magnetic field thus can be written in the form

$$\begin{equation} {\bf B}=\nabla \times ({\bf \hat{e}}_z\times \nabla P) + {\bf \hat{e}}_z \times \nabla T. \end{equation} \tag{ A4 }$$

In terms of the defining scalars P and T, Equation (A2) separates into two pieces,

$$\begin{equation} \frac{\partial P}{\partial t}=\eta \nabla P \qquad {\rm and} \qquad \frac{\partial T}{\partial t}=\eta \nabla T \,. \end{equation} \tag{ A5 }$$

Equation (A3) thus becomes

$$\begin{equation} \nabla ^2 P=0 \qquad {\rm and} \qquad T=0 \,. \end{equation} \tag{ A6 }$$

The condition of continuity of B implies that P, T and ∂P/∂z must be continuous on z = ±z₀. The general solutions of Equations (A5) and (A6) that are regular at the origin can be expressed in terms of Bessel functions of the first kind with the arguments kϖ and k_zz, where k and k_z are the radial and vertical wavenumbers, respectively.

$${\fontsize{8.3}{10}\selectfont{\begin{eqnarray*} && P(\varpi,\varphi,z,t)\nonumber\\ && \quad =\left\lbrace \begin{array}{@{}ll} \int _0^\infty [a_1(k)\sin (k_z z)+b_1(k)\cos (k_z z)]J_m(k\varpi)e^{im\varphi -\gamma t}\, dk & \mbox{if $|z|<z_0$} \\[6pt] \nonumber \int _0^\infty c(k)e^{-k|z|} J_m(k\varpi)e^{im\varphi -\gamma t}\, dk & \mbox{if $|z|>z_0$}, \nonumber \end{array} \right. \end{eqnarray*}}}$$

and

$${\fontsize{8.3}{10}\selectfont{\begin{eqnarray*} &&T(\varpi,\varphi,z,t)\nonumber\\ &&\quad =\left\lbrace \begin{array}{@{}ll} \int _0^\infty [a_2(k)\sin (k_z z)+b_2(k)\cos (k_z z)]J_m(k\varpi)e^{im\varphi -\gamma t}\, dk & \mbox{if $|z|<z_0$} \\[6pt] \nonumber 0 & \mbox{if $|z|>z_0$}, \nonumber \end{array} \right. \end{eqnarray*}}}$$

where k and k_z are the radial and vertical wavenumbers, respectively, J_m(kϖ) are Bessel functions of the first kind and a₁, a₂, b₁, b₂, and c are functions determined by the initial field distribution and the conditions of continuity on z = ±z₀. The decay rate is

$$\begin{equation} \gamma =\eta \big(k_z^2+k^2\big), \end{equation} \tag{ A7 }$$

proportional to the square of the modulus of the vector wavenumber, ${\bf k}=k\hat{\bf e}_\varpi +k_z\hat{\bf e}$. The condition of continuity of ∂P/∂z at z = ±z₀ leads to the relation

$$\begin{equation} \big(k_z^2-k^2\big)\tan (2k_z z_0)=2 k_z k. \end{equation} \tag{ A8 }$$

The limit of an infinitesimally thin disk is recovered under the following ordering: z₀ ≪ k⁻¹_z ≪ k⁻¹ ≪ ϖ (the disk thickness is much smaller than the vertical wavelength, that in turn is much smaller than the radial wavelength). With these approximations, Equation (A8) then gives k ≈ k²_zz₀ and Equation (A7) becomes γ ≈ ηk²_z ≈ k/z₀, as found in our case.

Note that in a thin disk, the size of axially symmetric magnetic perturbations is limited vertically by the disk thickness, and, horizontally by the radial extent of the perturbation. The perturbative calculation shows that the diffusion rate is not determined by the smallest of these two scales (z₀, the disk thickness) but rather by the geometric mean of the two. In a disk of finite size, the relevant diffusion scale is determined, as expected, by the modulus of the vector sum of the vertical and horizontal wavenumbers; boundary conditions impose a relation between the two wavenumbers that, in the thin disk limit, gives the derived result.