On the numerical scheme employed in gyrotron interaction simulations

We report on the influence of the numerical scheme employed in gyrotron interaction simulations. Results obtained with the Crank-Nicolson scheme are compared with those obtained with the Backward Time – Centred Space (BTCS) fully implicit scheme. We present realistic cases where, for discretisation parameters in the range usually used in gyrotron simulations, the results can be very different. Hence, the numerical scheme used can be responsible for obscuring the underlying physics if its convergence is not tested carefully.


Introduction
Gyrotrons are the leading microwave sources for high-frequency waves needed for plasma heating and current drive in magnetic confinement fusion experiments.The European Gyrotron Consortium (EGYC) is developing a 170 GHz gyrotron in support of Europe's contribution to ITER [1][2][3].In a gyrotron, a helical electron beam guided by a magnetostatic field delivers energy to a high-frequency (RF) wave through electron cyclotron resonance.Simulations of the beam-wave interaction are important for guiding gyrotron design and in support of the experiments.In order for the interaction codes to be used as designing and/or supportive tools, fast simulations are necessary.
There are four European interaction codes, available at EGYC, based on the slow-time scale approximation and thus capable of fast simulations: EURIDICE [4], SELFT [5], COAXIAL [6], and TWANG [7].The results of the codes are in agreement in the majority of cases that have been tested.Some discrepancies that appear are understood and attributed to the difference in the interaction models.However, in some cases of interest, the results of COAXIAL have been significantly different from those of the other codes.This triggered an extensive investigation which indicated that this discrepancy was a consequence of the difference in the employed numerical scheme: COAXIAL uses a scheme which is first-order in time, whereas the rest of the codes use schemes that are secondorder in time.To clarify the situation, EURIDICE was modified in order to provide the option to use either of the two schemes.The results confirmed that the previously observed discrepancies were due to the difference in the numerical scheme.In the next sections we present the beam-wave interaction equations and the numerical schemes for their integration, together with realistic cases where, for discretisation parameters in the range usually used for efficient simulations, the two schemes give results that are different not only quantitatively, but also qualitatively.

Interaction model and numerical schemes
Many variants of the slow-time-scale, self-consistent model for the beam-wave interaction in the gyrotron cavity exist in the literature ( [8] and references therein).The model consists of the equations for the electron motion and the equations for the field profiles of the TE modes of the cavity, which are interacting with the electrons.For the purpose of the present discussion, we use the model in the code COAXIAL [9]: ( ) The physical process described by ( 1)-( 2) is the interaction, at the fundamental cyclotron frequency, of a weakly relativistic electron beam having no spreads in electron energy, velocity and guiding centre, with a set of TE modes that are close to cut-off and in resonance with the beam.It is assumed that the cavity walls are perfectly conducting and that the axial electron momentum is constant.The latter approximation is valid for weak tapering of the magnetostatic field.Detailed definition of the normalised quantities in ( 1)-( 2) can be found in [9]: p(ζ) is the electron complex transverse momentum normalised to its initial absolute value, f s (ζ, τ) is a complex function describing the envelope of the RF field axial profile of the s-th mode, ζ (∝ z) and τ (∝ t) denote axial distance along the cavity, and time, respectively.The quantities p and f s are average quantities over the cyclotron period (gyro-averaged, slow-time-scale equations).The parameter ∆ B (ζ) describes the magnetic field tapering along the cavity axis and the parameter δ s (ζ) accounts for the mild axial inhomogeneity of the cavity wall.The normalised current I s (ζ, τ) incorporates the coefficient describing the coupling between the s-th mode and the electron beam, along with the stored energy coefficient.The parameter ∆ s is the frequency mismatch of the s-th mode and ψ s (τ, φ) is the mode phase, where φ is the azimuthal position of the electron guiding centre.The initial condition for (1) is p(0) = exp(iθ), 0 ≤ θ <2π, and the boundary conditions for (2) are reflectionless boundary conditions at the carrier frequency of the s-th mode, i. e. at the frequency of the fast oscillations of the mode.
The models used in EURIDICE, SELFT and TWANG are more complicated than ( 1)-( 2), accounting in a more elaborate way for several physical effects by dropping some of the assumptions of ( 1)-( 2).However, regarding the scope of this paper, these differences do not play an important role.As a result, what we state by using ( 1)-( 2) as a reference model applies also to the other models.
The equation of motion ( 1) is solved in the codes using some Runge-Kutta or equivalent scheme.Here we will concentrate on the numerical solution of the equation (2) for the field profile, which is a two-dimensional partial differential equation of the form: For the discretisation of (3) we assume a grid for the z-axis with points z k = k∆z, k = 0, 1, 2… N z -1 and a grid for the t-axis with points t n = n∆t, n = 0, 1, 2…, where ∆t and ∆z are the time-and space-04017-p.2 discretisation step, respectively.Using the notation ( , ) A z t A ≡ , we can come up with the following two finite difference schemes [10] for the numerical integration of (3): Backward Time -Centred Space scheme (BTCS): Crank-Nicolson scheme (C-N): The BTCS scheme, also known as the fully implicit scheme, is second-order in space and first-order in time, whereas the C-N scheme is second-order in both space and time.The code COAXIAL uses the BTCS for the integration of the field equation, whereas EURIDICE and SELFT use the C-N.TWANG uses a finite element method, which is also second order in time and space, like the C-N.
3 Numerical results

