Nonlinear stability of the ideal magnetohydrodynamic interchange mode at marginal conditions in a transverse magnetic field

Nonlinear stability of the ideal magnetohydrodynamic interchange mode at marginal conditions in a transverse magnetic field

Jupiter Bagaipo    P. N. Guzdar    A. B. Hassam Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, MD 20742-3511
18 July 2011

The stability of the ideal magnetohydrodynamic (MHD) interchange mode at marginal conditions is studied. A sufficiently strong constant magnetic field component transverse to the direction of mode symmetry provides the marginality conditions. A systematic perturbation analysis in the smallness parameter, , is carried out, where is the critical transverse magnetic field for the zero-frequency ideal mode, and is the deviation from . The calculation is carried out to third order including nonlinear terms. It is shown that the system is nonlinearly unstable in the short wavelength limit, i.e., a large enough perturbation results in instability even if (linearly stable). The normalized amplitude for instability is shown to scale as . A nonlinear, compressible, MHD simulation is done to check the analytic result. Good agreement is found, including the critical amplitude scaling.

I Introduction

It is well known that the magnetohydrodynamic (MHD) magnetized plasma interchange instability can be stabilized by a transverse magnetic field. For a given wavenumber, allowing a magnetic field component in the direction of the wavenumber introduces Alfvénic stabilizing tension such that beyond a critical transverse field (transverse to the direction of mode symmetry) that wavenumber is linearly stable.Freidberg (1987) The nonlinear evolution of the magnetized plasma interchange instability is less well understood. In particular, the state of the system for when the transverse B-field is marginally subcritical (or, equivalently, the plasma beta is slightly above critical) is an important question for magnetized fusion energy applications: does the mode saturate at low amplitude and how does the marginal convection and resulting transport scale with deviation from marginality? The question is an important consideration for stellarators, for example, since these fusion devices are engineered for very high precision magnetic fields and one of the precision constraints arises from ideal MHD linear stability results.Kulsrud (1963) If the nonlinear consequence of a slightly subcritical B-field were better understood, it may be possible to optimize over the MHD design constraints. It was also recently shown that the linear growth rate of ideal interchanges in a reversed-field pinch for a slightly subcritical B-field is weaker than expected and may be overcome by nonlinear effects.Gupta, Callen, and Hegna (2002)

The interchange mode is a pressure-driven mode that is characterized by the interchange of magnetic flux tubes so that the overall free energy of the system is lowered.Rosenbluth and Longmire (1957) The instability occurs when the equilibrium has a density gradient unfavorable to the direction of a “gravitational” force. In systems with curvature, this force comes from a centrifugal force generated by thermal motion in field curvature. The mode can be stabilized by introducing a strong enough field, transverse to the flutes, that prevents the flux tubes from being able to freely interchange. The strength of the stabilizing field can be determined using linear theory and will depend on the steepness of the density gradient and the magnitude of the gravitational force.

There have been a few studies done on nonlinear growth of interchange instabilities at marginal stability in tokamaks.Rutherford, Furth, and Rosenbluth (1971); Waelbroeck (1989); Beklemishev (1991); Cowley and Artun (1997); Zhu, Hegna, and Sovinec (2006); Zhu et al. (2007) Although a Lagrangian approach has been attempted,Pfirsch and Sudan (1993) the general approach is to expand the equation of motions of the unstable mode about marginal stability and thus the nonlinear terms in the system can be evaluated.Waelbroeck (1989); Beklemishev (1991) We can determine the overall stability of the system by comparing the behaviour of the nonlinear effects to the linear driving term. In Ref. Waelbroeck, 1989 the author found that, for the profiles investigated, the nonlinear effects were stabilizing. Similarly, in Ref. Beklemishev, 1991 the author showed nonlinear saturation at marginal stability. Both authors considered a system with a sheared magnetic field. In studying the line-tied mode, the authors in Ref. Cowley and Artun, 1997 showed that near the marginally stable point the system was nonlinearly unstable. However, Refs. Zhu, Hegna, and Sovinec, 2006; Zhu et al., 2007 showed that the nonlinear growth transitions through an initial regime where the nonlinear growth dominates the linear response, as shown in Ref. Cowley and Artun, 1997, but a secondary regime takes over when the amplitude is sufficiently large and so the mode amplitude remains bounded.

