Analysis of equations arising in gyrotron theory∗

Abstract. The gyrotron is a microwave source whose operation is based on the stimulated cyclotron radiation of electrons oscillating in a static magnetic field. Powerful gyrotrons can be used to heat nuclear fusion plasma. In addition, they have found a wide utility in plasma diagnostics, plasma chemistry, radars, extra-high-resolution spectroscopy, high-temperature processing of materials, medicine, etc. However, the main application of gyrotrons is in electron cyclotron resonance heating in tokamaks and stellarators. Equations describing gyrotron operation are ordinary differential equations and Schrödinger type partial differential equations. The present paper provides a survey of the analytical and numerical results concerning these equations obtained by our group in the last decade.


Introduction
First gyrotrons were built about three decades ago [1]. Various aspects of these microwave tubes are described in the monograph [2]. The radio frequency (RF) radiation is excited by gyrating electrons bunched near the phase in which they yield their energy to the highfrequency field. The phase bunching is due to the relativistic dependence of the electron mass on its velocity. This makes rotation of the decelerated electrons faster and that of the accelerated ones slower. The typical frequency range of gyrotrons working at the fundamental harmonic is between 20-200 GHz.
Operation of a gyrotron can be described by the system of equations derived by physicists [3]. This system consists of equations of two types -ordinary differential equation for the electron motion and the Schrödinger type partial differential equation for RF field. The simplest illustrative case deals with a single oscillating mode, but in realistic cases several RF modes have to be taken into account.
The first attempt to give a general mathematical analysis of the basic gyrotron equation which describes the electron interaction with the high-frequency field in a gyrotron resonator and to classify possible electron trajectories in the phase space was undertaken in the paper [4]. In that paper it was proved that in the case when the RF field is represented by a Gaussian type function, the solutions of the gyrotron equation are asymptotically equal to the solutions of the corresponding unforced equations. In the paper [5] this analysis was extended to the case of an infinitely long resonator in which the RF field is represented by a propagating wave. A completely realistic case -a finite length resonator and a true longitudinal distribution of RF field, was mathematically analyzed in the article [6]. In the paper [7] the complicated alternating sequences of regions of stationary, automodulation and chaotic oscillations in the plane of the generalized gyrotron variables -cyclotron resonance mismatch and dimensionless current -were determined (Fig. 4). Results of various subsequent studies have been published in a number of papers [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23].
We believe that the results of our long lasting work on studying various mathematical aspects of equations arising in gyrotron theory are useful for gyrotron theoreticians in their efforts to further develop these important microwave generators.

Basic equations describing single mode gyrotron oscillations
The gyrotron resonator is formed by a section of a cylindrical waveguide and tapered input and output sections. Taking into account the weak temporal dependence of the RF wave amplitude, the wave equation can be transformed into the Schrödinger type partial differential equation where ζ ∈ (0, ζ out ), τ ≥ 0 are the normalized axial and temporal coordinates, ζ out is the length of the interaction space, the complex function f (ζ, τ ) is the normalized amplitude of RF field, p(ζ, τ, θ 0 ) is the dimensionless complex transverse orbital momentum of the electron, p = 1 2π 2π 0 p(ζ, τ, θ 0 ) dθ 0 is the averaged value of p(ζ, τ, θ 0 ) with respect to the parameter θ 0 in the initial condition, I 0 is the dimensionless beam current, δ is the frequency mismatch, and i is the imaginary unit i = √ −1. The function p(ζ, τ, θ 0 ) is determined by the equation of electron motion where the real constant ∆ is the cyclotron resonance mismatch. The initial condition for (2) is where θ 0 ∈ [0, 2π) is the initial angle. Equation (2) can be transformed also in the form wherep = exp(i∆ζp). In this case the right hand side of equation (1) can be written in the following form: Equation (1) should also be supplemented with initial and boundary conditions: where f 0 (ζ) is a given function. At the output cross section ζ = ζ out two kinds of boundary conditions have been proposed: 1. The condition for the radiation of the Fourier components of the RF field is of Here γ is a positive parameter. If γ = ∞, then we have boundary conditions of first kind.

Reflectionless (nonlocal) boundary condition in the integral form
The electron perpendicular efficiency η ⊥ , which describes the extraction of the electron orbital momentum from the beam, is given by the expression For stationary oscillations the field amplitude f can be represented as f =f exp(iγ 2 t) with boundary condition (5) [3].

Basic equations describing multimode gyrotron oscillations
Multimode gyrotron oscillations are described by the nonlinear system of complex partial differential equations: where s = 1, . . . , m, f s (ζ, τ ) and p(ζ, τ, φ, θ 0 ). Here δ s = δ s (ζ) describes variation of critical frequencies, ψ s = a s τ + m s φ is the phase of the s-mode, a s and m s are physical constants, φ is the azimuthal coordinate, and 0 ≤ φ < 2π. The equation (9) has to be supplemented with the initial condition where 0 ≤ θ 0 < 2π is the initial angle, and equations (8) with the initial and boundary conditions where f s0 (ζ) are given functions and Here γ s are positive parameters. Parameters of equations (δ s , ∆ s , γ s , ζ out ) depend on different physical quantities [2].

