ON HYDROMAGNETIC STRESSES IN ACCRETION DISK BOUNDARY LAYERS

Basic arguments suggest that the angular frequency Ω(r) of an accretion disk surrounding a weakly magnetized star must attain a maximum value, Ω_max ≡ Ω(r_b), and decrease inward (or at least remain constant; see Medvedev 2004) to match the angular frequency of the star at the stellar radius (Ω_⋆(r_⋆); see, e.g., Frank et al. 2002; Hartmann 2009; Armitage 2010). The inner disk region, where r < r_b and dΩ/dr ⩾ 0, is referred to as the accretion disk boundary layer. Standard accretion disk theory (Shakura & Sunyaev 1973) predicts that half of the energy released in the accretion process takes place in this region, estimated to be a fraction of the stellar radius. The spectrum of the radiated energy depends on the detailed properties of this layer (Narayan & Popham 1993; Popham & Narayan 1995; Popham et al. 1996; Popham & Sunyaev 2001). Thus, understanding the various processes that determine its properties (see, e.g., Piro & Bildsten 2004; Balsara et al. 2009; Inogamov & Sunyaev 2010) is of fundamental importance.

Most detailed calculations for determining the structure of the boundary layer rely on effective models for turbulent angular momentum transport. These models are usually built as a turbulent version of the Newtonian viscous stress between fluid layers in a differentially rotating laminar flow (Landau & Lifshitz 1959), and thus assume a linear relationship between the stress and the angular frequency gradient (Lynden-Bell & Pringle 1974). This assumption, however, seems at odds with the properties of magnetohydrodynamic (MHD) turbulence revealed by numerical simulations of accretion disks (Abramowicz et al. 1996; Armitage 2002; Steinacker & Papaloizou 2002; Pessah et al. 2008) which show that angular momentum transport is inefficient in regions of the disk where dΩ/dr > 0, which are stable to the standard magnetorotational instability (MRI; see Balbus & Hawley 1991, 1998).

Motivated by the need of a deeper understanding of the behavior of an MHD fluid in a differentially rotating background that deviates from a Keplerian profile, we study the dynamics of MHD waves in configurations that are stable to the standard MRI. Employing the shearing-sheet framework, we show that transient amplification of shearing MHD waves can generate magnetic energy without leading to a substantial generation of hydromagnetic stresses. We discuss the implications of these findings.

We focus our attention on a subsonic, weakly magnetized fluid for which the ram pressure and magnetic pressure remain small compared to their thermal counterpart. As a first approximation, we thus consider a differentially rotating fluid with angular frequency $\boldsymbol{\Omega } = \Omega (r) \check{\boldsymbol{z}}$ and constant background density ρ₀. This is a reasonable assumption in light of the results presented by Armitage (2002) and Steinacker & Papaloizou (2002), who carried out numerical simulations of boundary layers of unstratified accretion disks and found that the density fluctuations throughout the simulations are in general quite small.

We work in the framework of the shearing-sheet approximation (Hill 1878; Goldreich & Tremaine 1978; Narayan et al. 1987; Hawley et al. 1995), where the equations describing an incompressible MHD fluid in a corotating frame are given by

$$\begin{eqnarray} \partial _t \boldsymbol{v} + (\boldsymbol{v}\cdot \nabla)\boldsymbol{v} &=& - 2\Omega _0\check{\boldsymbol{z}}\times \boldsymbol{v} + 2 q\Omega _0^2 x\check{\boldsymbol{x}} - \frac{\nabla P}{\rho _0} \nonumber \\ &&+ \frac{(\boldsymbol{B}\cdot \nabla)\boldsymbol{B}}{4\pi \rho _0} + \nu \nabla ^2\boldsymbol{v} , \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \partial _t \boldsymbol{B} + (\boldsymbol{v}\cdot \nabla)\boldsymbol{B} &= (\boldsymbol{B}\cdot \nabla)\boldsymbol{v} + \eta \nabla ^2\boldsymbol{B}. \end{eqnarray} \tag{ 2 }$$