We simplify our system to a slab geometry where we use an effective gravitational field, , to model centrifugal force due to field line curvatureRosenbluth and Longmire (1957) and we assume a constant transverse field. This reduces the complexity of the system so that the focus of the analysis can be on how nonlinear terms get introduced into the equations of motion. The idealized system is described in Sec. II along with the derivation of nonlinear time evolution equation. The goal is to have a simpler methodology in a simple system that can be generalized into more complicated systems, e.g. sheared fieldSuydam (1958), ballooningDrake and Antonsen (1984); Connor, Hastie, and Taylor (1978), etc. In Sec. III we verify our result using a dissipative numerical simulation. The results are summarized in Sec. IV.

Ii Theory

Consider a slab system with constant gravity and very strong magnetic field in the -direction such that and is much larger than , the typical flow in the system. This system is incompressible and can be described by the two-dimensional MHD reduced equationsStrauss (1976) given by,


where and and we have defined, in the usual way,

Variations in are suppressed since the fastest interchange has . The nonlinear system of equations (1)-(3) can be solved for the variables , the density, and and , the flow and magnetic streamfunctions, respectively.

We consider a static equilibrium with a density gradient that’s unstable to interchange and a constant, transverse magnetic field. More explicitly, we have


where henceforth the primes denote differentiation with respect to . We also add the assumption that at the boundaries and has even parity.

Small perturbations about this equilibrium yield the WKB dispersion relation


where is the Rayleigh-Taylor growth rate and is the wavenumber in the direction. The Rayleigh-Taylor growth rate represents the effective “gravitational” acceleration and is the driving force in an interchange instability. With a transverse magnetic field, field line bending results in Alfvénic restoring forces with frequency . In this paper, we consider the dynamics of the magnetized interchange mode when the magnetic field strength is strong enough to just stabilize interchanges, i.e., the system is near marginal stability. In particular, for a given , suppose everywhere in except for a single small region where it is very close to zero, positive or negative. In that case, weakly growing perturbations are possible in the vicinity of where is close to zero. The time rate of change of the perturbations will be very small compared to the local . Thus, we order


This implies that any deviations in away from criticality must be small. In particular, if


then, according to (5), must be of .

We allow small perturbations about this marginal point such that the amplitude of the magnetic perturbation, , while small, is large enough that the nonlinear magnetic tension forces can influence the growth time. This results in the optimal ordering


We represent the perturbation by expanding and in a series


where the order in is denoted by the subscript. The continuity equation (1) can be satisfied by letting and using (3). With this change of variable we can expand in terms of to get


Substituting the expansions (9)-(11) into (2) and (3) we can solve for the nonlinear evolution of the perturbations order by order.

ii.1 First order equations

Matching terms to lowest, non-vanishing order, we obtain from (2) and (3) the equations


where (11) gives


Substituting into (12) the equation becomes


where we have defined the operator

Writing as


we obtain the eigenvalue equation


that can be solved to get an eigenvalue for . The boundary condition for implies that decays exponentially close to the boundary.

In writing (16) we assumed a cosine perturbation in the density which implies from (11). If we also assume that this initial perturbation results in a pure mode for the lowest order flow then


is the solution to (13).

To the lowest order we have found that given the mode of the density perturbation, , and the equilibrium density gradient profile, , then the marginally stable field strength can be solved for using (17). This result is consistent with the prediction from linear theory for the existence of the marginally stable value.

ii.2 Second order equations

In order to solve for the time evolution of , it is necessary to proceed to higher order in the expansion. We now match terms in equations (2) and (3) to get


where (11) to the same order gives


Using from (16), can be solved for in (20) to obtain


where is a constant of integration.

Substituting (21) into (19) results in an equation for ,


where we have used (15) to simplify the Laplacian. This has a solution of the form


where is the homogeneous solution to (23). Substituting (24) into (23), we find satisfies


To fully analyze the stability of our system we still have to resolve the time-evolution of . It is also important to solve for and to make sure that those terms are well-behaved.

ii.3 Third order equations

As was done previously in lower orders, we match terms of in (2) and (3). The resulting higher order equations are


along with


from (11).

Integrating (II.3) over one period in we find that is not driven by so we can set


without loss of generality. No zonal flows are generated in the system when creating a periodic perturbation in the density. However, averaging (27) over we find that zonal fields are generated according to


For a given and the system is now solved up to second order with the exception of the time-evolution . The variables , , , and are defined by (16), (18), (24), and (22), respectively. We can solve for using (17) and then for using (25). The -independent terms and are given by (29) and (30).

To solve for we need to simplify (II.3) by making use of (15), (20), and (28). After some algebra (II.3) takes the form


where exact details of the functional is suppressed here for clarity but is shown in Appendix A. We can extract a time-evolution equation by substituting (16), (24), and (30) into the above equation and applying the operator evaluated over all space. This operation will annihilate the term and any higher order harmonics.