Gyrotron equation describing electron motion
In [4] we presented a detailed mathematical analysis of the simplest gyrotron equation in the so-called cold-cavity approximation when the high-frequency field in a gyrotron resonator depends only on the geometry of the resonator but not on the electron motion, i.e., g does not depend on p. In this case the high-frequency field in a cavity is usually represented by a Gaussian type function where F is the beam-high-frequency field coupling and µ is the dimensionless length of the resonator. We assume that g is a general scalar function decaying fast enough, such that the integral ∞ ζ|g(ζ)| dζ converges. Note that the Gaussian g satisfies this condition. We prove that (13) and the corresponding unforced equation are asymptotically equal, which means that each solution of (13) corresponds to a solution of (14) with changed initial condition such that difference of these two solutions approaches zero when ζ → +∞. As a result, we conclude that there are four different patterns of behaviour of the solution of (13) in the complex plane C.
In [5] we consider gyrotron equation (13) in the case of an infinitely long idealized resonator in which RF field is represented by a propagating wave. In this case g can be approximated as a wave moving in the ζ → +∞ direction Introducing the new function q = p exp(ikζ) in (13), we obtain the equation Another important property of (15) should be mentioned. It is an autonomous equation, i.e., its right-hand side does not depend on ζ. This makes it possible to obtain illustrative trajectories in the phase space. By introducing the notation q = x + iy in (15), we obtain the following dynamical system: The trajectories of (16) are fourth order algebraic curves: singular points, cycles and homoclinic trajectories. In principle we could integrate the system of equations (16) using the Jacobi elliptic functions and in this way obtain the parametric equations for the electron motion. Introducing the notation D = 27F 2 − 4(1 − ∆ − k) 3 , we can distinguish three cases.
1. If D > 0, then the equation (16) has one stationary point -center. Other trajectories are cycles around this center.
2. If D = 0, then the phase portrait is similar to the previous case. One new stationary point is born on one cycle and instead of the corresponding cycle homoclinic trajectory appears.
3. If D < 0, then there are three stationary points, one of which is saddle and two others are centers. Corresponding homoclinic trajectories are Pascal limacons.
Taking into account what has been said above, we can state that any solution of (13) can be represented in the form which is a product of two functions with different periods (except the function corresponding to a homoclinic trajectory).
Similar analysis is performed in [6] employing physically transparent Hamiltonian method using the angle-action variables in (15). The dimensionless coordinate ζ can be regarded as time and we obtain the time-independent integrable Hamiltonian with one degree of freedom, where E is the total energy. This Hamiltonian satisfies the following equation A detailed numerical study is performed of the behaviour of electron trajectories in a completely realistic case when a resonator has a finite length and the perturbation (RF field) is an aperiodic function for small ζ and an approximately periodic function for large ζ.
For computational purposes we have interpolated the function g by cubic splines and integrated the corresponding dynamical system by means of a Runge-Kutta method with high accuracy. We have found that the motion of the electrons in the vicinity of some particular initial angles is very sensitive with respect to the exact value of the initial angle. This sensitivity allows us to conclude that the trajectories of the electrons for these particular initial conditions may become chaotic.

Numerical experiments
For solving the Schrödinger type differential equations (1) and (8) we normally use the implicit finite difference scheme with the uniform spatial grid with the step ∆h and the temporal grid with the step ∆τ . We obtained a linear algebraic three diagonal system of equations from which we can find f for a fixed value of p. The differential equations (2) and (9) for fixed value of f were solved with the Runge-Kutta method using the routines "mgear" from the package MAPLE. Calculations show that the results strongly depend on the spatial and temporal steps. We consider the multimode system with seven modes. The character of solutions of system (11) can be seen in Fig. 3. The voltage step is ∆U = 0.1 (kV) and the temporal step is ∆τ = 0.1 (ns). The values of |f (ζ out )| show that in the time interval (0, 800) for U ≤ 75 (kV) the gyrotron oscillates in the second mode, in the interval (900, 1300) for U ≤ 85 (kV) in the fourth mode, and in the interval (1300, 1700) for U ≤ 90 (kV) in the first mode.