Here, $\boldsymbol{v}(\boldsymbol{x},t)$ and $\boldsymbol{B}(\boldsymbol{x},t)$, with $\nabla \cdot \boldsymbol{v} = \nabla \cdot \boldsymbol{B} = 0$, stand for the velocity and magnetic fields; ν and η denote the kinematic viscosity and resistivity; and Ω₀ ≡ Ω(r₀) is the corotating angular frequency at a fiducial radius r₀. The first and second terms on the right-hand side of Equation (1) correspond to the Coriolis and tidal forces, respectively. The local pressure P can be found by the divergence-less condition of the velocity field. Recalling that the local density ρ₀ is assumed to be constant, and in order to simplify notation, hereafter we redefine the symbols denoting the pressure and magnetic field in such a way that and P/ρ₀ → P and $\boldsymbol{B}/(4\pi \rho _0)^{1/2} \rightarrow \boldsymbol{B}$.

We decompose the flow into mean and fluctuations as

$$\begin{eqnarray} \boldsymbol{v}(\boldsymbol{x}, t) &\equiv \boldsymbol{U}_1(x) + \boldsymbol{u}(\boldsymbol{x}, t), \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \boldsymbol{B}(\boldsymbol{x}, t) &\equiv \boldsymbol{B}_0(t) + \boldsymbol{b}(\boldsymbol{x}, t). \end{eqnarray} \tag{ 4 }$$

The leading order background velocity is $\boldsymbol{U}_1(x) \equiv -q\,\Omega _0 x \check{\boldsymbol{y}}$, where the shear parameter q is given by

$$\begin{eqnarray} &&q \equiv \left.-\frac{d\ln \Omega }{d\ln r}\right|_{r_0} . \end{eqnarray} \tag{ 5 }$$

The homogeneous background magnetic field is, in general, a function of time and it evolves according to the induction Equation (2), i.e., $\partial _t\boldsymbol{B}_0 = - q\,\Omega _0 B_{0x} \check{\boldsymbol{y}}$.

The substitution of Equations (3) and (4) into (1) and (2) leads to a nonlinear system for the dynamical evolution of the perturbations $\boldsymbol{u}(\boldsymbol{x}, t)$ and $\boldsymbol{b}(\boldsymbol{x}, t)$. However, as pointed out in Goodman & Xu (1994), all the nonlinear terms in the resulting equations vanish identically if we consider the evolution of a single Fourier mode.³ In this case, we obtain

$$\begin{eqnarray} \mathcal {D}_t\boldsymbol{u} - q\,\Omega _0 u_x \check{\boldsymbol{y}} &= \boldsymbol{B}_0\cdot \nabla \boldsymbol{b} + \nu \nabla ^2\boldsymbol{u} - 2 \Omega _0 \check{\boldsymbol{z}}\times \boldsymbol{u} - \nabla P, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \mathcal {D}_t\boldsymbol{b} + q\,\Omega _0 b_x \check{\boldsymbol{y}} &= \boldsymbol{B}_0\cdot \nabla \boldsymbol{u} + \eta \nabla ^2\boldsymbol{b}. \end{eqnarray} \tag{ 7 }$$

The "semi-Lagrangian" time derivative $\mathcal {D}_t \equiv \partial _t \,{+}\, \boldsymbol{U}_1\cdot \nabla$ accounts for advection by the shearing background. The shearing component in the Coriolis term cancels out the tidal force. We remark that Equations (6) and (7) are not just linearized equations, they remain valid even if the amplitude of the perturbations is not small compared to the background values, and they are exact as long as a single Fourier mode is considered. Under these conditions, it is sensible to study the evolution of $\boldsymbol{u}(\boldsymbol{x}, t)$ and $\boldsymbol{b}(\boldsymbol{x}, t)$ for a long time.

In order to solve Equations (6) and (7), it is convenient to work in Fourier space. The x-dependence of the "semi-Lagrangian" time derivative can be removed by employing a shearing coordinate system (x', y', z', t') ≡ (x, y + q Ω₀xt, z, t) in which $\mathcal {D}_t = \partial _{t^{\prime }}$ (Goldreich & Lynden-Bell 1965). A single mode with a fixed "shearing" wavenumber $\boldsymbol{k}^{\prime }$ is thus given by