After simplification (see Appendix A), we arrive at the equation for


where the angled brackets are defined as

with evaluated at . We can simplify this further by letting

in order to introduce dimensionless versions of the variables and , and have with dimensions of . Applying this normalization to (II.3) we get


where and


where the primes and brackets now denote derivatives and integrals in . Using the same normalization on (17) and (25) we get the following equations,


for the dimension-free and .

The time-evolution equation (33) closes the system and we can fully determine the first and second order perturbations, and defined by (16) and (24), for a given , , and . This is achieved by first solving the eigenvalue problem (36) and using the solution for and to solve for using (37), and finally determining the coefficients (34), (II.3) and solving for in (33).

The coefficient is a positive number for , and so the linear stability of the system is determined by the sign of . This result agrees with the linear theory. However, the overall nonlinear stability of the system is going to be determined largely from the sign of compared to the sign of .

ii.4 Short wavelength limit

We can analytically solve (36) for the case in which regime the cells are elongated in -direction but still shorter than the scale of the gradient, i.e.,

With this scaling we can approximate to be


Assuming that , then from scaling arguments we find that (36) has the familiar form of a quantum harmonic oscillator. This has the well-known solution


for the ground state. This solution is correct only for and the solution for the “energy” adds a small correction to the initial assumption. Using the same scaling, to lowest order, (37) has the solution


The time-evolution equation (33) can be simplified in the limit by substituting the solutions (39)-(41) in the coefficients (34) and (II.3). After simplification we arrive at the following values for the coefficients


where we only kept the largest terms and have assumed that .

The above result implies that even if , if the initial amplitude is such that


then the system will be nonlinearly unstable and the amplitude will increase without bound. Furthermore, for the instability grows faster than predicted from linear theory and any small perturbation will continue to grow larger without saturation.

With the solution for the eigenmode we can check the ratio between the spatial scale of the perturbation, characterized by the displacement in the -direction , and the width of the eigenmode


given by (39). The displacement is related to the velocity such that , and from (20) we get that , which implies that . Substituting for using (43) gives us a scale for the displacement,


which yields


for the ratio of the two scale lengths. As should be expected, the spatial size of the amplitude required to be nonlinearly unstable is much smaller than the width of the eigenmode.

Figure 1: The equilibrium profiles for the background field , the density , and the magnetic streamfunction along with the difference from constant field.

Iii Numerical Simulation

To confirm this result, we used a two-dimensional code that solves the fully compressional equations (see Appendix B). The variables , , , and are solved numerically and stepped in time. We set so the equations are effectively reduced. The code is dissipative so we introduced source terms in the density in order to maintain a steady state profile suitable for our model. The sourcing, although weak, results in a profile for . To compare with analytic theory, we wish to keep approximately constant. Thus, we allowed to resistively relax at a somewhat slower rate than in the equilibrium.

The system is normalized so that initially and , where is the height of the box. We used hardwall, free-slip boundary conditions for the top and bottom walls and periodic boundary conditions for the sides. The periodic boundary conditions discretize the system so that the only wavenumbers allowed are integer multiples of , where is the width of the box. From (5) we know that the lower modes are the most unstable, so to study the case with , i.e. short wavelength, we selected such that the minimum value for satisfies this condition. By choosing we can satisfy the marginality condition by adjusting and/or such that for the minimum mode. We set , and from the density profile we have and so we satisfy the condition

which is necessary to compare with the analytical result from Section II.4. We attempted to run tests with a larger value of by decreasing , but the code was numerically unstable for smaller box widths.

To generate the equilibrium we initialize to and let the system reach an equilibrium which is steady state. The density source term results in a weak flow in the -direction. This flow scales with the diffusion, so a minimal, numerically-stable value for the diffusion is chosen to minimize its effect. The equilibrium profiles for the density and the background field generated are shown in Fig. 1. It is important to note that the equilibrium profile for the density does not have at the boundaries. The boundary conditions imply that at the wall.

After the equilibrium is made, a density perturbation is introduced with . From (14) we can relate the density perturbation amplitude, , to the perturbation amplitude of , i.e. . In Fig. 2 we show the resulting unstable eigenmode developing for the density. For tests done with far away from marginality, i.e. , there was excellent agreement for the growth rate/frequency in the simulation with (5). The theory predicts that there will be nonlinear coupling to the mode with wavenumber , so it is important that this mode and higher modes are allowed. Since the diffusivity is weak, it is ensured that this is the case.