Single mode gyrotron equations with reflectionless boundary conditions
The nonstationary oscillations in gyrotrons were investigated on the basis of the system of equations (1) and (2) with conditions (3), (4), and reflectionless boundary condition in the integral form (6) cf. [7]. By noting that the singularity in the boundary condition (6) can be isolated if the function ∂f (ζ, τ )/∂ζ satisfies the Lipschitz condition in the vicinity of τ , we can use a difference approximation of (6). In the calculations the implicit finite difference scheme was used with uniform spatial grid for δ = 0, If the problem (1)-(4), (6) admits a stationary solution, the electron efficiency η ⊥ can be computed as in [16]. There we use a spectral method (Chebyshev differentiation matrices) to approximate the spatial derivatives and terminate the iterations over time when a stopping criterion is met. η ⊥ was found to depend sensitivetely on the frequency mismatch δ = δ(ζ). The maximum value of the electron efficiency (7) η ⊥ = 0.7486 was found for ∆ = 0.6, δ = 0, and I = 0.01 which is in a good agreement with the results obtained in [7].
The stability of various finite difference schemes for the same problem on uniform spatial grids was studied in [11,24]. It was shown that certain three layer schemes (Richardson's and Simpson's) are absolutely unstable for the Robin boundary condition; conditional stability was established in the limiting case γ = ∞. The weighted two-layer finite difference scheme was shown to be stable for value of weight coefficient σ ≥ 1/2 [11].
In [10] the initial-boundary value problem (1), (4), and (5) (I 0 = 0) was solved by using the Fourier, finite-differences and degenerate matrix methods. The degenerate matrix method is based on using differentiation matrix for derivatives with respect to ζ with the Chebyshev and Legendre nodes. Finite-difference and degenerate matrix methods give the same order of accuracy, but for this degenerate matrix method needs approximately three times smaller number of grid points than finite-difference (especially for δ(ζ) = const). For δ = δ(ζ) = tanh(7ζ − 3.5ζ out ) degenerate matrix method is more effective.
Assuming |p| 2 = 1 + C (C ≤ 0), I 0 = 0 the system (1), (2) can be reduced to an integro-differential equation that was considered in [13,14]. The temporal evolution of the amplitude of the quasi-stationary solutions f (τ, ζ) = h(ζ) exp(iατ ) depends on the imaginary part of α, which can be calculated either from a transcendental equation arising from an eigenvalue problem for the continuous version or by discretization of the continuous problem, the comparison allows to better understand the algorithm. 80 eigenvalues for the discrete case are shown in Fig. 7; the tail consisting of 20 rightmost eigenvalues are of parasitic character. The method of lines with nonuniform spatial grid (Chebyshev's points) is shown to be a suitable numerical approach.  In [15] the problem (1)-(4) was solved with the local condition (5). Here η ⊥ was obtained using (7), also f | ζ=ζout was calculated. Using the method of stationarity for the system of equations (1) and (2), the term ∂f /∂τ was approximated with finite differences of first order. The limiting solutions, when τ → ∞, were obtained with finite difference method on nonuniform grid with grid points as the roots of Chebyshev polynomials. Using these grid points, the derivatives in equations (1) and (2) were approximated with corresponding matrixes by means of the Lagrange interpolation. These results were obtained for the following values: ζ out = 15, γ ∈ [0, 4], and I 0 ∈ [0.01, 0.1].

Natural upper limit of power for gyrotrons
In [8,12] nonstationary oscillations in gyrotrons were investigated both in the longitudinal and azimuthal coordinates. Similarly as in [18], the electron motion and the temporally and spatially dependent high-frequency field in the gyrotron resonator were calculated by means of the following system of partial differential equations Here ζ as before is the dimensionless longitudinal coordinate, but ξ is the dimensionless azimuthal coordinate. Using the new variable w = (4β 2 /β 2 ⊥ )(τ − ξ), where β and β ⊥ are parallel and perpendicular components of the electron velocity normalized to the speed of light, from the system of equations (18) and (19) we obtain the system Equation (21) can be solved by means of the method of characteristics taking into account the periodicity condition where W is the dimensionless period in the azimuthal direction, w ∈ (0, W ). If the function p = p(ζ, w, τ, θ 0 ) is given, the function f can be found from the equation (20) by means of finite difference scheme [7]. The initial condition is in the form where f 0 (ζ) = sin(πζ/ζ out ) is RF field profile obtained in the cold cavity approximation.
The calculations resulted in specific critical values of m (azimuthal index of the mode), beyond which gyrotron oscillations become irregular in the azimuthal direction for particular values of ∆ and I 0 (calculations were performed for five specific coaxial gyrotrons). It was argued that there exists a natural upper limit on the power (about 4 MW) that can be generated by a single gyrotron. It should be noted that coaxial gyrotrons have proven excellent potential to become useful high-power microwave sources. The theory of coaxial gyrotrons is somewhat more complex (see review article [25] but has significant synergy with conventional cavity gyrotron analysis).
In Fig. 8 we can see the temporal dependence of the efficiency. Summarizing the results of numerical experiments, we can list the following items: • Depending on the beam current I 0 , three kinds of the asymptotic temporal behaviour of the modulus of the RF field amplitude are observed: stationary, periodic and non-periodic. • For the single mode case the uniform grid schemes (Richardson's, Simpson's etc.) turn out to be absolutely unstable. • Complexity of differential operator eigenvalues yield a large number of terms in Fourier series. • In the case of a uniform grid, imaginary parts of eigenvalues asymptotically tend to zero, causing additional errors in finite difference schemes. • In the case of a nonuniform grid (grid points are roots of Chebyshev polynomial) the strongly oscillating parasitic eigenvalues must be rejected. www.mii.lt/NA