$$\begin{eqnarray} \boldsymbol{u}(\boldsymbol{x},t) &= 2\mathrm{Re} [\hat{\boldsymbol{u}}_{\boldsymbol{k}^{\prime }}(t)\exp (i\boldsymbol{k}^{\prime }\cdot \boldsymbol{x}^{\prime })], \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \boldsymbol{b}(\boldsymbol{x},t) &= 2\mathrm{Re} [\hat{\boldsymbol{b}}_{\boldsymbol{k}^{\prime }}(t)\exp (i\boldsymbol{k}^{\prime }\cdot \boldsymbol{x}^{\prime })], \end{eqnarray} \tag{ 9 }$$

where $\boldsymbol{k}^{\prime }\cdot \boldsymbol{x}^{\prime } = \boldsymbol{k}(t)\cdot \boldsymbol{x} = (k_x^{\prime } + q\,\Omega _0 t k_y^{\prime }) x + k_y^{\prime } y + k_z^{\prime } z$.

Substituting the ansatz (8) and (9) into Equations (6) and (7) leads a set of equations for the Fourier amplitudes

$$\begin{eqnarray} d_t\hat{\boldsymbol{u}} - q\,\Omega _0 \hat{u}_x \check{\boldsymbol{y}} &= i\omega _\mathrm{A}\hat{\boldsymbol{b}} - \nu k^2\hat{\boldsymbol{u}} - 2 \Omega _0 \check{\boldsymbol{z}}\times \hat{\boldsymbol{u}} - i\boldsymbol{k} \hat{P}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} d_t\hat{\boldsymbol{b}} + q\,\Omega _0 \hat{b}_x \check{\boldsymbol{y}} &= i\omega _\mathrm{A}\hat{\boldsymbol{u}} - \eta k^2\hat{\boldsymbol{b}}. \end{eqnarray} \tag{ 11 }$$

Here, we have replaced $\partial _{t ^{\prime }}$ by d_t and omitted the subscripts $\boldsymbol{k}^{\prime }$ in the Fourier coefficients in order to simplify the notation. We have also introduced the (time-independent) Alfvén frequency $\omega _\mathrm{A}\equiv \boldsymbol{B}_0(t)\cdot \boldsymbol{k}(t)$ (see also Balbus & Hawley 1992a).

The pressure term can be eliminated from Equation (10) using the solenoidal character of the velocity field, which implies $d_t i\boldsymbol{k} = i\boldsymbol{k} d_t + iq\,\Omega _0 k_y\check{\boldsymbol{x}}$. This leads to

$$\begin{eqnarray} &&\hat{P} = -\frac{2i\Omega _0}{k^2}[(q-1)k_y\hat{u}_x + k_x\hat{u}_y], \end{eqnarray} \tag{ 12 }$$

which is independent of $\hat{u}_z$. Because we are interested in the transport of angular momentum along the radial direction, this decoupling allows us to solve separately for the x- and y-components.

The dynamical evolution of the modes with $\boldsymbol{k} \equiv k_z\check{\boldsymbol{z}}$ is quite simple; they grow exponentially if k²_zω_A² ⩽ 2qΩ²₀ (Balbus & Hawley 1992b; Pessah et al. 2006). Thus, a Keplerian disk (with q = 3/2) can exhibit exponential growth but a shear profile with q < 0 only supports stable oscillations. In order to isolate the interesting dynamics that could arise from modes that are not associated with the MRI, we thus focus on modes with k_z = 0. Taking the curl of the momentum equation, it is easy to verify that the Coriolis term does not play a role in the equation for the vorticity, and hence it has no effect on the dynamics of the system (Lithwick 2007). This shows explicitly that the standard MRI is absent in our analysis.

We choose the origin of time so that k_x(t) is initially zero. In other words, we use k_x(t) to define our time coordinate,

$$\begin{eqnarray} &&\tau \equiv k_x(t)/k_y \equiv q\,\Omega _0 t , \end{eqnarray} \tag{ 13 }$$

so the divergence-less conditions become

$$\begin{eqnarray} &&\tau \,\hat{u}_x + \hat{u}_y = \tau \,\hat{b}_x + \hat{b}_y = 0. \end{eqnarray} \tag{ 14 }$$

Assuming ν, η ≪ ω_A/k², we can work in the ideal limit and neglect viscosity and resistivity. The x-components of the MHD equations then become

