Quasi-optical theory of microwave plasma heating in open magnetic trap
A. G. Shalashov, A. A. Balakin, E. D. Gospodchikov, T. A. Khusainov
QQuasi-optical theory of microwave plasma heating in open magnetic trap
A. G. Shalashov,
1, 2, a) A. A. Balakin,
1, 2
E. D. Gospodchikov,
1, 2 and T. A. Khusainov
1, 2 Institute of Applied Physics of the Russian Academy of Sciences, Ulyanova str. 46, 603950 Nizhny Novgorod,Russia Budker Institute of Nuclear Physics, Siberian Branch of the Russian Academy of Sciences,Akademika Lavrentieva ave. 11, 630090 Novosibirsk, Russia (Dated: July 31, 2018)
Microwave heating of a high-temperature plasma confined in a large-scale open magnetic trap, including allimportant wave effects like diffraction, absorption, dispersion and wave beam aberrations, is described for thefirst time within the first-principle technique based on consistent Maxwell’s equations. With this purpose, thequasi-optical approach is generalized over weakly inhomogeneous gyrotrotropic media with resonant absorp-tion and spatial dispersion, and a new form of the integral quasi-optical equation is proposed. An effectivenumerical technique for this equation’s solution is developed and realized in a new code
QOOT , which isverified with the simulations of realistic electron cyclotron heating scenarios at the Gas Dynamic Trap at theBudker Institute of Nuclear Physics (Novosibirsk, Russia).
I. INTRODUCTION
The absorption of electromagnetic waves under theelectron cyclotron resonance (ECR) conditions is widelyused for heating of high-temperature plasmas in toroidalmagnetic traps (tokamaks and stellarators). However,for many years the use of this method in open magneticconfigurations has been limited either by plasma heat-ing in relatively compact laboratory installations , or byMHD stabilization of low-density plasma . The onlyexception was the TMX-U experiment at Lawrence Liv-ermore National Laboratory where electron temperaturesup to 0.28 keV were obtained with ECR heating, butthese studies were concluded soon . And only recently anefficient ECR heating of a dense (comparable to toroidaldevices) plasma has been demonstrated at a large-scalemirror trap. We mean successful experiments on a com-bined plasma heating by neutral beams and microwaveradiation performed at the Gas Dynamic Trap (GDT)device in the Budker Institute of Nuclear Physics . Inparticular, direct and highly-localized heating of the ther-mal electron component allows to achieve the GDT ex-periment the record for open traps electron temperatureof 1 keV at the plasma axis . These studies convincinglydemonstrate good prospects for the use of simple axi-ally symmetric open magnetic traps as a powerful sourceneutron for fusion applications .Implementing the effective ECR heating of a denseplasma at a large open trap requires a revision of theprevailing ideas about the physics of cyclotron absorp-tion as well as the subsequent transport of energy andMHD stabilization of a plasma column. None of the well-understood techniques developed for the toroidal plasmaheating works well in this case . Numerical modelingof the propagation and absorption of electromagneticwaves in an inhomogeneous plasma plays an important a) Author to whom correspondence should be addressed. Electronicmail: [email protected] role in this occasion. Until recently, such modeling wasonly possible in the framework of the geometric opticsapproximation also known as ray-tracing. In particu-lar, the basic ECR heating scenario used in the GDThas been originally proposed and justified with the ray-tracing calculations . However, in this scenario thereare areas of reflection and strong absorption of waves,in which the medium is not smooth as compared to awavelength. So, detailed understanding of the involvedprocesses requires to go beyond the geometric optics ap-proximation.The main effects that violate the geometric optics areassociated with the spatial dispersion in a strongly inho-mogeneous region of resonant absorption, the diffractionof a wave beam, and the caustic formation in the vicinityof internal reflection points of a wave beam. Straightfor-ward simulation of these effects for large devices within acomplete set of Maxwell’s equations is very complicated,in particular, because of the smallness of a wavelength. Agood alternative is the consistent quasi-optical approachbased on an asymptotic expansion of Maxwell’s equa-tions in the paraxial approximation in the vicinity of theselected Wentzel–Kramers–Brillouin (WKB) mode .In this paper, the quasi-optical approximation isadopted to describe the propagation of wave beams ina high-temperature plasma in an open magnetic trap. Asimilar approach was previously developed for toroidalmagnetic traps . However, description of the mi-crowave heating in modern open traps needs substantialmodification of the quasi-optical theory due to a require-ment of more accurate accounting of the spatial disper-sion inside the resonant absorption zone. The physicsof this difference originates from the fact that the mag-netic field in open configuration typically changes alongits direction, while in the toroidal geometry the magneticfield varies presumably in orthogonal direction. As a re-sult, the factorizing approach to the quasi-optical Hamil-tonian, that works well in the modeling of wave prop-agation in a toroidal trap, fails completely in case of alarge open trap. This motivates us to find a new, moregeneral and accurate, non-factorized formulation of the a r X i v : . [ phy s i c s . p l a s m - ph ] N ov theory resulted in the substantial modification of a nu-merical algorithm, which eventually led to developmentof an entirely new code.The remainder of this paper is organized as follows:Sec. II derives the basic equations of quasi-optical ap-proximation and their applicability conditions; Sec. IIIdefines a general ansatz for a dielectric response of inho-mogeneous dispersive media and derive the quasi-opticalevolution operator; peculiarities related to wave dissipa-tion are considered, and modification of the dissipativepart of the early defined evolution operator is proposed inSec. IV; a numerical method for solving the quasi-opticalequation is developed in Sec. V; the model is adopted to aparticular case of hot magnetized plasmas confined in anopen magnetic trap in Sec. VI; first results of simulationsof real scenarios of ECR plasma heating at the GDT de-vice are discussed in Sec. VII; and Sec. VIII summarizesthe results. II. BASIC QUASI-OPTICAL EQUATION
A rigorous derivation of the quasi-optical approx-imation for a weakly inhomogeneous gyrotropic andanisotropic media with spatial dispersion may be foundin Refs. 14 and 18. The main steps are described below.Let us start with Maxwell’s equations for e − i ωt processeswritten in the “operator” form ˆ L ij [ E j ( r )] = 0 , ˆ L ij ≡ δ ij ˆ k − ˆ k i ˆ k j − k ˆ ε ij , (1)where a summation over double indexes is implied. Theoperator ˆ L ij [ E j ( r )] acts on j -th Cartesian component E j ( r ) of the electric field vector, k = ω/c is the vacuumwave number, ˆ k is the wave number operator definedas the differentiation in the coordinate space, and ˆ ε ij isthe linear operator obtained from the dielectric tensor ε ij ( r , k ) by a formal substitution k → ˆ k , ˆ k = − i ∇ , ˆ ε ij = ε ij ( r , ˆ k ) . The latter operation is not uniquely defined because r and ˆ k do not commute; in the next section we describethe “best way” of regularizing the expression for ˆ ε ij pro-viding that the kernel ε ij is known. The result may beexpressed in terms of Fourier-integral: ˆ ε ij [ E j ( r )] = (cid:90) ε ij ( r , k (cid:48) ) + ε ij ( r (cid:48) , k (cid:48) )2 ×× e i k (cid:48) ( r − r (cid:48) ) E j ( r (cid:48) ) d r (cid:48) d k (cid:48) (2 π ) . (2)For the smoothly inhomogeneous medium, the approx-imate solution of Maxwell’s equations can naturally befound in the form of quasi-optical beam correspondingto some chosen WKB mode. To do this, we define thepolarization vector e for the interested WKB mode as asolution of algebraic Maxwell’s equations L ij e j = 0 in a locally homogeneous medium. Next, using the sametechnique as in the transition ε ij → ˆ ε ij , the polarizationvector e ( r , k ) can be associated with the polarization op-erator ˆ e ( r , ˆ k ) . Then, the approximate solution of Eq. (1)can be found in the following form E j ( r ) = ˆ e j ( r , ˆ k )[ U ( r )] . Here, the polarization operator acts on the scalar ampli-tude U ( r ) of the wave beam. The equation for U takesthe form ˆ H ( r , ˆ k )[ U ( r )] = 0 , ˆ H ≡ ˆ e ∗ i ˆ L ij ˆ e j . (3)Note that the corresponding function H GO = R e H ( r , k ) is usually used as a ray Hamiltonian in the geometricaloptics. Similar, we will call the operator ˆ H as a quasi-optical Hamiltonian.Presently, we consider axisymmetric open traps andmagnetic configurations with weakly broken axial sym-metry. In this case, we can use either Cartesian coordi-nate system or a cylindrical coordinate system with thepolar axis aligned along the trap axis. Let us introducethe coordinates and the canonically conjugated momentaas r ≡ ( x , z ) , ˆ k ≡ (ˆ q , ˆ k z ) , (4)where z is the coordinate along the trap axis, x = r ⊥ isthe set of two coordinates in the transverse plane, ˆ q = − i ∇ ⊥ and ˆ k z = − i ∂/∂z .Assuming the smoothness of beam parameters along z , we can introduce the scalar amplitude of U ( r ) in theform U ( r ) = ˜ u ( x , z ) exp (cid:18) i (cid:90) κ ( z ) dz (cid:19) , where ˜ u is the complex envelope of a wave beam, andfunction κ ( z ) defines the dependence of the “carrierphase” of the wave field along the axis. Substituting thisfield in Eq. (3), we obtain the equation for the envelope ˜ u as ˆ H (cid:16) r , ˆ k + κ ( z ) z (cid:17) [˜ u ( x , z )] = 0 . (5)In contrast to geometrical optics, here we have certainfreedom in the choice of κ ( z ) related to the phase varia-tion over the transverse aperture of the wave beam. Afterseveral tries, we found that the most convenient and re-liable way in numerical calculations is to define κ ( z ) asa solution to a local dispersion equation with transversewave vector ˜q ( z ) corresponding to the center of mass ofa wave spectrum in the cross-section z : H GO ( ˜r , ˜q ( z ) + κ ( z ) z ) = 0 , ˜r = ( ˜x , z ) , where ˜x is the center of mass for a wave field in the cross-section z .Assuming the complex envelope ˜ u to be a smooth func-tion on the longitudinal coordinate, | ∂ ˜ u/∂z | (cid:28) κ | ˜ u | , weexpand Eq. (5) to the first order in powers of ˆ k z (or,equivalently, in powers of ∂ ˜ u/∂z ): ˆ H [˜ u ] ≈ ˆ H [˜ u ] − i ∂ ˆ H ∂κ (cid:20) ∂ ˜ u∂z (cid:21) − i2 ∂ ˆ H ∂κ∂z [˜ u ] = 0 , (6)where ˆ H = ˆ H ( r , ˆ q + κ z ) , i.e. it is ˆ H with omittedderivatives over z . With the limit of geometrical optics,the last term in Eq. (6) results in a pre-exponential factorresponsible for the variation of the wave amplitude duethe effect of “group” slowing down. This term may beeliminated with formal substitution ˜ u = ˆ A [ u ] , ∂ ˆ H ∂κ ∂ ˆ A∂z + 12 ∂ ˆ H ∂κ∂z ˆ A = 0 . As a result, we obtain a more simple equation ∂u∂z = i k ˆ H [ u ] , ˆ H = − (cid:18) k ∂ ˆ H ∂κ ˆ A (cid:19) − ˆ H ˆ A. (7)This is our basic equation that describes the evolutionof the scalar wave beam amplitude u along z -axis, in-cluding the effects of diffraction, spatial dispersion anddissipation, as long as condition | ∂u/∂z | (cid:28) κ | u | of theparaxial approximation holds. Except few very specialcases, operator ∂ ˆ H /∂κ ≈ ∂ H GO /∂κ is local, i.e. it maybe approximated with multiplication by known functionas long as the media is smooth compared to the wavelength. In this case, the new amplitude is also defined bylocal transformation u = ˜ u (cid:112) k ∂ H GO /∂κ .The simple form of Eq. (7) suggests a clear physicalinterpretation of the evolution operator ˆ H . Indeed, thisis a longitudinal wave number expressed as a function ofthe transverse wave vector in an operator form: k ˆ H = k z ( r , − i ∇ ⊥ ) − κ ( z ) . (8)Below we will show how to restore this operator fromthe WKB-solution k z ( r , k ⊥ ) of an algebraic dispersionrelation in a locally homogeneous medium.Our approach may be adopted to highly curved mag-netic configurations, e.g. typical for radiation belts inEarth ionosphere or Solar flares, in a straightforwardmanner—just by introducing new coordinates with acurvilinear axis z and taking into account the curvaturewhen calculating a conjugated momenta operator. Thesimilar approach was previously used for toroidal mag-netic traps in which the axis z was chosen along the ref-erence geometric optics ray representing the center of thequasi-optical wave beam . However, using geometricoptics rays as a reference for the quasi-optical equationin open traps is not optimal because such rays may bestrongly curved and even divergent inside the plasma col-umn in most interesting cases . III. EVOLUTION OPERATOR
The main difficulty in the practical use of the quasi-optical equation (7) is the uncertainty of the dielectricpermittivity operator ˆ ε , which is necessary for definitionof the evolution operator ˆ H . Obviously, the operator ˆ ε for any inhomogeneous linear media can be generated bysome kernel function ε ij ( r , k ) , and such kernel can beintroduced in infinite different ways. The most optimal(for our needs) way is proposed in quite rigorous mannerin Ref. 19, where the dielectric operator is represented as ˆ ε ij = ε ij + 12 (cid:18) ∂ε ij ∂k α ∆ˆ k α + ∆ˆ k α ∂ε ij ∂k α (cid:19) ++ 12 12! (cid:18) ∂ ε ij ∂k α ∂k β ∆ˆ k α ∆ˆ k β + ∆ˆ k α ∆ˆ k β ∂ ε ij ∂k α ∂k β (cid:19) + . . . == 12 ∞ (cid:88) n =0 n ! (cid:18) ∂ n ε ij ∂ k n ∆ˆ k n + ∆ˆ k n ∂ n ε ij ∂ k n (cid:19) . (9)Here ∆ˆ k = ˆ k − k = − i ∇ − k . For such operator, aHermitian kernel, ε ij = ε ∗ ji , always generates a Her-mitian operator ˆ ε in a sense of natural scalar product ( a, b ) = (cid:82) ab ∗ d r . Correspondingly, an anti-Hermitiankernel, ε ij = − ε ∗ ji , generates an anti-Hermitian operator.This important non-trivial feature allows to preserve theenergy conservation law when constructing the dielectricresponse of inhomogeneous media. The kernel function ε ij is obtained as the sum of the amendments to thedielectric permittivity tensor over the powers of charac-teristic inhomogeneity scale l : ε ij = ∞ (cid:88) m =0 ε ( m ) ij , ε ( m ) ij = O (( k l ) − m ) . As a reasonable approximation in a smoothly inhomo-geneous media with k l (cid:29) , one can only account thelowest term, ε ij ≈ ε (0) ij , which is indeed the dielectricpermittivity tensor obtained for a “locally homogeneous”media, i.e. in the geometric optics approximation.To find a Fourier representation, let us substitute E j ( r ) = (cid:90) e i k (cid:48) ( r − r (cid:48) ) E j ( r (cid:48) ) d r (cid:48) d k (cid:48) (2 π ) ,E j ( r ) ∂ n ε ij ( r , k ) ∂ k n = (cid:90) e i k (cid:48) ( r − r (cid:48) ) E j ( r (cid:48) ) ∂ n ε ij ( r (cid:48) , k ) ∂ k n d r (cid:48) d k (cid:48) (2 π ) . into Eq. (9). Then ˆ ε ij [ E j ( r )] = (cid:90) ∞ (cid:88) n =0 n ! ∂ n ∂ k n ( ε ij ( r , k ) + ε ij ( r (cid:48) , k )) ×× ( k (cid:48) − k ) n e i k (cid:48) ( r − r (cid:48) ) E j ( r (cid:48) ) d r (cid:48) d k (cid:48) (2 π ) . Noting that summation over n is just a Taylor seriesfor ε ij ( . . . , k (cid:48) ) we obtain Eq. (2). It should be stressedthat expressions (2) and (9) remain valid for any ker-nel ε ij ( r , k ) , and for every linear inhomogeneous problemsuch kernel exists. The problem is, however, that exceptsome trivial model cases, see e.g. Ref. 19, this kernelis known only with the accuracy provided by geometricoptics.The dielectric permittivity operator unambiguouslydefines the evolution operator in Eq. (7). However, wecan cut the computations and obtain the same result byjust noting that the evolution operator may be found es-sentially in the same way as the dielectric operator ˆ ε .So, the operator ˆ H is generated by the scalar function H ( x , q )ˆ H [ u ] = (cid:90) H ( x , q (cid:48) ) + H ( x (cid:48) , q (cid:48) )2 e i q (cid:48) ( x − x (cid:48) ) u ( x (cid:48) ) d x (cid:48) d q (cid:48) (2 π ) , (10)and as a practical approximation to H ( x , q ) in a weaklyinhomogeneous media, one can use a solution of thelocally homogeneous problem. The only difference weshould have in mind is that coordinate z is a param-eter here—there is no differentiation over z in opera-tor ˆ H , therefore instead of 6-dimensional phase-space d r (cid:48) d k (cid:48) / (2 π ) we use 4-dimensional space d x (cid:48) d q (cid:48) / (2 π ) .Due to the term exp(i q (cid:48) ( x − x (cid:48) )) , the evolutionary opera-tor (10) becomes local either in real space for an inhomo-geneous medium without spatial dispersion with H ( x ) ,or in q -space for a homogeneous dispersive medium with H ( q ) , or for combination H ( x ) + H ( q ) . The latter hasbeen used in simulations for toroidal configurations ,but we find this approximation being not suitable foropen traps. Thus, we must deal with a general form H ( x , q ) below.As already mentioned, symmetric form of the integralrepresentation allows to introduce the dielectric responseof an inhomogeneous medium that is compatible withthe energy conservation. In the present context, Eq. (10)always results in Hermitian operators for real kernel H and anti-Hermitian operators for purely imaginary H .However, this feature is insufficient to provide a correctdescription of inhomogeneous dissipative medium withspatial dispersion. As shown in the next section, equa-tion (10) may give absurd results and, therefore, needsmodification. IV. DISSIPATIVE CORRECTION
Let us decompose the linear operator ˆ H over Hermitianand anti-Hermitian parts ˆ H = ˆ H H + i ˆ H A , where ˆ H H ≡ ( ˆ H + ˆ H † ) / and ˆ H A ≡ ( ˆ H − ˆ H † ) / are bothHermitian operators in the sense of the standard scalarproduct (cid:82) ab ∗ d x . It can be shown that with the sameaccuracy as Eq. (7), (cid:82) | u | d x is proportional to the totalwave energy flux in the cross-section z . Then, using the quasi-optical equation (7) one can find the linear densityof the absorbed power as P ≡ − ∂∂z (cid:90) | u | d x = 2 k (cid:90) u ∗ ˆ H A [ u ] d x . (11)Hermitian operator in Eq. (7) conserves the energy flux,thus describing non-dissipative medium. As expected,the anti-Hermitian part of the evolution operator is re-sponsible for the dissipation or pumping of the wave field.From the geometric optics limit, we know thatenergy dissipation corresponds to a positive sign of I m H ( x , q ) . Let us now consider the eigenvalues λ A ofthe “dissipative” operator ˆ H A with the dissipative kernel I m H ( x , q ) > according to Eq. (10). One may expectthat all λ A > , however this is not true. Procedure (10)guarantees only λ A > since ˆ H A is Hermitian, but notthe positive sign of λ A . In particular, we had negative λ A with a first try to apply Eq. (10) to model the resonantwave absorption in magnetized plasma.Fortunately, we have found a simple practical way tocorrect the evolution operator, which guarantees a def-inite sign of λ A . To do this, we define ˆ H A a square ofsome other Hermitian operator: ˆ H A = ˆ G (cid:2) ˆ G [ u ] (cid:3) , (12)where operator ˆ G is constructed similar to Eq. (10) butwith kernel G ( x , q ) = (cid:112) I m H ( x , q )ˆ G [ u ] = (cid:90) G ( x , q (cid:48) ) + G ( x (cid:48) , q (cid:48) )2 e i k q (cid:48) ( x (cid:48) − x ) u ( x (cid:48) ) d x (cid:48) d q (cid:48) (2 π ) . New operator ˆ H A provides the same WKB limit as theold one. If the operator kernel is defined within the ge-ometric optics approximation, then formal difference be-tween both operators is of higher order than their accu-racy. The modified quasi-optical equation (7) then takesthe following form ∂u∂z = i k ˆ H H [ u ] − k ˆ G [ ˆ G [ u ]] . (13)The discussed problem of the absorbed power in aninhomogeneous medium with spatial dispersion may beillustrated with a simple one-dimensional example. Con-sider the following Hamiltonian kernel with positiveimaginary part H A ( x, k x ) = α x k x ≥ . Corresponding non-modified evolution operator (10) is ˆ H A = − α (cid:0) x ∂ xx + 2 x∂ x + 1 (cid:1) . According to Eq. (11), the absorbed power density alongthe z axis is P = 2 α k (cid:90) (cid:0) | x∂ x u | − | u | (cid:1) dx. The r.h.s. of this expression does not conserve the sign.For example, the Gaussian wave beam u = u e − ( x − a ) /w corresponds to P = (cid:112) π/ α u k (4 a /w − w ) . In particular, non-shifted beam with a = 0 , regardless ofits width w , always results in a negative absorbed power.Shifted beam with a = w/ propagates without dissi-pation. These strange solutions disappear for modifiedevolution operator (12). In this case ˆ G = i α (cid:0) x∂ x + (cid:1) , ˆ H A = − α (cid:0) x ∂ xx + 2 x∂ x + (cid:1) , and the absorbed power density of the new operator pos-sesses a positive sign: P = 2 α k (cid:90) (cid:0) | x∂ x u | − | u | (cid:1) dx == 2 α k (cid:90) (cid:12)(cid:12) x∂ x u + u (cid:12)(cid:12) dx > . Note that the dissipation-free solution is still present, butin this case it corresponds to unbounded distribution u ∝| x | − / with an infinite power flux. V. NUMERICAL SOLUTION
The linear integro-differential equation (7) involvestime consuming integration in 4-dimensional space d x (cid:48) d q (cid:48) . Therefore, an algorithm economic both in cal-culation of the r.h.s. and in number of z -steps is of ma-jor importance. We develop such algorithm having inmind essential features of a beam propagation physics inresonant dispersive media.As a starting point one can consider an explicit schemein which a small but finite integration step ∆ over theevolution coordinate z is defined as u ( x , z + ∆) ≈ ˆ S H [ u ( x , z )] , where ˆ S H [ u ( x , z )] ≡ (cid:90) exp (cid:18) i k ∆ H ( x , q (cid:48) ) + H ( x (cid:48) , q (cid:48) )2 (cid:19) ×× e i q (cid:48) ( x − x (cid:48) ) u ( x (cid:48) , z ) d x (cid:48) d q (cid:48) (2 π ) . (14)This is a generalization of the operator exponent tech-nique studied in Refs. 16 and 20. One can find that thisscheme results in an exact solution with any finite ∆ ei-ther for a homogeneous dispersive medium, H ( q ) , or foran inhomogeneous medium without spatial dispersion, H ( x ) . For a general case of a dispersive inhomogeneousmedium, H ( x , q , z ) , this scheme converges to the exactsolution for ∆ → .The operator exponent method is rigorously justifiedonly for dissipation-free media described by a Hermitian evolution operator. In this case, there is the effective nu-merical scheme for calculating the operator (14) requiringno more than O ( N ) operations with N being the num-ber of mesh points in 2-dimensional space x . In contrastto the traditional finite difference methods , there areno formal limitations on the step ∆ imposed by grid sizein x . Although the scheme is explicit, the integrationstep is limited only by physical “inhomogeneities” of amedium, namely, by the value of commutator between ˆ H ( x , q ) and ˆ H ( x (cid:48) , q ) : ∆ (cid:28) k (cid:82) | u | d x (cid:90) (cid:12)(cid:12)(cid:12)(cid:12) ∂∂x i ˆ H [ x i u ] − x i ∂∂x i ˆ H [ u ] (cid:12)(cid:12)(cid:12)(cid:12) d x . (15)However, operator (14) may result in a solution withan exponentially growing energy in inhomogeneous dis-sipative media. One of the reasons, already describedin the previous section, is that the dissipation in termsof geometric optics, I m H ( x , q ) > , does not guaran-tee the positive-definiteness of the corresponding oper-ator term ˆ H A . More generally, this is a known issueof finding a proper approximation for physically correctwave absorption in media with a spatial dispersion .In our case, errors in the approximation of the Hermitianpart (e.g. numerical errors due to violation of condition(15)) disturb only the wave beam phase, therefore, theirinfluence is manifested over the large diffraction length.On the contrary, errors in the approximation of the anti-Hermitian part directly affect the beam amplitude, re-sulting in the worst case in an exponential grows of anerror. Our experience tells that this problem is especiallyactual for microwave beams in the electron cyclotron fre-quency range, when the resonant dissipation varies dra-matically both in coordinate and wave vectors spaces.These problems are essentially solved with the modi-fication of anti-Hermitian part of the evolution operatorsuggested with Eq. (12). For numerical solving of the cor-responding quasi-optical equation (13), we propose to usethe “split-step” method, i.e. the sequential solutions forthe separate Hermitian and anti-Hermitian parts. Theintegration step is then u ( x , z + ∆) ≈ exp( − k ∆ ˆ G ˆ G ) (cid:104) ˆ S (cid:48) H [ u ( x , z )] (cid:105) with ˆ S (cid:48) H ≈ exp(i k ˆ H H ) is an approximation to the op-erator exponent of the Hermitian part, it is given byEq. (14) with ˆ H = ˆ H H . Opposite to the Hermitianpart, the dissipation operator does not fit the form suit-able for an efficient computation of the operator ex-ponent (just because it requires two iterations of ˆ G ).Calculation of the operator exponent in this case re-quires about O ( N max λ G k ∆) number of steps, where λ is an eigenvalue of matrix ˆ G , which is typically toomuch for practical applications. Another way is to com-pute exp( − k ∆ ˆ G ˆ G ) approximately by using conventionalmethods such as Adams or Runge–Kutta . Such calcu-lations are also very time demanding, mainly because oflarge eigenvalues that may appear after a double itera-tion of ˆ G and the need to reduce the integration step ∆ to provide a numerical stability.On the other hand, large positive eigenvalues corre-spond to quickly damped or even evanescent modes thattypically do not represent physical interest. This makes itpossible to replace the operator exponent exp( − k ∆ ˆ G ˆ G ) with some approximate operator ˆ D , which correctly de-scribes all weakly damped modes and ensures exponentialdecay for heavily damped and non-propagating modeswith large eigenvalues. To fulfill the first condition werequire that both operators have the same asymptoticbehavior at ∆ → : exp( − k ∆ ˆ G ˆ G ) = ˆ1 − k ∆ ˆ G ˆ G + O (∆ ) , ˆ D = ˆ1 − k ∆ ˆ G ˆ G + O (∆ ) . To ensure that the second requirement, it is sufficient to“cut” the integral kernel of ˆ D at large ∆ . The simplestoperator satisfying the above conditions is ˆ D [ u ] = u − ˆ T [ ˆ T [ u ]] , (16)where operator ˆ T is constructed similar to Eq. (10) butwith kernel T ( x , q ) = φ ( (cid:112) k ∆ G ( x , q )) ,φ ( x ) is an arbitrary one-dimensional function with thefollowing properties: φ ( x ) → x when x → , φ (cid:48)(cid:48) (0) = 0 and φ ( x ) → for x → ∞ . In our code we use the hyper-bolic tangent φ = tanh( x ) . Calculation of the operator(16) reduces to a consistent application of the Fouriertransform to the factorized kernel, so the computationalcomplexity for ˆ D is only by a factor of two greater thanthat for ˆ S (cid:48) H .Now, we are ready to formulate the final explicit nu-merical scheme for quasi-optical equation (13): one inte-gration step is defined as u ( x , z + ∆) = ˆ D [ ˆ S (cid:48) H [ u ]] . (17)In explicit form this reads S ( x , q ) = exp (i k ∆ R e H ( x , q , z ) / ,T ( x , q ) = tanh (cid:112) k ∆ I m H ( x , q , z ) ,u = (cid:90) S ( x , q (cid:48) ) S ( x (cid:48) , q (cid:48) ) e i q (cid:48) ( x − x (cid:48) ) u ( x (cid:48) , z ) d x (cid:48) d q (cid:48) (2 π ) ,u = (cid:90) T ( x , q (cid:48) ) + T ( x (cid:48) , q (cid:48) )2 e i q (cid:48) ( x − x (cid:48) ) u ( x (cid:48) , z ) d x (cid:48) d q (cid:48) (2 π ) ,u = (cid:90) T ( x , q (cid:48) ) + T ( x (cid:48) , q (cid:48) )2 e i q (cid:48) ( x − x (cid:48) ) u ( x (cid:48) , z ) d x (cid:48) d q (cid:48) (2 π ) ,u ( x , z + ∆) = u − u . Verification of this scheme on a set of model exactly solv-able problems is published separately . VI. DESCRIPTION OF HOT MAGNETIZED PLASMA INOPEN MAGNETIC TRAP
In this section, we define the kernel for quasi-opticalevolution operator used in the numerical simulation ofmicrowave plasma heating in open magnetic trap. Ourmodel is based on the dispersion relation for electromag-netic waves propagating in a locally homogeneous hotmagnetized plasma with Maxwellian distribution func-tions. This approach is justified for modern experimentsin which the microwave energy is deposited into relativelydense plasma. Exactly in this case the electrodynamicquasi-optical effects become important, while the kineticeffects related to the rf-driven modification of the elec-tron distribution function are not essential. For example,the Maxwellian distribution formed during ECR heatingof bulk electrons was proved at the GDT device by di-rect laser scattering measurements . Kinetic effects areimportant in the opposite case of low density plasma,such as GAMMA-10 experiment or ECR assisted start-up at the GDT , in which electrodynamics is usuallywell described by standard geometric optics or beam trac-ing, but dielectric response of non-Maxwellian electronsis calculated self-consistently by means of a quasi-linearFokker–Planck code .In a vicinity of the fundamental electron cyclotron har-monic the dispersion relation of the Maxwellian plasmacan be represented as follows : n ⊥ (cid:2) ( ε + − ε (cid:107) )( ε − − n ) + ( ε − − ε (cid:107) )( ε + − n )) (cid:3) == 2 ε (cid:107) ( ε + − n )( ε − − n ) . (18)where n (cid:107) = ck (cid:107) /ω and n ⊥ = cn ⊥ /ω are parallel andperpendicular components of the index of refraction withrespect to the external magnetic field, n = n ⊥ + n (cid:107) , ε − = 1 + XZ ( ζ ) n (cid:107) β e , ε + = 1 − X Y , ε (cid:107) = 1 − X,X = ω p /ω and Y = ω H /ω are standard Stix parametersthat determine the plasma density and the strength of theconfining magnetic field through, respectively, Langmuirand cyclotron frequencies of electrons, β e = (cid:112) T e /mc is the normalized thermal velocity for electrons. All ef-fects associated with the spatial dispersion and resonantdissipation at the fundamental harmonic are defined inthe “warm” part in ε − by so-called plasma dispersionfunction Z ( ζ ) = exp( − ζ ) (cid:32) i √ π sign R e( n (cid:107) ) − (cid:90) ζ exp( t ) dt (cid:33) with argument ζ = (1 − Y ) / ( n (cid:107) β e ) . These formulas areobtained for a weakly relativistic plasma β e (cid:28) , a smallLarmor radius n ⊥ β e (cid:28) Y , and a quasi-longitudinal wavepropagation n (cid:107) (cid:38) Y β e that allows omitting the relativis-tic resonance broadening compared to the non-relativisticDoppler shift. These conditions are definitely met, atleast, in current experiments at the GDT device.Having in mind applications to a long and axially sym-metric devices, such as GDT, in the current implementa-tion of the code we consider the external magnetic fielddirected along the trap axis z . This simplifies the recov-ery of the quasi-optical Hamiltonian from the dispersionequation (18), however this assumption is not essentialfor our technique and may be bypassed if necessary.For strictly longitudinal magnetic field, the quasi-optical variables (4) are mapped to Stix variables as k z = k n (cid:107) , q = k n ⊥ . With this substitution, for eachpoint in the real space and each transverse momentum q we determine a complex solution n (cid:107) ( x , q , z ) of the localdispersion relation (18). Away from the EC resonancethis solution is close to the roots of the biquadratic dis-persion relation corresponding to a cold plasma limit, β e → . This allows us to choose the “right” root, match-ing the studied electromagnetic mode, and use it as astarting point for an iterative search for solution n (cid:107) ofthe transcendent equation (18). Additional control forthe right choice of the mode is provided by requirementof a smoothness of n (cid:107) over x , in case of occasional sharpjumps we use the values at neighboring grid cells to re-cover the correct mode.Following Eq. (8), the complex kernel of the evolutionoperator ˆ H may be defined as H = n (cid:107) ( x , q , z ) − κ ( z ) /k , (19)where κ ( z ) = k R e n (cid:107) ( ˜x , ˜q , z ) is the carrier wave vectorthat corresponds to centers of mass for the wave fielddistribution and its spectrum in the cross-section z , ˜x = (cid:82) x | u ( x , z ) | d x (cid:82) | u ( x , z ) | d x , ˜q = (cid:82) q (cid:12)(cid:12)(cid:82) u ( x , z )e i qx d x (cid:12)(cid:12) d q (2 π ) (cid:82) | u ( x , z ) | d x . Another problem specific for magnetic plasma confine-ment, is evaluation how the absorbed microwave power isdistributed over a trap volume. Here one should keep inmind that this power rapidly redistributes over the mag-netic flux surfaces. In this case, the most relevant char-acteristic of the power deposition is a one-dimensionaldistribution over the magnetic surfaces, usually referredto as a power deposition profile. This distribution maybe found taking into account Eq. (11), as P ( ρ ) = 2 k l ( ρ ) (cid:90) R e( u ∗ ˆ H A [ u ]) δ ( ρ − ρ ( x , z )) d x dz, (20)where the integral is taken over the whole volume occu-pied by a wave beam, equation ρ ( x , z ) = ρ defines themagnetic surface, ρ without arguments is the label of amagnetic surface, l ( ρ ) = (cid:90) δ ( ρ − ρ ( x , z )) d x is the effective perimeter of a magnetic surface at somefixed cross-section z = z , for an axially symmetric mag-netic configuration l = 2 πρ . The code QOOT (Quasi-Optics for Open Traps) is de-veloped to solve the quasi-optical equation (13) with iter-ative procedure (17) for the evolutionary operator definedby kernel (19). The code calculates the power depositionprofile (20) in the geometry of an open magnetic trap andvisualizes the results. We use the Cartesian coordinatesfor the wave amplitudes and output of the results, andthe internal cylindrical coordinates for the Hamiltoniantaking advantage of the axial symmetry of the plasmaconfiguration. As mentioned in the introduction, the de-velopment is based on our previous quasi-optical code
LAQO , created for toroidal magnetic traps . However,since conditions for the wave propagation in open mir-ror traps and toroidal traps are significantly different, wehave eventually developed an entirely new code. VII. SIMULATION OF ECR PLASMA HEATING AT THEGDT DEVICE
The basic ECR heating scheme at the GDT relies onradiation trapping by a non-uniform plasma column .This effect is caused by a dependence of the plasma re-fraction of the magnetic field strength. The radiation islaunched through a side of the plasma column at a highmagnetic field (close to the magnetic mirror). As a mi-crowave beam propagates in plasma in the direction ofthe trap center, the magnitude of the magnetic field de-creases resulting in conditions for an internal reflectionfrom the plasma-vacuum boundary. The plasma columnacts as kind of waveguide, heterogeneous in both trans-verse and longitudinal directions, whereby the radiationis delivered to the ECR surface. The geometric opticsmay not be valid in the vicinity of the wave reflectionsurface due to caustic formation, and near ECR surfacedue to sharp inhomogeneous damping of the wave field.This casts doubt on the applicability of geometrical op-tics approximation.Therefore, the first task for the quasi-optical code isto check our previous results obtained with ray-tracing.In particular, we repeated simulations aimed at optimiza-tion of the ECR heating efficiency at the GDT conditions.Previously, it has been found that the efficient ECR heat-ing can be achieved in two distinguished regimes char-acterized by very different distributions of the absorbedmicrowave power . Switching between these regimes,as well as a fine-tuning of the microwave power deposi-tion profile can be controlled, at specific conditions of theGDT experiment, by a relatively small adjustment of theexternal magnetic field in the close proximity of the ECresonance.The results are shown in Figures 1–5. For two-dimensional visualization of a quasi-optical beam in ( x, y ) plane (face view along the trap axis) and ( y, z ) plane (sideview) we use the wave intensity integrated along the z and x axis, correspondingly: J z ( x, y ) = (cid:90) χ | u | dz, J x ( y, z ) = (cid:90) χ | u | dx. Figure 1. (color online) Simulations of the “narrow power deposition” regime before the ECR heating. Distributions of thewave beam intensity (a – side view, b – face view) and the power deposition profile related to the trap center (c). The solidlines in panels (a) and (b) show the geometric optics rays, dashed and dotted curves indicate the EC resonance and the plasmaboundary, respectively. The efficiency of absorption is almost 100 %.Figure 2. (color online) Simulations of the “narrow power deposition” regime after the ECR heating. The efficiency of absorptionis almost 100 %.
Here | u | characterizes the density of the energy fluxalong the z axis, and χ = | v gr | / ( v gr z ) projects thatflux to the direction of the the group velocity v gr . Inexplicit form χ = 1 + (cid:18) ∂ R e k z ∂k x (cid:19) q = ˜q + (cid:18) ∂ R e k z ∂k y (cid:19) q = ˜q . To improve the contrast, the logarithmic color scale isused corresponding to the levels of ln J , which allowedus to visualize the wave field in caustics. The powerdeposition profiles P ( ρ ) are calculated using Eq. (20),applied at the central cross-section where the confiningmagnetic field has its minimum.Figures 1 and 2 show the results of simulation, the wavebeams and the power deposition profiles, in the “narrow Figure 3. (color online) Simulations of the “broad power deposition” regime before the ECR heating. The efficiency of absorptionis about 80 %.Figure 4. (color online) Simulations of the “broad power deposition” regime after the ECR heating. The efficiency of absorptionis about 80 %. power deposition” regime. In this regime the record val-ues of the electron temperature (up to 1 keV) have beenachieved in the GDT experiments. First and second fig-ures correspond, respectively, to states before and afterthe ECR heating modeled for the experimentally mea-sured plasma density and electron temperature profiles.For comparison, we show the ray-tracing results. In thiscase, the power deposition profiles obtained by two dif-ferent methods coincide quite well. Some discrepancy isobserved on the periphery of the plasma column. This corresponds to the absorption of radiation after reflec-tion in the vicinity of the caustic surface where geomet-ric optics is violated (caustics are indicated as crossing ofneighboring ray trajectories). Note, that the geometricoptics rays reproduce the quasi-optical beam satisfacto-rily both before and inside the caustic region.Figures 3 and 4 show the same plots for the “broadpower deposition profile” regime. In the experiments inthis mode, we observe a pronounced increase of a to-tal plasma energy related to an improved confinement of0
Figure 5. (color online) Simulations of the “improved broad power deposition” regime before the ECR heating. The efficiencyof absorption is almost 100 %. the hot ions. There is a better agreement between thequasi-optical and geometro-optical power deposition pro-files since the ratio of microwave power deposited in thevicinity of the caustic is much lower. As well as the ray-tracing, the quasi-optical simulations predict incompleteabsorption of the microwave power in this regime.Figure 5 shows predictions for the new regime with“improved broad power deposition”. This mode has notbeen yet demonstrated experimentally; its implementa-tion will be eventually possible after up-grade of the mir-ror magnetic coil of the GDT system which is currentlyon progress. All plots correspond to the stage beforethe microwave heating because the experimental data onplasma profiles after the heating is not available. In thismode, an extended region of the caustic is formed, andthe ratio of the power deposition after the caustic is muchhigher than in the regimes discussed previously. There-fore, the power deposition profiles predicted by the quasi-optics and the ray-tracing vary considerably. We areaware that the quasi-optical modeling provides a moreadequate description, but confirmation of this statementrequires further studies and experiments, the results ofwhich will be published separately.Finally, it can be concluded that geometric opticsseems to be a reasonable agreement with the more accu-rate quasi-optical simulations for the most heating sce-narios that are already realized in the GDT experiment.However, this conclusion does not exclude a possible im-pact of spatial dispersion on the resonance dissipationand of diffraction losses near the caustic surfaces in fu-ture experiments with more optimized heating scenarios.
VIII. SUMMARY
The quasi-optical model of propagation of wave beamsin high-temperature magnetically confined plasmas, de-veloped earlier for toroidal traps, is generalized over openmagnetic systems. The specifics of microwave heatingin modern open traps require substantial improvementof the early quasi-optical theory, associated with moreaccurate description of the effects of spatial dispersionin the region of resonant wave dissipation. As a result,a new form for the quasi-optical equation is proposed,see Eq. (13). Basing on this equation the universal code
QOOT is developed for simulation of electron cyclotronplasma heating, which allows to resolve the diffraction,dispersion and aberration effects in the propagation andabsorption of electromagnetic wave beams in open traps.The code is used to verify the results of optimizing theefficiency of the ECR heating in the large mirror trapGDT, previously obtained by using ray-tracing withingeometric optics approximation. First quasi-optical sim-ulations justify the possibility to control the radial distri-bution of the deposited microwave power by local mod-ification of the magnetic field in near the EC resonanceand, in particular, the ability of effective on-axis heatingthe electron component at the GDT conditions. The in-fluence of wave caustics on a localization of the depositedmicrowave power is demonstrated.
ACKNOWLEDGMENTS
This work was supported by the Russian Science Foun-dation (grant No 14-12-01007). The authors would like1to thank Dr. Alexander Solomakhin from the Budker In-stitute for his support with magnetic configurations andray-tracing modeling for the GDT.
REFERENCES A.V. Vodopyanov, S.V. Golubev, V.G. Zorin, A.Yu. Kryachko,A.Ya. Lopatin, V.I. Luchin, S.V. Razin, and A.N. Smirnov, Tech.Phys. Lett. (12), 1075 (2000). T. Cho, V. P. Pastukhov, W. Horton, T. Numakura, M. Hirata,J. Kohagura, N. V. Chudin, and J. Pratt, Phys. Plasmas ,056120 (2008). T. Tamano, Phys. Plasmas , 2321 (1995). T. Saito, K. Ishii, A. Itakura, M. Ichimura, Md. K. Islam,I. Katanuma, J. Kohagura, Y. Tatematsu, Y. Nakashima, T.Numakura, H. Higaki, M. Hirata, H. Hojo, M. Yoshikawa, K.Sakamoto, T. Imai, T. Cho, and S. Miyoshi, J. Plasma FusionRes. (4), 288 (2005). T. C. Simonen and R. Horton, Nucl. Fusion , 1373 (1989). P.A. Bagryansky, Yu.V. Kovalenko, V.Ya. Savkin, A.L. Solo-makhin, and D.V. Yakovlev, Nuclear Fusion , 082001 (2014). P.A. Bagryansky, A.G. Shalashov, E.D. Gospodchikov, A.A.Lizunov, V.V. Maximov, V.V. Prikhodko, E.I. Soldatkina, A.L.Solomakhin, and D.V. Yakovlev, Physical Review Letters ,205001 (2015). P.A. Bagryansky, E.D. Gospodchikov, Yu.V. Kovalenko, A.A.Lizunov, V.V. Maximov, S.V. Murakhtin, E.I. Pinzhenin, V.V.Prikhodko, V.Ya. Savkin, A.G. Shalashov, E.I. Soldatkina, A.L.Solomakhin, and D.V. Yakovlev, Fusion Science and Technology , 87 (2015). P.A. Bagryansky, A.V. Anikeev, G.G. Denisov, E.D. Gospod-chikov, A.A. Ivanov, A.A. Lizunov, Yu.V. Kovalenko, V.I. Ma-lygin, V.V. Maximov, O.A. Korobeinikova, S.V. Murakhtin, E.I.Pinzhenin, V.V. Prikhodko, V.Ya. Savkin, A.G. Shalashov, O.B.Smolyakova, E.I. Soldatkina, A.L. Solomakhin, D.V. Yakovlev,and K.V. Zaytsev, Nuclear Fusion , 053009 (2015). T.C. Simonen, Journal of Fusion Energy , 63 (2016). A. Shalashov, E. Gospodchikov, O. Smolyakova, P. Bagryansky,V. Malygin, and M. Thumm, Physics of Plasmas , 052503(2012). А. G. Shalashov, E. D. Gospodchikov, O. B. Smolyakova, P.A.Bagryansky, V. I. Malygin, and M. Thumm, Problems of AtomicScience and Technology, Series: Plasma Physics , 49 (2012). P.A. Bagryansky, S.P. Demin, E.D. Gospodchikov, Yu.V. Ko-valenko, V. I. Malygin, S.V. Murakhtin, V.Ya. Savkin, A.G. Sha-lashov, O.B. Smolyakov, A.L. Solomakhin, M. Thumm, and D.V.Yakovlev. Fusion Science and Technology , no 1T, 40 (2013). A.A. Balakin, Radiophysics and Quantum Electronics , 502(2012). A.A. Balakin, M.A. Balakina, A.I. Smirnov, and G.V. Permitin,Plasma Phys. Reports , 302 (2007). A.A. Balakin, M.A. Balakina, and E. Westerhof, Nuclear Fusion , 065003 (2008). N. Bertelli , A.A. Balakin, E. Westerhof, and M.N. Buyanova,Nucl. Fusion , 115008 (2010). A.A. Balakin, Radiophysics and Quantum Electronics , 472(2012). A.A. Balakin and E.D. Gospodchikov, Journal of Physics B:Atomic, Molecular and Optical Physics , 215701 (2015). G.M. Fraiman, E.M. Sher, A.D. Yunakovsky, and W.Laedke,Physica D. , 325 (1995). W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery,
Numerical Recipes: The Art of Scientific Computing (Cam-bridge University Press, 2007). A.A. Balakin, M.A. Balakina, A.I. Smirnov, and G.V. Permitin,Plasma Phys. Reports , 533 (2008). M. Brambilla,
Kinetic Theory of Plasma Waves: HomogeneousPlasmas (Clarendon Press, Oxford, 1998). V.L.Ginzburg,
Propagation of electromagnetic waves in plasma (Pergamon Press, 1970). M.D. Tokman, E. Westerhof, and M.A. Gavrilova, Plasma Phys.and Contr. Fusion , 91 (2000). M.D. Tokman, E. Westerhof, and M.A. Gavrilova, Journal ofExperimental and Theoretical Physics (6) 1141 (2000). M. D. Tokman, E. Westerhof, and M.A. Gavrilova, Nuclear Fu-sion, , 1295 (2003). A.G. Shalashov, A.A. Balakin, E.D. Gospodchikov, T.A. Khu-sainov, and A.L. Solomakhin,
Quasi-optical description mi-crowave plasma heating , accepted by JETP, 2016. D. V. Yakovlev, A. G. Shalashov, P. A. Bagryansky, E.D. Gospodchikov, V. Ya. Savkin, and A. L. Solomakhin.Electron-cyclotron plasma startup in the GDT experiment.arXiv:1607.01051 [physics.plasm-ph] ( accepted by Nucl. Fusion,2016 ). B.W. Stallard, Y. Matsuda and W.M. Nevins, Nucl. Fusion (2), 213-223 (1983). E. D. Gospodchikov and E. V. Suvorov, Radiophysics and Quan-tum Electronics48