Dynamic after-cavity interaction
Recently, in some cases where the simulated interaction region is not confined to the gyrotron cavity but it also incorporates the non-linear uptaper after the cavity, dynamic after-cavity interaction (ACI) has been observed [11].In dynamic ACI, a mode is interacting with the beam in the cavity at some frequency and, at the same time, the mode is also interacting with the beam at the non-linear uptaper at a different frequency.The existence of two frequencies results in a strong beating which is apparent both in the field profile and in the output power.An example of dynamic ACI in the 1 MW, 140 GHz gyrotron for W7-X [12], as simulated by EURIDICE, is shown in Figure 1.Along with the operating frequency at 140 GHz, parasitic frequencies at ~131 GHz are excited due to interaction at the uptaper, deforming the RF field profile at this area (z > 80 mm).Dynamic ACI is predicted, with good agreement, by EURIDICE, SELFT and TWANG in several cases.However, COAXIAL did not predict such interaction.An extensive investigation excluded the differences in the employed interaction models from being the reason for the discrepancy, and indicated the difference in the numerical schemes used.This became apparent from a convergence Fig. 1.Simulated dynamic after-cavity interaction in the W7-X gyrotron.Left: Resonator wall, magnetic field, and absolute field profile of the operating TE 28,8 mode versus the resonator axis.Right: Spectrum at the end of the interaction region.
study performed with EURIDICE, after its modification to be able to use both C-N and BTCS.Results are presented in Figure 2. The test case is the design of a 1 MW, 170 GHz gyrotron, operating at the TE 32,9 mode, which is the EU fall-back solution for ITER [1].The results reveal that, compared to the C-N, the BTCS needs a three orders of magnitude smaller time-discretisation step ∆t to show the dynamic ACI.This results in prohibitively long computation time.Following this result, COAXIAL was finally able to simulate dynamic ACI by using an adequately small ∆t.
In order to illustrate more the dissipative action of the BTCS which suppresses dynamic ACI, we performed a similar convergence study by enforcing a field profile with ACI as an initial condition Test case as in Figure 1 with a field profile exhibiting ACI enforced as initial condition.
In the graphs, the output power of the operating TE 32,9 mode versus time is shown.
for the field profile f(z, t) and let it evolve over time.The results (in terms of output power) shown in Figure 3, confirm that BTCS needs a much smaller time-discretisation step ∆t to converge to the results of C-N.The power fluctuations have a period of ~0.1 ns, which is the beating period between the operating (170 GHz) and ACI (~160 GHz) frequencies.At this point we should stress that dynamic ACI is an effect that is currently under extensive analysis regarding both the experimental findings and the capability of the present interaction models and codes to describe it correctly [11][12].Discussion on the validity of the models and codes in treating dynamic ACI is beyond the scope of this paper, which focuses on the influence of the numerical scheme on the calculations.