$$\begin{eqnarray} &&d_\tau {\left(\begin{array}{@{}l@{}}\hat{u}_x \\ \hat{b}_x\end{array}\right)} = {\left(\begin{array}{@{}cc@{}}-\Gamma & i\omega \\ i\omega & 0 \\ \end{array}\right)} {\left(\begin{array}{@{}c@{}}\hat{u}_x \\ \hat{b}_x\end{array}\right)}, \end{eqnarray} \tag{ 15 }$$

where all the temporal dependence is contained in the factor

$$\begin{eqnarray} &&\Gamma (\tau) \equiv 2 \tau / (\tau ^2 + 1), \end{eqnarray} \tag{ 16 }$$

and

$$\begin{eqnarray} &&\omega \equiv \omega _\mathrm{A}/ q\,\Omega _0 \end{eqnarray} \tag{ 17 }$$

is the dimensionless Alfvén frequency. The linear system (15) can be recast into one equation as

$$\begin{eqnarray} &&d_\tau ^2 \hat{b}_x + \Gamma (\tau)\,d_\tau \hat{b}_x + \omega ^2 \hat{b}_x = 0, \end{eqnarray} \tag{ 18 }$$

where the dependence on the parameters q, Ω₀, and ω_A is only through the combination (ω_A/q Ω₀)². This second-order differential equation for $\hat{b}_x$ is identical to Equation (2.20) in Balbus & Hawley (1992a) when k_z = 0. In this case, the perturbations in the z-coordinate decouple from the perturbations in the perpendicular direction (see also their Equation (2.19)).

Unfortunately, Equation (18) does not have an analytical solution. However, if we consider the limit τ² ≫ 1, it reduces to a spherical Bessel equation, which possesses as solutions

$$\begin{eqnarray} \hat{u}_x &= S j_1(\omega \,\tau) + C y_1(\omega \,\tau), \nonumber \\ &= - \frac{C + S \omega \,\tau }{\omega ^2 \tau ^2}\cos (\omega \,\tau) + \frac{S - C \omega \,\tau }{\omega ^2 \tau ^2}\sin (\omega \,\tau), \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} \hat{b}_x &= -i S j_0(\omega \,\tau) - i C y_0(\omega \,\tau) \nonumber \\ &= \frac{iC}{\omega \,\tau }\cos (\omega \,\tau) - \frac{iS}{\omega \,\tau }\sin (\omega \,\tau), \end{eqnarray} \tag{ 20 }$$

where j_n(x) and y_n(x), with n = 0, 1, are spherical Bessel functions of the first and second kind, respectively; and S and C are complex constants determined by the initial conditions. Using the divergence-less conditions in Equation (14), the y-components are simply

$$\begin{eqnarray} \hat{u}_y &= - \tau \left[S j_1(\omega \,\tau) + C y_1(\omega \,\tau)\right]\nonumber \\ &= \frac{C + S\omega \tau }{\omega ^2\tau }\cos (\omega \,\tau) - \frac{S - C\omega \tau }{\omega ^2\tau }\sin (\omega \,\tau), \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \hat{b}_y &=\,\,i\tau \left[S j_0(\omega \,\tau) + C y_0(\omega \,\tau)\right]\nonumber \\ &= - \frac{iC}{\omega }\cos (\omega \,\tau) + \frac{iS}{\omega }\sin (\omega \,\tau). \end{eqnarray} \tag{ 22 }$$

Because the pressure in Equation (12) is independent of $\hat{u}_z$, the exact solutions (for all time) are Alfvén waves,

$$\begin{eqnarray} \hat{u}_z &=\,C^{\prime }\cos (\omega \,\tau) +\;S^{\prime }\sin (\omega \,\tau), \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} \hat{b}_z &= iC^{\prime }\sin (\omega \,\tau) - iS^{\prime }\cos (\omega \,\tau), \end{eqnarray} \tag{ 24 }$$

where C' and S' are some other complex constants. We note that, although the direction of the dimensionless time τ depends on the sign of the shear parameter q, the combination ω τ ≡ ω_At is insensitive to it. Hence, given the same initial conditions, $\hat{u}_x$ and $\hat{b}_x$ are symmetric, while $\hat{u}_y$ and $\hat{b}_y$ are antisymmetric, in the shear parameter q.