Since we can adjust both and to achieve marginal stability, we decided to fix the value of at , and adjust . With this value of we expect that based on (40). However, we found that an equilibrium with is stable to perturbations as large as in the simulation. We decreased the strength of the transverse field until it became unstable to perturbations with . This value was at and we took this to be the critical value of the transverse field for the numerical simulation. Since the critical amplitude scales like the square root of the deviation from marginality, we are limited to perturbations only as small as otherwise smaller perturbations would have meant having deviations that are close to the limits of our computational power.

Figure 2: The linear growth of an unstable localized mode cut at and for . Time traces separated by are shown.
Figure 3: Result of stability test for a range of deviations from and magnitude of perturbation, . Stable and unstable results are denoted by a circle or a cross, respectively. The solid line is the theoretical boundary.

We created multiple equilibria with different transverse field strength within 10% of the numerical critical field strength. These equilibria were then perturbed with of different orders of magnitude. The result of the test is shown in Fig. 3 where circles and crosses mark stable and unstable points, respectively, and the solid line is for , from our theory, using the parameters from the numerical simulation. The slope of the theory line seems consistent with the numerical data, however, the theory requires larger for nonlinear instability. This inconsistency could be due to the diffusion in the code and, in particular, the resistivity may allow for slippage in the magnetic field lines which can shift the stability boundary at marginal stability. We can calculate the scale size of this shift based on the values used in the simulation (see Appendix B),


This implies that there could be a shift in of order . At marginal stability, even small diffusion can cause significant shifts in stable-unstable boundaries. However, this implies a shift in ; it is harder to explain why resistivity results in a nonlinear instability at large amplitude of perturbation. It is possible that diffusive effects may affect the critical amplitude for nonlinear instability, but the existence of a nonlinear instability phenomenon is harder to explain as a diffusive effect.

In addition to checking the perturbations for a growing linear mode, we also check the time trace of the amplitude for nonlinear effects. In Fig. 4 we show a time trace of the amplitude of , , for the same but different . We can see that the behaviours are different for the two cases. In the unstable case, Fig. 4a, the density perturbations become very large quickly and eventually dissipate after it hits the boundaries (). The time trace of shows that the density profile flattens out () after reaching a peak. So, even though our analysis in Sec. II is only valid as long as we can see from the trace that it continues beyond this limit until the profile collapses. The stable case, Fig. 4b, has an initial growth eventually hitting a peak and then has stable oscillations. Even though the amplitude increases some, it is still small and the density profile holds. This can be seen from the fact that is staying constant the entire time. We can see in Fig. 5 that as we increase further from marginality, this initial growth decreases in magnitude. It also develops faster and has more noise that is indicative of a transient oscillatory mode.

Figure 4: Time trace of the amplitude of density perturbations (solid line) and derivative of the density (dashed line) for with (a) and (b) .

Iv Summary and Conclusions

In this paper, we studied the nonlinear behaviour of a marginally stable interchange system. We used the reduced equations to find an analytic solution near marginality given a density profile, , deviation from marginality, , and wavenumber of perturbation, , of the B-field. The result is a nonlinear differential equation for the amplitude of the density perturbations as a function of time. The threshold for nonlinear instability is dependent on the above quantities, along with . The principal finding of this paper is that marginally stable interchange modes in a magnetized plasma can be nonlinearly unstable for large enough initial perturbations. We arrived at this result from a systematic asymptotic expansion about marginality in the smallness parameter, , carried out to third order. The first order solution can be found using the linear eigenvalue problem. This solution is then used as a source for the second order problem. The third order analysis yields the equation for the time-dependence of the perturbation. We found that the stability of the solution can be determined by calculating the coefficient of the nonlinear term in the differential equation. This is a nontrivial task for a general perturbation, but we could analytically solve this in the short wavelength limit. In this limit we found that the nonlinear coefficient had a positive sign. This meant that in the linearly stable case () it was possible to be nonlinearly unstable if the initial perturbation was large enough. We found the critical amplitude to be proportional to .

Figure 5: Time trace of the amplitude of density perturbations (solid line) and derivative of the density (dashed line) for with .

A nonlinear numerical MHD simulation fully confirms the analytic result. We have used a numerical simulation of the nonlinear, full, compressible, MHD equations with small dissipation to verify our analytical result. We showed very good agreement between the simulation and the theory for deviations, , from of up to 10%. The numerical results show that in the short wavelength limit the system is nonlinearly unstable. There is some disagreement in the time evolution of the density with the analytical result, but this is possible since the analytic calculation is for an ideal system with no dissipation. We also discussed why a shift in for the linear instability threshold, due to dissipation, is possible at marginal stability and how it is harder to explain why the nonlinear result has an amplitude dependent stability. Furthermore, the dependence is cubic so the mode grows without bound once it is unstable. This is even harder to explain as a resistive effect.