Frequency-tunable gyrotron
In contrast to the dynamic ACI case, the case addressed in this section lies well into the region of validity of the interaction models and codes.It refers to a frequency-tunable gyrotron for Nuclear Magnetic Resonance spectroscopy, enhanced by Dynamic Nuclear Polarisation (DNP-NMR) [13].
The gyrotron operates at the TE 7,2 mode at the frequency of about 263 GHz and frequency tuning can be realised by changing the magnetic field.The frequency tunability is achieved by the successive excitation of higher axial harmonics of the operating mode.The output power remains at a level above 10 W.
In Figure 4 we show the output power of this gyrotron as a function of time, where the magnetic field is increasing linearly with time from 9.65 T up to 9.94 T. The operating frequency changes correspondingly from 262.93 GHz up to 269.82 GHz.The power fluctuations correspond to the successive excitation of higher axial harmonics of the TE 7,2 mode.In this example ohmic losses are not taken into account.The results are obtained by EURIDICE using both the BTCS scheme and the C-N.Again the BTCS results converge to the C-N results for adequately small time-discretisation step ∆t.For larger time-discretisation step BTCS predicts much narrower frequency tunability.The essential conclusion is that the results of C-N converge to a solution for relatively large timediscretisation step ∆t, while the BTCS results obviously depend on the chosen ∆t and converge to the result of C-N only for very small ∆t.

Discussion
We have compared results of simulations of the beam-wave interaction in gyrotrons for two realistic cases.We use two different numerical schemes -the BTCS fully implicit scheme and the Crank-Nicolson scheme.It is shown that the BTCS needs about three orders of magnitude smaller timediscretisation step ∆t in order to converge to the results obtained by the C-N.This makes the BTCS In the graphs, the output power of the operating TE 7,2 mode vs. time is shown.The magnetic field is increasing linearly with time 04017-p.5

EC-17 Workshop
EPJ Web of Conferences very inefficient in the simulated cases, as the required simulation time becomes very long.We should expect such behaviour, since the BTCS is first-order in time whereas the C-N is second-order in time.However, in the vast majority of practical cases, such behaviour is not encountered; both schemes give the same results using a time-discretisation step of the order of ∆t ~ 10 -1 to 10 -2 ns.Such time-steps are usually used in interaction simulations of mm-wave gyrotrons.This agreement has led to a confidence in both schemes, which, following the results of the present study, needs to be reassessed as far as the BTCS scheme is concerned.The issue requires some special attention because the convergence of the BTCS with the time-discretisation step ∆t may be misleading.For example, in Figure 1, BTCS gives identical results for time-discretisation down to ∆t = 10 -4 ns, which is already two orders of magnitude below the "usual" time-discretisation.Without trying even smaller ∆t, one would believe that the BTCS scheme has already converged.
In [14] it was suggested that, for gyrotron time-dependent simulations, the BTCS may be preferable to the C-N on the grounds that C-N could be numerically unstable over long simulated time intervals.A discussion on this topic is beyond the scope of this paper.However, in the cases we have studied there are no traces of numerical instability for the C-N, i. e., no "exploding" solutions.Furthermore, from our experience using the C-N for a large variety of simulations, we have never faced instability issues.

Fig. 2 .Fig. 3 .
Fig. 2. Study of convergence, vs. the time discretisation step ∆t, for BTCS and C-N schemes.Test case: 1 MW, 170 GHz gyrotron (EU fall-back solution for ITER), ∆z = 0.1 mm, ∆t parameter.In the graphs, the absolute field profile of the operating TE 32,9 mode vs. the resonator axis is shown.

Fig. 4 .
Fig. 4. Study of convergence, vs. the time discretisation step ∆t, for BTCS and C-N schemes.Test case: 10 W, 263 GHz frequency-tunable gyrotron for DNP-NMR, ∆z = 0.05 mm, ∆t parameter.In the graphs, the output power of the operating TE 7,2 mode vs. time is shown.The magnetic field is increasing linearly with time