We demonstrate the accuracy of these analytical approximations in Figure 1, which shows both the numerical and analytical solutions for $\mathrm{Im}[\hat{b}_x]$ and $\mathrm{Im}[\hat{b}_y]$ with q = −3/2 and ω_A = Ω₀ as an example. The initial conditions are set at ω τ = ω_At = −20 by choosing C₋ = 1 and S₋ = 0. The numerical solutions, shown with thick solid lines, result from integrating Equation (18) with the definition (16) and using the divergence-less conditions (14). The dotted lines in the two panels are obtained by setting C = C₋ and S = S₋ in the analytical approximations (20) and (22). These solutions are indistinguishable for τ ≲ −1. As expected, the approximations break down for τ ≃ 0. This is precisely where the numerical solutions change their amplitudes significantly. The analytical expressions (20) and (22) are again in excellent agreement with the numerical solutions for τ ≳ 1, provided that their amplitudes are given by C = C₊ = 0.285 and S = S₊ = −1.896. These constants are found by requiring that both the numerical and analytical solutions match for ω τ = ω_At ≫ 1 (in practice, we set ω τ = 20).

Zoom In Zoom Out Reset image size

Even though our analytical approach cannot predict the change in amplitude close to τ ≈ 0, the solutions that we obtain are a very good approximation to the numerical results as long as τ² > 1. We could, in principle, obtain the coefficients C₊ and S₊ by an asymptotic matching technique similar to the one employed in Heinemann & Papaloizou (2009). However, the analytical solution near τ ≃ 0 contains special functions that are too complicated to be useful. More importantly, as we show below, the most interesting features of the solutions are independent of the precise values of these constants.

Given the solutions (19)–(24) for the Fourier amplitudes, we obtain the (mean) total stress T_xy ≡ 〈u_xu_y − b_xb_y〉 and (mean) energy density E ≡ 〈u² + b²〉/2 of the fluctuating fields,⁴ where the brackets stand for the spatial average (see, e.g., Pessah et al. 2006). Because these solutions are only valid for early/late times, we can approximate the total stress and energy density up to first order in 1/ω τ as

$${\fontsize{9.7}{12}\selectfont{\begin{eqnarray} &&T_{xy} \approx -\frac{2}{\omega ^2\tau} [(|S|^2 - |C|^2)\cos (2\omega \,\tau) +(S^*C + SC^*) \sin (2\omega \,\tau)],\nonumber\\ \end{eqnarray}}} \tag{ 25 }$$

$$\begin{eqnarray} E &\approx& -\frac{1}{\omega ^3\tau } [(|S|^2 \!-\! |C|^2)\sin (2\omega \,\tau) -(S^*C \!+\! SC^*)\cos (2\omega \,\tau)]\nonumber \\ &&+ \frac{1}{\omega ^2} (|S|^2 - |C|^2) + |S^{\prime }|^2 + |C^{\prime }|^2 , \end{eqnarray} \tag{ 26 }$$

where the asterisk denotes complex conjugation. Using these expressions, it is easy to see that the energy balance equation d_tE = qΩ₀T_xy is also satisfied up to order 1/ω τ.

In Figure 2, we illustrate the numerical solutions for the stress T_xy(t) and energy E(t), given by the thin and thick solid lines, together with the analytical approximation for the energy. The latter has been obtained by substituting the two pairs of constants, C₋ = 1 and S₋ = 0, and C₊ = 0.285 and S₊ = −1.896, into Equation (26). It is thus clear that the late-time stress oscillates around zero with decreasing amplitude, while the energy density asymptotes to a non-vanishing, time-independent value. The expression for the energy density at early/late times in terms of the constants S_± and C_± is given by⁵

$$\begin{eqnarray} &&E_\pm \equiv \lim _{t\rightarrow \pm \infty } E(t) = \frac{|S_\pm |^2 + |C_\pm |^2}{\omega ^2}. \end{eqnarray} \tag{ 27 }$$

Therefore, the energy gain via swing amplification, E₊/E₋, is in general a function of the ratio ω = ω_A/q Ω₀ and the initial conditions. However, it is possible to obtain conclusions that are independent of the latter.

Zoom In Zoom Out Reset image size

