- Email: [email protected]
Two-step control of wall mode and the monodromy matrix
Physics Letters A 307 (2003) 304–312 www.elsevier.com/locate/pla
Two-step control of wall mode and the monodromy matrix ✩ H. Tasso a,∗ , G.N. Throumoulopoulos a,b a Max-Planck-Institut für Plasmaphysik, Euratom association, D-85748 Garching bei München, München, Germany b University of Ioannina, Association Euratom-Hellenic Republic, Department of Physics, Section of Theoretical Physics, GR 451 10 Ioannina, Greece 1
Received 21 October 2002; received in revised form 10 December 2002; accepted 12 December 2002 Communicated by F. Porcelli
Abstract The Floquet stability of systems of differential equations with piecewise constant periodic coefficients is considered. In the “two-step” case the monodromy matrix is the product of two matrix exponentials. It can be evaluated either by reducing the matrix exponentials to polynomials on the eigenvalues of the corresponding matrices or, without knowledge of the eigenvalues, by using the Baker–Campbell–Hausdorff formula. The two methods are discussed and applied to several examples including the “two-step” dissipative Hill’s equation introduced recently H. Tasso, G.N. Throumoulopoulos [Phys. Plasmas 9 (2002) 2662] to stabilize the “resistive wall mode” in magnetically confined plasmas. Application to large systems and matrices is not only possible but needed for the full understanding of the dynamic stabilization of the resistive wall mode. Its practical implementation is discussed. 2002 Elsevier Science B.V. All rights reserved.
1. Introduction In a recent publication [1] we proposed to stabilize the “resistive wall mode” by modulating in time the resistivity of the wall. An actual implementation of this modulation can be produced by, e.g., a second nearby wall having cuts and electronic switches between the sides of the cuts as explained in the last paragraph of Section 4 of [1]. To study the stabilizing effect of this modulation we introduced the concept of an “effective potential” depending upon the instantaneous value of the resistivity. This led us to concentrate on the rather simple dissipative Mathieu–Hill equations, which resulted in very encouraging results related to the modulation strength of wall resistivity and the dissipation in the plasma. The physical situation can be well understood by assuming viscous magnetohydrodynamics in the plasma, then a vacuum region inside a resistive wall with Ohm’s law, and a vacuum region outside the wall extending to infinity or to a hypothetic perfectly conducting wall at a large but finite distance. The equations for linearized perturbations ✩
Presented at the 44th Annual Meeting of the Division of Plasma Physics of the APS, 11–15 November 2002, Orlando, FL, USA.
* Corresponding author.
E-mail address: [email protected] (H. Tasso). 1 Permanent address.
0375-9601/02/$ – see front matter 2002 Elsevier Science B.V. All rights reserved. doi:10.1016/S0375-9601(02)01769-3
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
can be written in a compact way as follows ξ˙ −µ∇ × ∇× −FMHD ∂ ξ = 1 0 ∂t 0 β A
ξ˙ α ξ , 0 −η(t)∇ × ∇× A
305
(1)
are the Lagrangian displacement vector inside the plasma and the potential vector inside the wall, where ξ and A respectively. µ is the plasma viscosity taken here as a scalar, η(t) is the time dependent resistivity of the wall. FMHD is the magnetohydrodynamic stability operator, α and β are nonlocal operators connecting the displacement vector at the plasma boundary with the vector potential at the surface of the wall, the gauge being chosen as the one of zero electrostatic potential. For details concerning FMHD and viscosity see, e.g., [2]. It is obvious that the Floquet stability of Eq. (1) is not easy to solve in real geometry of magnetic confinement devices like tokamaks, reversed field pinches or stellarators. As suggested in [1] this is a formidable task, which should be preceded by a through investigation of simple and relevant models. Discretization of Eq. (1) allows, however, to define a monodromy matrix as introduced in the next section. Here, preliminary investigation of the monodromy of simplified models, more advanced than those of Ref. [1], is carried through, especially for a “twostep” dependence of the modulated parameter.
2. Monodromy matrix Let Y be a n-vector and A(t) a n × n matrix depending periodically upon time in a piecewise constant manner A = A1 ,
0 < t < ts ,
(2)
A = A2 ,
ts < t < T ,
(3)
where ts is the time at which a jump in A occurs, and T is the period. A1 and A2 are time independent n × n matrices. Define now the equation with piecewise constant periodic coefficients as ˙ = A(t)Y. Y
(4)
The solution of (4) can be written, in general, as Y(t) = M(t)Y(0),
M(0) = I
(5)
with the matrix M(t) obeying ˙ = A(t)M. M
(6)
For the two-step case defined by (2) and (3) we have M(t) = exp tA1 , 0 < t < ts , M(t) = exp (t − ts )A2 exp(ts A1 ),
(7) ts < t < T .
(8)
Since the monodromy matrix is defined by M(T), it can be written [3] after integrating (6) over a period as M(T ) = exp (T − ts )A2 exp(ts , A1 ). (9) Since A1 and A2 do not commute, in general, it is not easy to find explicitly the product of two matrix exponentials, especially for matrices of large order. In case the eigenvalues of A1 and A2 are known, it is possible to replace each exponential by a matrix polynomial on the eigenvalues of the corresponding matrix by applying the Cayley–Hamilton (CH) theorem as done implicitly in [3] or, more generally, as explained in [4]. This method
306
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
gives exact results at least for low order ( 4) matrices, though the calculations may become cumbersome and then should be done symbolically by computer. Finally, let us recall that the system is stable if and only if all the eigenvalues of the monodromy matrix (9) have an absolute value 1.
3. Baker–Campbell–Hausdorff formula For large matrices, for which the eigenvalues are not known explicitly, it may be convenient to use the Baker– Campbell–Hausdorff (BCH) formula (see, e.g., [5–7]), which states that M(T ) = exp C
(10)
with
ts (T − ts ) t 2 (T − ts ) [A2 , A1 ] + s [A2 , A1 ], A1 2 12
ts (T − ts )2 + [A1 , A2 ], A2 + · · · , 12 which, in the “symmetric” case ts = T /2, becomes C = (T − ts )A2 + ts A1 +
(11)
T T2 T 3 (12) (A1 + A2 ) + [A2 , A1 ] + [A2 , A2 ], A1 + [A1 , A2 ], A2 + · · · , 2 8 96 where [x, y] denotes the commutator of x and y. The matrices A1 and A2 are identical up to entries depending of the two-step control function. For a vanishing jump the two matrices commute, which means that their commutator must be a polynomial in the jump η2 − η1 . Since C is the log of the monodromy matrix M(T ), the stability of the system demands that the real parts of the eigenvalues of the matrix C be negative or zero. It is difficult to answer this question for finite jumps since C is given by an infinite series. Throughout this Letter the BCH formula is, however, taken as in (11) and (12), i.e., to the third order. C=
4. Two-step dissipative Hill equation For small jumps we keep in (12) the three first terms, for an example, the “two-step” dissipative Hill equation [1] x¨ + d x˙ + b(t)x = 0, or in the form −d y˙ = 1 with
−b(t) y, 0
x˙ y= . x
A1 and A2 are given by −d −b1 A1 = , 1 0 −d −b2 A2 = . 1 0
(13)
(14)
(15)
(16) (17)
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
Inserting (16) and (17) in formula (12) we obtain 2 2 2 + d T8 (b1 − b2 ) − −dT + T8 (b1 − b2 ) −T b1 +b 2 C= 2 T − T8 (b1 − b2 )
T3 48 (b1
− b 2 )2
307
.
(18)
The eigenvalues of C as given by (18) are
4 −dT ± d 2 T 2 − 2T 2 (b1 + b2 ) − T48 (b1 − b2 )2 λ= (19) . 2 We see that, if (b1 + b2 ) > 0, the system is stable with and without dissipation. This shows that the corresponding unstable regions found in the charts displayed in Figs. 3 and 5 of Ref. [1] need higher order terms in the BCH expansion for their description. For (b1 + b2 ) < 0, however, we can find the “parabolic” beginning near the origin for the lowest curve of the same chart by requiring that the real part of the unstable eigenvalue be zero as follows (b1 + b2 ) T 2 (b1 − b2 )2 = . (20) 2 48 4 Let us notice that (b1 + b2 )/2 is to be set equal to K of Ref. [1] and (b1 − b2 )/2 should be set equal to 2k of Ref. [1], T being equal to 2π/Ω. It will not be easy to obtain the full charts of Ref. [1] since we should be able to take into account the higher order commutators in the BCH formula either by symbolic computation or by summing up the series. −
5. Asymmetric case or ts = T /2 The extension of the previous calculation to the time-asymmetric case can be done on the same line. Since the algebra is rather cumbersome the elements of matrix C and its eigenvalues were obtained via symbolic computation. Explicit expressions are given below. −12 dT − (b1 − b2 )(T − ts )ts (−6 + dT − 2 dts ) , 12 1 −12b1ts − (T − ts ) 12b2 + (−b1 + b2 ) d(6 − dT ) + 2T b2 ts + 2(b1 − b2 ) −d 2 + b1 + b2 ts2 , C12 = 12 (b1 − b2 )(T − 2ts )(T − ts )ts , C21 = T + 6 (b1 − b2 )(T − ts )ts (−6 + dT − 2dts ) C22 = 12 C11 =
and λ=
√ 1 6 dT ± ∆ , 12
where ∆ = 36T 2 d 2 − 4b2 + 144T (−b1 + b2 )ts + T 2 (b1 − b2 )2 −12 − d 2 T 2 + 4T 2 b2 ts2 − 2T (b1 − b2 )2 −3 4 + d 2 T 2 + 2T 2 (b1 + 5b2) ts3 + (b1 − b2 )2 −12 − 13 d 2T 2 + 4T 2 (4b1 + 9b2) ts4 + 4T 3d 2 − 5b1 − 7b2 (b1 − b2 )2 ts5 − 4 d 2 − 2b1 − 2b2 (b1 − b2 )2 ts6 .
308
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
Fig. 1. The lowest marginal curves for two-step Hill’s system obtained analytically by means of the BCH formula (11) and numerically on the basis of the CH method associated with the monodromy matrix (9). (a) Shows curves for symmetric in time modulations (ts = 1) while (b) and (c) are associated with asymmetric modulations, i.e., ts = 1.3 and ts = 0.7, respectively.
Introducing physically plausible parameters, i.e., the modulation average b1 ts + b2 (T − ts ) T and the deviation of b1 from K ts b1 ts + b2 (T − ts ) ts , k = (b1 − K) = b1 − T T T K=
and requiring that the real part of the unstable eigenvalue be zero yields K=
k 2 T 2 (12 + (d 2 + 8k)T 2 − 4(d 2 + 4k)T ts + 4d 2 ts2 ) . 4(−6 + kT 2 − 2kT ts )(6 + kT 2 − 2kT ts )
(21)
The lowest marginal curves determined by (21) are shown in Fig. 1 along with the respective exact ones. The latter were obtained numerically on the basis of the CH method; in particular, to determine the marginal curves the eigenvalues λ of the monodromy matrix (9) were calculated numerically and then the marginal b values via the relation |λ| = 1. The results were checked by means of the dispersion relation (17) of Ref. [1]. BCH and CH symmetric curves are also presented in Fig. 1(a) which compare well each other for values k < 1 but they deviate for larger ones. The deviation is affected rather strongly by the asymmetry as can be seen in Fig. 1(b) and (c) thus indicating that higher order terms are required in the BCH expansion in particular for large asymmetric modulations.
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
309
6. Jumps on the diagonal For diagonal matrices A1 and A2 all the commutators in (11) and (12) vanish, which means that there is no proper “dynamic stabilization” in contrast to the case of the Hill’s equation. For nondiagonal matrices having jumps on the diagonal, it may be expected that the jumps have to be large enough to be able to stabilize. It is instructive, as a crude model of Eq. (1), to consider the following problem: −d −b α y, y˙ = 1 0 0 (22) β1 β2 −η(t) with x˙ y= x . (23) u A1 and A2 are given by −d −b α 0 0 , A1 = 1 (24) β1 β2 −η1 −d −b α 0 0 . A2 = 1 (25) β1 β2 −η2 As mentioned above the jump is now in η, which is the third element of the diagonal. Let us now calculate C as given by Eq. (12) T 2 (η2 −η1 ) T 3 (η2 −η1 )2 C=
−dT
−bT
T
0
αT +
0
+
T 2 (η2 −η1 ) T 3 (η2 −η1 )2 T β1 − β1 + β1 8 96
T 2 (η2 −η1 ) T 3 (η2 −η1 )2 T β2 − β2 + β2 8 96
η +η − 12 2
8
96
. (26)
The eigenvalues of matrix (26) are given by a cubic equation, in which the jump appears with even powers, the coefficients of the second power being the strongest as displayed in Eq. (27).
T (η1 + η2 ) λ(λ + dT ) + bT 2 λ+ 2 αβ1 T 2 (η2 − η1 )2 αβ2 T 2 (η2 − η1 )2 − λ αβ1 T 2 + (27) + · · · − T αβ2 T 2 + + · · · = 0. 192 192 Eq. (27) shows that the stability will not be affected by small or moderately large jumps in contrast with the “two-step” Hill’s equation. Marginal curves for symmetric modulations respective to the lowest two-step Hill ones (which can be recovered for α = 0) were also obtained numerically for β1 = 0 and d = 1. Inspection of (27) and the numerical results reveal the following characteristics of the marginal curves: (1) The marginal points are determined by given product K = αβ2 regardless of the particular values of α and β2 . (2) For a “perfectly conducting wall”, η = 0, the marginal b take larger values as either α or β2 increase. (3) The larger the resistivity η1 = η2 ≡ η the smaller the marginal b as can be seen in Fig. 2(a) and (b) for η ≡ η2 − η1 = 0, i.e., for η = 0.6 and η = 1 the respective marginal b values are nearly 0 and −1. It is also noted that for η = 0 the marginal b is 1. (4) There are parametric regimes for which the marginal b take negative values and decrease with increasing η (Fig. 2(a)).
310
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
Fig. 2. BCH and CH marginal curves for a wide variation of the symmetric in time resistivity modulation parameter ∆η = η2 − η1 and different combinations of η1 and the wall parameters α and β2 . The system is stable just above the curves and unstable below.
(5) The monotonicity of the curves depends on the sign of K as follows: • For K < 0 the curves can vary nonmonotonically as η increases. For example, as shown in Fig. 1(a) there is a sharp turning point at η ≈ 1 at which the marginal (b) is nearly −1. For α = −1 and β2 = 1 such a turning point appears only for η1 < 1. For η1 1 the marginal b in the CH case increases monotonically with η remaining, however, negative. The particular case of η1 = 1 is shown in Fig. 2(b). • For K > 0 the CH curves go down monotonically as η increases (Fig. 2(c)). Even for small values of the parameters α and β2 , however, the curves do not involve negative values of b. (6) Similar to the case of the two-step Hill’s system the BCH curves are reliable for η 1 while for larger η they deviate significantly from the CH ones. In order to have a proper interpretation of the above results, let us first state the main features of the actual resistive wall mode described by the full physical equation (1): (1) The mode is an “external kink” with a reduced growth rate due to the resistivity of the wall. (2) The mode should be stable if the wall is perfectly conducting (η = 0) and the wall is not too far from the plasma boundary. (3) The mode should become unstable if the same wall is resistive (η = 0). Resistivity makes the wall “permeable” to the perturbed magnetic field in the full vacuum in whole space or up to a perfectly conducting second wall at large distances. For the 3 by 3 model the parameters α and β2 stand for the importance of the vacuum region while the parameter b represents the perturbation of plasma potential energy, the wall mode being associated with negative
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
311
values of b. If α = 0 and η = 0, the system is marginal for b = 0 (see Fig. 5 of Ref. [1]). Taking α = −1 and keeping zero resistivity makes the system marginal at b = 1. This is representative of the external kink which becomes unstable as soon as the vacuum space is increased which has to be compensated by a higher value of b. So the features (1) and (2) are qualitatively supported by the model. The calculations show, however, that for η1 = η2 = 0.6 the value of b at marginality decreases from 1 to almost zero and for η1 = η2 = 1.0 it decreases further to −1. This would suggest that resistivity of the wall helps stabilizing. In fact, this would be plausible for a wall of zero thickness (“skin” thickness is larger than wall thickness). This may explain also why the agreement between CH and BCH curves in Fig. 2(a) is so good even for η/η1 = 1. At this stage the dynamic stabilization is not visible; this was confirmed by taking in the BCH formula the first term only (no commutators). The existence of a turning point is particularly interesting since the marginal curve has a behaviour similar to what one would expect for dynamic stabilization though quite crude. Indeed, if for example, b is fixed at −0.5 in Fig. 2(a) and η is increased from zero to larger values the system becomes first stable and then for η increasing further becomes unstable again. This behaviour ceases, however, for b < −1. In contrast, the monotonic curves of Fig. 2(b) and (c) do not indicate that a dynamic stabilization of the mode is taking place. The only way to cure the model is to go to larger matrices. This may, of course, be limited by numerical feasibility. The most significant progress is that the method proposed seems to work and could, in principle, be extended in the future to larger problems.
7. Conclusions The CH and BCH methods proposed here to handle Floquet stability of linear systems of differential equations with periodic piecewise constant coefficients are very powerful as demonstrated in the previous sections. For instance, the first marginal curve for “negative energy” modes of the “two-step” Hill equation (13) [1] can be partly reproduced analytically. As mentioned in the introduction, it would be good to assess the results of this model by a more accurate treatment of the dynamic stabilization of the resistive wall mode through time-modulation of the wall resistivity as stated in system (1). In principle, CH and BCH methods can be applied to large systems of differential equations. As an example, the 3 × 3 system (21), a crude mock up of system (1), has been successfully treated in Section 6. The discussion there shows, however, that model (21) is too crude to give all qualitative features of the resistive wall mode. Though dynamical stabilization does not show up as expected because the “wall” has, so to say, zero thickness, the results are very encouraging. This means that we have to go to much larger systems approaching more realistically the plasma and a resistive wall having a thickness larger or comparable to the “skin” thickness. An accurate discretization of system (1) with subsequent application of CH or BCH methods would lead to much higher numerical effort. Since the BCH series does not need the knowledge of the eigenvalues of the matrix obtained from the discretized model, it has an advantage over the CH procedure for large systems. In this respect let us mention a short “M ATHEMATICA” code [8] able to calculate any order of the BCH formula. The BCH method should be compared, however, with alternative methods of calculation of the monodromy matrix, which go beyond the two-step models, such as described in [9] and references therein.
Acknowledgements Part of this work was conducted during a visit of one of the authors (G.N.T.) to the Max-Planck Institut für Plasmaphysik, Garching. The hospitality of that Institute is greatly appreciated. The present work was performed
312
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
under the Contract of Association ERB 5005 CT 99 0100 between the European Atomic Energy Community and the Hellenic Republic.
References [1] [2] [3] [4] [5] [6]
H. Tasso, G.N. Throumoulopoulos, Phys. Plasmas 9 (2002) 2662. H. Tasso, P.P.J.M. Schram, Nucl. Fusion 6 (1966) 284. Yu.F. Golubev, J. Appl. Math. Mech. 63 (1999) 197. I.E. Leonard, SIAM Review 38 (1996) 507. G.H. Weiss, A.A. Maradudin, J. Math. Phys. 3 (1962) 771. J.G. Belinfante, B. Kolman, A Survey of Lie Groups and Lie Algebras with Applications and Computational Methods, SIAM, Philadelphia, 1972, p. 48. [7] M. Suzuki, Commun. Math. Phys. 57 (1977) 193. [8] M.W. Reinsch, J. Math. Phys. 41 (2000) 2434. [9] A.P. Seyranian, F. Solem, P. Pedersen, J. Sound, Vibration 229 (2000) 89.
Two-step control of wall mode and the monodromy matrix ✩ H. Tasso a,∗ , G.N. Throumoulopoulos a,b a Max-Planck-Institut für Plasmaphysik, Euratom association, D-85748 Garching bei München, München, Germany b University of Ioannina, Association Euratom-Hellenic Republic, Department of Physics, Section of Theoretical Physics, GR 451 10 Ioannina, Greece 1
Received 21 October 2002; received in revised form 10 December 2002; accepted 12 December 2002 Communicated by F. Porcelli
Abstract The Floquet stability of systems of differential equations with piecewise constant periodic coefficients is considered. In the “two-step” case the monodromy matrix is the product of two matrix exponentials. It can be evaluated either by reducing the matrix exponentials to polynomials on the eigenvalues of the corresponding matrices or, without knowledge of the eigenvalues, by using the Baker–Campbell–Hausdorff formula. The two methods are discussed and applied to several examples including the “two-step” dissipative Hill’s equation introduced recently H. Tasso, G.N. Throumoulopoulos [Phys. Plasmas 9 (2002) 2662] to stabilize the “resistive wall mode” in magnetically confined plasmas. Application to large systems and matrices is not only possible but needed for the full understanding of the dynamic stabilization of the resistive wall mode. Its practical implementation is discussed. 2002 Elsevier Science B.V. All rights reserved.
1. Introduction In a recent publication [1] we proposed to stabilize the “resistive wall mode” by modulating in time the resistivity of the wall. An actual implementation of this modulation can be produced by, e.g., a second nearby wall having cuts and electronic switches between the sides of the cuts as explained in the last paragraph of Section 4 of [1]. To study the stabilizing effect of this modulation we introduced the concept of an “effective potential” depending upon the instantaneous value of the resistivity. This led us to concentrate on the rather simple dissipative Mathieu–Hill equations, which resulted in very encouraging results related to the modulation strength of wall resistivity and the dissipation in the plasma. The physical situation can be well understood by assuming viscous magnetohydrodynamics in the plasma, then a vacuum region inside a resistive wall with Ohm’s law, and a vacuum region outside the wall extending to infinity or to a hypothetic perfectly conducting wall at a large but finite distance. The equations for linearized perturbations ✩
Presented at the 44th Annual Meeting of the Division of Plasma Physics of the APS, 11–15 November 2002, Orlando, FL, USA.
* Corresponding author.
E-mail address: [email protected] (H. Tasso). 1 Permanent address.
0375-9601/02/$ – see front matter 2002 Elsevier Science B.V. All rights reserved. doi:10.1016/S0375-9601(02)01769-3
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
can be written in a compact way as follows ξ˙ −µ∇ × ∇× −FMHD ∂ ξ = 1 0 ∂t 0 β A
ξ˙ α ξ , 0 −η(t)∇ × ∇× A
305
(1)
are the Lagrangian displacement vector inside the plasma and the potential vector inside the wall, where ξ and A respectively. µ is the plasma viscosity taken here as a scalar, η(t) is the time dependent resistivity of the wall. FMHD is the magnetohydrodynamic stability operator, α and β are nonlocal operators connecting the displacement vector at the plasma boundary with the vector potential at the surface of the wall, the gauge being chosen as the one of zero electrostatic potential. For details concerning FMHD and viscosity see, e.g., [2]. It is obvious that the Floquet stability of Eq. (1) is not easy to solve in real geometry of magnetic confinement devices like tokamaks, reversed field pinches or stellarators. As suggested in [1] this is a formidable task, which should be preceded by a through investigation of simple and relevant models. Discretization of Eq. (1) allows, however, to define a monodromy matrix as introduced in the next section. Here, preliminary investigation of the monodromy of simplified models, more advanced than those of Ref. [1], is carried through, especially for a “twostep” dependence of the modulated parameter.
2. Monodromy matrix Let Y be a n-vector and A(t) a n × n matrix depending periodically upon time in a piecewise constant manner A = A1 ,
0 < t < ts ,
(2)
A = A2 ,
ts < t < T ,
(3)
where ts is the time at which a jump in A occurs, and T is the period. A1 and A2 are time independent n × n matrices. Define now the equation with piecewise constant periodic coefficients as ˙ = A(t)Y. Y
(4)
The solution of (4) can be written, in general, as Y(t) = M(t)Y(0),
M(0) = I
(5)
with the matrix M(t) obeying ˙ = A(t)M. M
(6)
For the two-step case defined by (2) and (3) we have M(t) = exp tA1 , 0 < t < ts , M(t) = exp (t − ts )A2 exp(ts A1 ),
(7) ts < t < T .
(8)
Since the monodromy matrix is defined by M(T), it can be written [3] after integrating (6) over a period as M(T ) = exp (T − ts )A2 exp(ts , A1 ). (9) Since A1 and A2 do not commute, in general, it is not easy to find explicitly the product of two matrix exponentials, especially for matrices of large order. In case the eigenvalues of A1 and A2 are known, it is possible to replace each exponential by a matrix polynomial on the eigenvalues of the corresponding matrix by applying the Cayley–Hamilton (CH) theorem as done implicitly in [3] or, more generally, as explained in [4]. This method
306
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
gives exact results at least for low order ( 4) matrices, though the calculations may become cumbersome and then should be done symbolically by computer. Finally, let us recall that the system is stable if and only if all the eigenvalues of the monodromy matrix (9) have an absolute value 1.
3. Baker–Campbell–Hausdorff formula For large matrices, for which the eigenvalues are not known explicitly, it may be convenient to use the Baker– Campbell–Hausdorff (BCH) formula (see, e.g., [5–7]), which states that M(T ) = exp C
(10)
with
ts (T − ts ) t 2 (T − ts ) [A2 , A1 ] + s [A2 , A1 ], A1 2 12
ts (T − ts )2 + [A1 , A2 ], A2 + · · · , 12 which, in the “symmetric” case ts = T /2, becomes C = (T − ts )A2 + ts A1 +
(11)
T T2 T 3 (12) (A1 + A2 ) + [A2 , A1 ] + [A2 , A2 ], A1 + [A1 , A2 ], A2 + · · · , 2 8 96 where [x, y] denotes the commutator of x and y. The matrices A1 and A2 are identical up to entries depending of the two-step control function. For a vanishing jump the two matrices commute, which means that their commutator must be a polynomial in the jump η2 − η1 . Since C is the log of the monodromy matrix M(T ), the stability of the system demands that the real parts of the eigenvalues of the matrix C be negative or zero. It is difficult to answer this question for finite jumps since C is given by an infinite series. Throughout this Letter the BCH formula is, however, taken as in (11) and (12), i.e., to the third order. C=
4. Two-step dissipative Hill equation For small jumps we keep in (12) the three first terms, for an example, the “two-step” dissipative Hill equation [1] x¨ + d x˙ + b(t)x = 0, or in the form −d y˙ = 1 with
−b(t) y, 0
x˙ y= . x
A1 and A2 are given by −d −b1 A1 = , 1 0 −d −b2 A2 = . 1 0
(13)
(14)
(15)
(16) (17)
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
Inserting (16) and (17) in formula (12) we obtain 2 2 2 + d T8 (b1 − b2 ) − −dT + T8 (b1 − b2 ) −T b1 +b 2 C= 2 T − T8 (b1 − b2 )
T3 48 (b1
− b 2 )2
307
.
(18)
The eigenvalues of C as given by (18) are
4 −dT ± d 2 T 2 − 2T 2 (b1 + b2 ) − T48 (b1 − b2 )2 λ= (19) . 2 We see that, if (b1 + b2 ) > 0, the system is stable with and without dissipation. This shows that the corresponding unstable regions found in the charts displayed in Figs. 3 and 5 of Ref. [1] need higher order terms in the BCH expansion for their description. For (b1 + b2 ) < 0, however, we can find the “parabolic” beginning near the origin for the lowest curve of the same chart by requiring that the real part of the unstable eigenvalue be zero as follows (b1 + b2 ) T 2 (b1 − b2 )2 = . (20) 2 48 4 Let us notice that (b1 + b2 )/2 is to be set equal to K of Ref. [1] and (b1 − b2 )/2 should be set equal to 2k of Ref. [1], T being equal to 2π/Ω. It will not be easy to obtain the full charts of Ref. [1] since we should be able to take into account the higher order commutators in the BCH formula either by symbolic computation or by summing up the series. −
5. Asymmetric case or ts = T /2 The extension of the previous calculation to the time-asymmetric case can be done on the same line. Since the algebra is rather cumbersome the elements of matrix C and its eigenvalues were obtained via symbolic computation. Explicit expressions are given below. −12 dT − (b1 − b2 )(T − ts )ts (−6 + dT − 2 dts ) , 12 1 −12b1ts − (T − ts ) 12b2 + (−b1 + b2 ) d(6 − dT ) + 2T b2 ts + 2(b1 − b2 ) −d 2 + b1 + b2 ts2 , C12 = 12 (b1 − b2 )(T − 2ts )(T − ts )ts , C21 = T + 6 (b1 − b2 )(T − ts )ts (−6 + dT − 2dts ) C22 = 12 C11 =
and λ=
√ 1 6 dT ± ∆ , 12
where ∆ = 36T 2 d 2 − 4b2 + 144T (−b1 + b2 )ts + T 2 (b1 − b2 )2 −12 − d 2 T 2 + 4T 2 b2 ts2 − 2T (b1 − b2 )2 −3 4 + d 2 T 2 + 2T 2 (b1 + 5b2) ts3 + (b1 − b2 )2 −12 − 13 d 2T 2 + 4T 2 (4b1 + 9b2) ts4 + 4T 3d 2 − 5b1 − 7b2 (b1 − b2 )2 ts5 − 4 d 2 − 2b1 − 2b2 (b1 − b2 )2 ts6 .
308
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
Fig. 1. The lowest marginal curves for two-step Hill’s system obtained analytically by means of the BCH formula (11) and numerically on the basis of the CH method associated with the monodromy matrix (9). (a) Shows curves for symmetric in time modulations (ts = 1) while (b) and (c) are associated with asymmetric modulations, i.e., ts = 1.3 and ts = 0.7, respectively.
Introducing physically plausible parameters, i.e., the modulation average b1 ts + b2 (T − ts ) T and the deviation of b1 from K ts b1 ts + b2 (T − ts ) ts , k = (b1 − K) = b1 − T T T K=
and requiring that the real part of the unstable eigenvalue be zero yields K=
k 2 T 2 (12 + (d 2 + 8k)T 2 − 4(d 2 + 4k)T ts + 4d 2 ts2 ) . 4(−6 + kT 2 − 2kT ts )(6 + kT 2 − 2kT ts )
(21)
The lowest marginal curves determined by (21) are shown in Fig. 1 along with the respective exact ones. The latter were obtained numerically on the basis of the CH method; in particular, to determine the marginal curves the eigenvalues λ of the monodromy matrix (9) were calculated numerically and then the marginal b values via the relation |λ| = 1. The results were checked by means of the dispersion relation (17) of Ref. [1]. BCH and CH symmetric curves are also presented in Fig. 1(a) which compare well each other for values k < 1 but they deviate for larger ones. The deviation is affected rather strongly by the asymmetry as can be seen in Fig. 1(b) and (c) thus indicating that higher order terms are required in the BCH expansion in particular for large asymmetric modulations.
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
309
6. Jumps on the diagonal For diagonal matrices A1 and A2 all the commutators in (11) and (12) vanish, which means that there is no proper “dynamic stabilization” in contrast to the case of the Hill’s equation. For nondiagonal matrices having jumps on the diagonal, it may be expected that the jumps have to be large enough to be able to stabilize. It is instructive, as a crude model of Eq. (1), to consider the following problem: −d −b α y, y˙ = 1 0 0 (22) β1 β2 −η(t) with x˙ y= x . (23) u A1 and A2 are given by −d −b α 0 0 , A1 = 1 (24) β1 β2 −η1 −d −b α 0 0 . A2 = 1 (25) β1 β2 −η2 As mentioned above the jump is now in η, which is the third element of the diagonal. Let us now calculate C as given by Eq. (12) T 2 (η2 −η1 ) T 3 (η2 −η1 )2 C=
−dT
−bT
T
0
αT +
0
+
T 2 (η2 −η1 ) T 3 (η2 −η1 )2 T β1 − β1 + β1 8 96
T 2 (η2 −η1 ) T 3 (η2 −η1 )2 T β2 − β2 + β2 8 96
η +η − 12 2
8
96
. (26)
The eigenvalues of matrix (26) are given by a cubic equation, in which the jump appears with even powers, the coefficients of the second power being the strongest as displayed in Eq. (27).
T (η1 + η2 ) λ(λ + dT ) + bT 2 λ+ 2 αβ1 T 2 (η2 − η1 )2 αβ2 T 2 (η2 − η1 )2 − λ αβ1 T 2 + (27) + · · · − T αβ2 T 2 + + · · · = 0. 192 192 Eq. (27) shows that the stability will not be affected by small or moderately large jumps in contrast with the “two-step” Hill’s equation. Marginal curves for symmetric modulations respective to the lowest two-step Hill ones (which can be recovered for α = 0) were also obtained numerically for β1 = 0 and d = 1. Inspection of (27) and the numerical results reveal the following characteristics of the marginal curves: (1) The marginal points are determined by given product K = αβ2 regardless of the particular values of α and β2 . (2) For a “perfectly conducting wall”, η = 0, the marginal b take larger values as either α or β2 increase. (3) The larger the resistivity η1 = η2 ≡ η the smaller the marginal b as can be seen in Fig. 2(a) and (b) for η ≡ η2 − η1 = 0, i.e., for η = 0.6 and η = 1 the respective marginal b values are nearly 0 and −1. It is also noted that for η = 0 the marginal b is 1. (4) There are parametric regimes for which the marginal b take negative values and decrease with increasing η (Fig. 2(a)).
310
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
Fig. 2. BCH and CH marginal curves for a wide variation of the symmetric in time resistivity modulation parameter ∆η = η2 − η1 and different combinations of η1 and the wall parameters α and β2 . The system is stable just above the curves and unstable below.
(5) The monotonicity of the curves depends on the sign of K as follows: • For K < 0 the curves can vary nonmonotonically as η increases. For example, as shown in Fig. 1(a) there is a sharp turning point at η ≈ 1 at which the marginal (b) is nearly −1. For α = −1 and β2 = 1 such a turning point appears only for η1 < 1. For η1 1 the marginal b in the CH case increases monotonically with η remaining, however, negative. The particular case of η1 = 1 is shown in Fig. 2(b). • For K > 0 the CH curves go down monotonically as η increases (Fig. 2(c)). Even for small values of the parameters α and β2 , however, the curves do not involve negative values of b. (6) Similar to the case of the two-step Hill’s system the BCH curves are reliable for η 1 while for larger η they deviate significantly from the CH ones. In order to have a proper interpretation of the above results, let us first state the main features of the actual resistive wall mode described by the full physical equation (1): (1) The mode is an “external kink” with a reduced growth rate due to the resistivity of the wall. (2) The mode should be stable if the wall is perfectly conducting (η = 0) and the wall is not too far from the plasma boundary. (3) The mode should become unstable if the same wall is resistive (η = 0). Resistivity makes the wall “permeable” to the perturbed magnetic field in the full vacuum in whole space or up to a perfectly conducting second wall at large distances. For the 3 by 3 model the parameters α and β2 stand for the importance of the vacuum region while the parameter b represents the perturbation of plasma potential energy, the wall mode being associated with negative
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
311
values of b. If α = 0 and η = 0, the system is marginal for b = 0 (see Fig. 5 of Ref. [1]). Taking α = −1 and keeping zero resistivity makes the system marginal at b = 1. This is representative of the external kink which becomes unstable as soon as the vacuum space is increased which has to be compensated by a higher value of b. So the features (1) and (2) are qualitatively supported by the model. The calculations show, however, that for η1 = η2 = 0.6 the value of b at marginality decreases from 1 to almost zero and for η1 = η2 = 1.0 it decreases further to −1. This would suggest that resistivity of the wall helps stabilizing. In fact, this would be plausible for a wall of zero thickness (“skin” thickness is larger than wall thickness). This may explain also why the agreement between CH and BCH curves in Fig. 2(a) is so good even for η/η1 = 1. At this stage the dynamic stabilization is not visible; this was confirmed by taking in the BCH formula the first term only (no commutators). The existence of a turning point is particularly interesting since the marginal curve has a behaviour similar to what one would expect for dynamic stabilization though quite crude. Indeed, if for example, b is fixed at −0.5 in Fig. 2(a) and η is increased from zero to larger values the system becomes first stable and then for η increasing further becomes unstable again. This behaviour ceases, however, for b < −1. In contrast, the monotonic curves of Fig. 2(b) and (c) do not indicate that a dynamic stabilization of the mode is taking place. The only way to cure the model is to go to larger matrices. This may, of course, be limited by numerical feasibility. The most significant progress is that the method proposed seems to work and could, in principle, be extended in the future to larger problems.
7. Conclusions The CH and BCH methods proposed here to handle Floquet stability of linear systems of differential equations with periodic piecewise constant coefficients are very powerful as demonstrated in the previous sections. For instance, the first marginal curve for “negative energy” modes of the “two-step” Hill equation (13) [1] can be partly reproduced analytically. As mentioned in the introduction, it would be good to assess the results of this model by a more accurate treatment of the dynamic stabilization of the resistive wall mode through time-modulation of the wall resistivity as stated in system (1). In principle, CH and BCH methods can be applied to large systems of differential equations. As an example, the 3 × 3 system (21), a crude mock up of system (1), has been successfully treated in Section 6. The discussion there shows, however, that model (21) is too crude to give all qualitative features of the resistive wall mode. Though dynamical stabilization does not show up as expected because the “wall” has, so to say, zero thickness, the results are very encouraging. This means that we have to go to much larger systems approaching more realistically the plasma and a resistive wall having a thickness larger or comparable to the “skin” thickness. An accurate discretization of system (1) with subsequent application of CH or BCH methods would lead to much higher numerical effort. Since the BCH series does not need the knowledge of the eigenvalues of the matrix obtained from the discretized model, it has an advantage over the CH procedure for large systems. In this respect let us mention a short “M ATHEMATICA” code [8] able to calculate any order of the BCH formula. The BCH method should be compared, however, with alternative methods of calculation of the monodromy matrix, which go beyond the two-step models, such as described in [9] and references therein.
Acknowledgements Part of this work was conducted during a visit of one of the authors (G.N.T.) to the Max-Planck Institut für Plasmaphysik, Garching. The hospitality of that Institute is greatly appreciated. The present work was performed
312
H. Tasso, G.N. Throumoulopoulos / Physics Letters A 307 (2003) 304–312
under the Contract of Association ERB 5005 CT 99 0100 between the European Atomic Energy Community and the Hellenic Republic.
References [1] [2] [3] [4] [5] [6]
H. Tasso, G.N. Throumoulopoulos, Phys. Plasmas 9 (2002) 2662. H. Tasso, P.P.J.M. Schram, Nucl. Fusion 6 (1966) 284. Yu.F. Golubev, J. Appl. Math. Mech. 63 (1999) 197. I.E. Leonard, SIAM Review 38 (1996) 507. G.H. Weiss, A.A. Maradudin, J. Math. Phys. 3 (1962) 771. J.G. Belinfante, B. Kolman, A Survey of Lie Groups and Lie Algebras with Applications and Computational Methods, SIAM, Philadelphia, 1972, p. 48. [7] M. Suzuki, Commun. Math. Phys. 57 (1977) 193. [8] M.W. Reinsch, J. Math. Phys. 41 (2000) 2434. [9] A.P. Seyranian, F. Solem, P. Pedersen, J. Sound, Vibration 229 (2000) 89.