It should be noted that the fully analytic calculation is facilitated by using a very simple form (a constant) for the transverse stabilizing magnetic field. So, while the conclusions of this paper seem to be on solid ground, the application of these findings to various systems, to the extent that the transverse field of this paper is very special, must be appropriately qualified. For example, in tokamaks and stellarators, the interchange mode arises on rational surfaces which corresponds to a slab model with a sheared magnetic field vanishing at . In the solar coronal case, line-tying is an important characteristic absent in our simple case. Nonetheless, the conclusions are sufficiently dissimilar as to indicate further investigation. Thus, for example, a neighboring nonlinear saturated state for the interchange mode was found in Refs. Waelbroeck, 1989; Beklemishev, 1991 – whereas the corresponding result in our case, for , indicates a robustly growing mode with no nonlinear saturation. Of course, the transverse magnetic field in these papers was a sheared field with a rational surface for the unstable wave mode. Attempting a marginal stability analysis for sheared field, similar to that used in the present paper, is not straightforward. The fact that the sheared field goes to zero as goes to zero means that a new inner ordering is required, which makes the calculation more involved.

Our results are more consistent with the nonlinear instability found in Ref. Cowley and Artun, 1997 where the authors were also in the parameter range with , , and . It should be noted that their analysis was for the 3D line-tied mode with no transverse field at marginal stability. Even so, the suprising result is that in both cases the system takes off once it becomes nonlinearly unstable. This occurs even when the linear term is stabilizing. The primary difference between the results is the amplitude dependence of the nonlinear term. In Ref. Cowley and Artun, 1997 the nonlinear term has a quadratic dependence, while our analysis yields a cubic dependence on amplitude. If we construct an effective potential, we observe that the result from Ref. Cowley and Artun, 1997 indicates a dependence on the sign of the perturbation at the metastable boundary, while our potential is symmetric in . Another difference is that the result in Ref. Cowley and Artun, 1997 was somewhat mitigated by Refs. Zhu et al., 2007; Zhu, Hegna, and Sovinec, 2006 in that the latter papers argued that the ordering giving nonlinear growth would break down at small amplitudes before the instability fully takes off. In our case, our numerical simulations seem to show, in agreement with analytic constraints, that the nonlinear instability growth continues without bound and the theory only fails when (as saturation is reached).

Our results could also be relevant to tokamak ballooning modes to the extent that these modes are stabilized by an “average minimum- well” and thus always have some parallel wavenumber. Work is in progress to quantify this better. Finally, our results also indicate a closer look at interchange stability in stellarators, presumably in average minimum- stabilized systems.

Further investigation is necessary to answer some questions regarding the results found in this paper. The transient initial growth in the time traces, mentioned in Sec. III, needs to be explained. The change in the growth rate as the system gets closer to marginal stability, with , needs to be investigated and compared to the results from Ref. Gupta, Callen, and Hegna, 2002.


This work was supported by the U.S. Department of Energy. JB and ABH dedicate this manuscript to the memory of Dr. Parvez Guzdar.


Appendix A

Derivation of (33)-(ii.3)

In simplifying (II.3), we found the functional, , to be


The above equation can be simplified by writing and a certain way. From (24) we can write




Writing in this way, we get the following results


where we have written and


The result (A5) can be derived by multiplying (25) with and recombining the terms. Similarly, if we multiply (17) by , we find that


Since for all orders, we get that


and therefore


Similarly, since


then it follows that


We can now simplify the last two terms in (Derivation of (33)-()). Using (A2) we have


where we used (A5), and took advantage of the fact that has no dependence, to remove the Laplacians. Similarly, we use (A8) to get


where we used (A4) to get the second term.

Combining (Derivation of (33)-()) and (Derivation of (33)-()) we can use (Derivation of (33)-()) and (Derivation of (33)-()) to further simplify the terms with a gradient operator. So finally we get


We can also rewrite the first term of (Derivation of (33)-()),


We can now substitute for , , , and using (16), (30), (A3), (A6), and (A7). As described in Section II.3, we use the operator on (II.3) in order to extract the terms that have a dependence. The other terms will be irrelevant since the integration will evaluate to zero if the dependence doesn’t match. And so we find that


Finally, we use the operator on the above equation to get


We made use of the fact that decays exponentially at the boundaries to combine the three terms proportional to in (Derivation of (33)-()) into one term through integration by parts.

To complete the derivation of (33)-(II.3) we still need to simplify the rest of the terms. It is easy to see that after using the annihilation operator then we get


Applying the same operator, we find that