The dependence of the energy gain on the initial conditions for ω² = 1 is shown in Figure 3. The horizontal axis describes how the initial energy is distributed between the j_n modes and the y_n modes, while the different lines show the phase difference in the corresponding initial amplitudes. When $\arg (S_-/C_-) = \pi /2$ or 3π/2, the j_n and y_n modes are completely out of phase and evolve independently. This results in an energy gain that is linear in the initial amplitudes (thick solid line). We have found that the dependence of the energy gain on the phase difference between the constants determining the initial conditions is weaker if ω decreases below unity. In this case, all the different curves converge to the thick line corresponding to $\arg (S_-/C_-) = \pi /2$. At the same time, as ω decreases below unity, this line gets steeper, providing thus a larger energy gain. This justifies referring to the y_n and j_n as the "growing" and "decaying" modes, respectively.

Zoom In Zoom Out Reset image size

We illustrate the dependence of the energy gain on the shear parameter in Figure 4 (because the results depend only on ω², we only show the positive domain in the horizontal axis). In the limit of weak shear, there are only pure Alfvén waves and there is no net energy gain. The dashed line shows that the energy gain tends to the value E₊/E₋ = 10/ω² as 1/ω ≫ 1, which provides a good description of the numerical results for strong shear. The asymptotic behavior is insensitive to the initial conditions as long as the growing mode is excited, i.e., C₋ ≠ 0.

Zoom In Zoom Out Reset image size

5.1. Summary and Connection to Previous Work

We have employed the shearing-sheet framework to study the dynamical evolution of MHD waves in weakly magnetized differentially rotating backgrounds which are stable to the MRI. While the fact that these waves can be transiently amplified is widely appreciated, our motivation to study them, as well as the results that we obtained, concern dynamical aspects that have not received as much attention. This is whether these shearing MHD waves can play a significant role in the transport of angular momentum in regions of the disk where the MRI is inefficient, such as the accretion disk boundary layer.

The equations that we have solved are similar to those presented in Balbus & Hawley (1992a), who provided numerical solutions showing that transient amplification of MHD waves is a general outcome for the modes with wavevectors that are not exactly aligned with the rotation z-axis.⁶ Analytical insight on the nonlinear dynamics of these waves has been usually hindered by the fact that the governing equations cannot be simplified beyond a set of coupled differential equations (see, e.g., Fan & Lou 1997; Kim & Ostriker 2000; Brandenburg & Dintrans 2006; and also Johnson 2007, where higher order WKB solutions for the linear evolution of the shearing waves are provided). In this paper, by isolating the modes that are unrelated to the standard MRI, i.e., by setting k_z = 0, we have been able to provide analytical solutions that are valid for all times, except close to the instant where the waves evolve from leading to trailing, i.e., when k_x(t) = 0. These solutions are exact when only one Fourier mode is considered (Goodman & Xu 1994).

Despite the fact that we do not analytically predict the amplification factor for these waves, the characteristics of the solutions that we found allowed us to draw important conclusions. We showed that the amplification factor is only a function of the (time-independent) dimensionless ratio (qΩ₀/ω_A)², with

$$\begin{eqnarray} &&\frac{E_+}{E_-} \approx 10 \left(\frac{q\Omega _0}{\omega _{\rm A}}\right)^2 \quad \hbox{for} \ q\Omega _0 \gg \omega _{\rm A} , \end{eqnarray} \tag{ 28 }$$

(see Equation (27) and Figure 4) and it is thus insensitive to the sign of the shear parameter q.

An important result of this study is that while the energy of these MHD waves can be significantly amplified, their net associated stresses oscillate around zero (see Equation (25) and Figure 2). This suggests that these shearing MHD waves are unlikely to play an important role in the transport of angular momentum in the accretion disk boundary layer region. These findings are consistent with the results of global MHD simulations of accretion disks with a rigid inner boundary carried out in Armitage (2002) and Steinacker & Papaloizou (2002). These simulations show that the inner disk regions, where dΩ/dr ⩾ 0, can develop strong toroidal magnetic fields, with associated magnetic energies that can easily reach a few tenths of the thermal energy, without leading to efficient angular momentum transport.

5.2. Implications

5.3. Final Remarks

M.E.P. is grateful to the Knud Højgaard Foundation for its generous support. C.K.C. is supported by a NORDITA fellowship. We thank Tobias Heinemann, John Wettlaufer, and Jim Stone for useful discussions.