The phase-space structure of resonances between fast ions and an Alfvénic mode and the associated modification of density profiles in tokamaks are studied as a function of particle collisions. Guiding-center simulations in a realistic tokamak equilibrium are employed to address the resonance broadening parametric dependencies with respect to changes in the pitch-angle scattering rate. The rate of collisional replenishment, along with resonance strength, given by the combination of eigenmode and resonance structures and equilibrium parameters, determines saturation amplitudes for a given damping rate. As seen from the distribution function flattening, collisions have an effect of broadening the resonances, while the absolute value of δf decreases with increasing collisionality. It is observed that the collisional broadening can be comparable to the collisionless resonance width due to the mode amplitude alone. The resonance broadening coefficients are compared with the existing theory based on analytically expected saturation levels, showing fair agreement. The results can be useful in assisting reduced kinetic models, such as quasilinear models, when prescribing the effective resonance phase-space width, i.e., the mode-particle interaction platform, due to collisional or turbulent processes.

In large tokamaks, the partial pressure due to energetic particles (EPs) is comparable to the pressure due to the thermal plasma population. The EP pressure gradient is capable of severely destabilizing magnetohydrodynamic modes through resonant interactions, which poses concerns for the operation of ITER.1,2 In the absence of collisions, the separatrix of the phase-space resonant island structure, although possessing a small stochastic width, is a robust barrier that delimits resonant and non-resonant particles. Collisions can, however, redistribute particles across the separatrix and possibly alter the mode energy extraction and the scalings of phase mixing dynamics.3–5 In this paper, we employ guiding-center-following ORBIT7 simulations to address a basic question regarding the degree to which the Alfvénic mode particle resonance is effectively broadened upon the introduction of collisional scattering.

The width of an isolated linear resonance, as observed with a kinetic Poincaré plot, is proportional to the square root of the perturbation amplitude. If the mode structure does not change considerably within the resonance8 and the perturbation amplitude is small enough to prevent significant higher-order Fibonacci resonances, the dynamics of a resonant particle can be described by a nonlinear pendulum framework.9 In this case, the maximum resonance width in terms of frequency ΔΩ can be shown to be ΔΩ = 4ωb, where ωb is the bounce (or trapping) frequency.9,10 Since full simulation of the dynamics of fast ions resonating with Alfvén waves has proven to be costly, we wish to be able to model the kinetic interaction by means of reduced approaches, such as a quasilinear (QL) theory,11,12 where other unexplored mechanisms can as well contribute to the broadening. In particular, we are interested in further developing the physics basis of the Resonance Broadened Quasi-linear (RBQ) code.6,29

Dupree13 showed that background turbulent fields act to diffuse particle orbits in such a way as to increase the effective resonance width. In that case, it has been shown that the resonance platform can be written as a window function, expressed in an integral form. In the limit of zero turbulent-induced diffusion and small mode amplitude, this window function approaches the usual delta function. In principle, any mechanism that leads to resonant particle phase randomization such as pitch angle scattering15,16 should contribute to the broadening.

Along the lines of Dupree's reasoning of stochasticity-induced broadening, it has been suggested17 that the resonance width in QL theory should also scale with the effective pitch-angle scattering frequency experienced by the resonant population. The underlying reason is that scattering leads to non-resonant particles being kicked into resonance while kicking resonant particles out of resonance. Therefore, intuitively, scattering increases the effective resonance extension and allows for enhanced release of energy from the particle distribution.

The dimensionality reduction of the full kinetic phase space (involving actions and their canonically conjugated angles) to a reduced QL diffusive approach along the actions of the unperturbed motion18 requires an appreciable level of phase randomization in order to destroy phase space correlations and to average over the ballistic response. This randomization is often associated with the overlap of resonances described by the Chirikov criterion.19 A QL approach was proposed17,20 in order to cover both regimes of isolated and overlapping resonances. The former case is only within the assumptions required by the underlying QL applicability if collisions are strong enough to prevent particles from staying in phase with a resonance for several bounce times, which in practice forbids coherent nonlinear phenomena such as chirping and avalanches from taking place.14,16 The proposed approach, called the line broadened QL model,17,20 relies on tuning the constant coefficients a and c of the resonance broadening recipe

(1)

in such a way that the evolutions of single modes are forced to meet the expected saturation levels resulting from analytical theory3–5,21 in both limiting cases of very close to and very far from linear marginal stability. In this paper, we use a guiding center following code to examine at a more fundamental level whether the collisional broadening implies an actual broadening, in the sense of visibly extending the flattened region of the distribution function. In the presence of collisions, Poincaré plots become blurred and therefore cannot be used as a diagnostic tool for the resonance width. For that purpose, instead, we perform the analysis by means of density plots, using the flattening as an indication of the resonant extension.

The unperturbed orbits in an axisymmetric tokamak are uniquely determined by values of toroidal canonical momentum Pζ, energy E, and magnetic moment μ. The guiding center dynamics formalism we employ exploits the particle actions, which are invariants of the unperturbed motion, as the relevant variables to describe the effects of the perturbation. For the interaction of EPs with low-frequency Alfvénic modes, μ is approximately conserved and the other two actions are the toroidal and poloidal canonical momenta, defined as Pζ=gρψp and Pθ=ψ+ρI,7,22 respectively, where ψ and ψp are the toroidal and poloidal fluxes and dψ/dψp=q(ψp) is the field line helicity. g(ψp) and I(ψp) are the toroidal and poloidal covariant components of the equilibrium magnetic field, and ρ is defined as the parallel velocity divided by the strength of the magnetic field. Hamilton's equations are

(2)

where the poloidal and toroidal angles are θ and ζ. The Hamiltonian is H=ρ2B2/2+μB.

In order to examine the collisional effects on a resonance, we searched through a number of equilibria and mode spectra to find a case exhibiting a single isolated resonance,23,24 well surrounded by good Kolmogorov–Arnold–Moser (KAM) surfaces,25 existing within the co-passing particle domain. Cases with well isolated resonances allow the study of large scale islands without interaction with other nearby resonances which complicate the particle transport. We use a tokamak equilibrium with a field on axis of 4.9 kG. The major radius was 215 cm, and the minor radius was 62 cm. The equilibrium and q profile are shown in Fig. 1. The perturbation was an Alfvén mode with a simple broad Gaussian radial form. The particle distribution consisted of 85 keV to 98 keV deuterium ions.

FIG. 1.

Equilibrium and q profile using Boozer coordinates.

FIG. 1.

Equilibrium and q profile using Boozer coordinates.

Close modal

The Alfvén polarization mode is chosen with a given toroidal mode number n and with an eigenfrequency ω without loss of generality. This is because in the present study, we do not specify detailed plasma parameters but rather adjust them to investigate the resonant ion broadening near the isolated resonance. The Alfvén velocity is a function of the plasma density, and our choice for the toroidal and poloidal mode numbers is to have the mode located near the q =1 surface. Such modes as toroidal Alfvén eigenmodes (TAEs) and reversed shear Alfvén eigenmodes (RSAEs) are well suited to such conditions.1 

Assuming the Hamiltonian to be a function of ψp, θ, and nζωt, we have from Ṗζ=ζH and Ė=tH that there is a constant of the motion given by22 

(3)

simply related to the particle energy in the frame rotating with the mode. The kinetic Poincaré recurrence plots are produced at constant values of μ and K. Those plots involve either Pζ or E as a function of θ on the orbit. Each point in the Poincaré plot is obtained when the recurrence condition nζωt=2πk,k integer, is satisfied.

A passing population was chosen with all particles having the same value of magnetic moment μ and K. The equations of motion are easily generalized to include flute-like perturbations of the form δB=×αB, with B being the equilibrium field and α=m,nαm,n(ψp)sin(nζmθωnt). The perturbation α has units of a length, simply related to the cross field ideal displacement produced by the mode. A well isolated resonance was found for μB0= 50 keV, with B0 being the on-axis field strength, with a perturbation with a single harmonic with ratio between poloidal and toroidal mode numbers of m/n =6/5 and a frequency of 130 kHz. The plane of E and Pζ in Fig. 2 shows a radially isolated resonance extending from 80 keV upwards in the co-passing domain, and the particle distribution is shown as a line at fixed K intersecting the resonance approximately at the mid radius, for Pζ = 0.31ψw and E =93 keV. The mode can move particles only along this line, not across it. Amplitude A is the maximum value of the ideal eigenfunction of the perturbation, normalized to the major radius of 210 cm. For a discussion of plots of this type, see Ref. 23.

FIG. 2.

The plane of E, Pζ for μB0 = 50 keV, showing an isolated resonance in the co-passing domain, extending from 80 keV upwards, for amplitude A =4 × 10−5, with a frequency of 130 kHz, and poloidal and toroidal mode numbers m/n =6/5. The red squares in a nearly vertical strip in the center are the resonance, and the particle distribution is the near horizontal line cutting across it near the plasma center. The right edge of this line is the magnetic axis, and the left end is near the plasma edge. ψw is the poloidal flux at the plasma edge.

FIG. 2.

The plane of E, Pζ for μB0 = 50 keV, showing an isolated resonance in the co-passing domain, extending from 80 keV upwards, for amplitude A =4 × 10−5, with a frequency of 130 kHz, and poloidal and toroidal mode numbers m/n =6/5. The red squares in a nearly vertical strip in the center are the resonance, and the particle distribution is the near horizontal line cutting across it near the plasma center. The right edge of this line is the magnetic axis, and the left end is near the plasma edge. ψw is the poloidal flux at the plasma edge.

Close modal

A density gradient in energy E or Pζ drives the mode amplitude. Particle rotation within the resonance at the bounce frequency partially flattens the distribution within the island, a process which delivers energy to the mode, causing growth. An example of this flattening is shown in Fig. 3. The resonance is at Pζ = 0.31ψw, where ψw is the poloidal flux at the last closed flux surface. Without particle collisions to replenish the density gradient within the resonance, mode growth ceases within roughly one average particle bounce time within the resonance. Mode saturation occurs when the rate at which the resonance is flattened, represented by the internal rotation (mixing) frequency ωb, becomes linearly proportional to the effective frequency of collisional replenishment of the density gradient, represented by the effective scattering frequency νeff, which was introduced theoretically to combine the effect of Coulomb collisions and the rate of particle decorrelation from the resonance location coming from other stochastic mechanisms such as turbulent scattering and RF diffusion. The exact proportionality constant between the two is determined by the initial linear drive and background damping rates.3–5,21,26,27

FIG. 3.

The flattening of the distribution due to a resonance without collisions. The positive gradient in Pζ, corresponding to a negative gradient in the minor radius, is destabilizing. The resonance is located at Pζ = 0.31ψw, and the mode amplitude is A =2 × 10−4.

FIG. 3.

The flattening of the distribution due to a resonance without collisions. The positive gradient in Pζ, corresponding to a negative gradient in the minor radius, is destabilizing. The resonance is located at Pζ = 0.31ψw, and the mode amplitude is A =2 × 10−4.

Close modal

The dynamics of this process can be observed by launching particles within the resonance, both with and without collisions, at single values of E, μ, and Pζ, but with different values of θ and ζ. In the unperturbed equilibrium, all these particles describe the same orbit—they are distinguished only by their phase relations with respect to the mode. The position of the initial particle launch to determine the mean internal rotation and the time for local diffusion to replenish the density gradient in the resonance is at E = 93 keV and Pζ/ψw = 0.31.

Without collisional scattering, one observes a simple damped oscillation, as the initial distribution is mixed by the internal rotation, and the energy and Pζ spread out to the bounding KAM surfaces. The oscillation can be observed through either the mean square energy or the mean modification of Pζ2 of the particle distribution. We choose to examine dPζ2, but these variables are simply related through the conservation of K. The magnitude of the final distribution width gives the mean width of the resonance. Since the rotation rate is a function of the distance between the O-point and the separatrix, dropping to zero at the separatrix (where the bounce period diverges logarithmically), as time goes on the mixing occurs at smaller and smaller scale lengths, and becomes an irreversible process, increasing the entropy of the system due to dissipation.

The collisionless simulation determines the mean bounce time of particles in the resonance and the width of the resonance. The collisional simulation determines the time it takes for particles to diffuse across the resonance, giving the time for replenishment of the local distribution. Note that the diffusion rate is given by the motion away from the isolated resonance, and its collisional diffusion across the surrounding good KAM surfaces. Collisions in the ORBIT code are given by energy conserving pitch angle scattering, in a Monte Carlo representation.30 

Figure 4 presents a determination of the effective mode width and the local diffusion time. The unit of time T in these plots is the unperturbed toroidal transit time of a particle near the major axis, equal to 4.7 μs in this case. Replenishment of the density gradient within the island is provided by collisions, with a characteristic time given by the time for nearby particles to diffuse across the width of the resonance Td, and is given by the square of the resonance width divided by the slope of the collisional plot S, which in this case is 10 transits or 47 μs. The period of rotation in the resonance, approximately given by the first crossing of the collisionless plot with its asymptotic value, is also 10 transits, indicating that this mode amplitude corresponds approximately to a saturation amplitude for this value of the collision frequency.

FIG. 4.

Collisionless evolution and collisional diffusion, A =2 × 10−4, ν=104/T,dPζ2=7.8×104,andS=8×105, with T being the on axis particle transit time. S is the slope of the diffusive plot. The mixing time Tm is given by the first crossing of the collisionless curve with the asymptotic value, Tm = 10, in units of toroidal transits. Td, the time to diffuse across the resonance, is given by dPζ2/S.

FIG. 4.

Collisionless evolution and collisional diffusion, A =2 × 10−4, ν=104/T,dPζ2=7.8×104,andS=8×105, with T being the on axis particle transit time. S is the slope of the diffusive plot. The mixing time Tm is given by the first crossing of the collisionless curve with the asymptotic value, Tm = 10, in units of toroidal transits. Td, the time to diffuse across the resonance, is given by dPζ2/S.

Close modal

For a fixed value of collision frequency ν, a plot of the mixing time Tm and the diffusion time Td versus mode amplitude A shows a crossing of these values. This occurs because for the increasing amplitude the mixing time decreases, whereas because of the increasing width of the resonance, the time to diffuse across the resonance increases. In Fig. 5 is shown an estimate of the saturation value using this method. Note that this estimation is best applicable for the strongly driven scenario as it ignores mode damping, which can modify the long time behavior and the saturation value. In the fast excitation regime, the mode is expected to saturate not far from the level instantaneously reached after the first bounce, i.e., on a timescale faster than those associated with damping and collisional effects. Large values of damping, i.e., in the marginal stability case, can produce very different behavior, with large amplitude mode oscillations and possibly intermittency,4–7 which is expected to break down the validity of the saturation level estimate provided in Fig. 5. In addition, we observe that in the marginally stable case, the slope of the distribution function is very slightly changed with respect to the equilibrium one. This prevents the identification of the resonance extension, and the marginal regime is not explored in the remainder of this paper.

FIG. 5.

An example of the estimation of the saturation amplitude A. A plot of the mixing period Tm and the diffusion period Td vs A, showing intersection giving an estimation of the saturation amplitude. Collisions are fixed at ν=104/T for each value of A.

FIG. 5.

An example of the estimation of the saturation amplitude A. A plot of the mixing period Tm and the diffusion period Td vs A, showing intersection giving an estimation of the saturation amplitude. Collisions are fixed at ν=104/T for each value of A.

Close modal

The scaling of the saturation level with collisions can be determined by equating Tm and Td. Since Tm1/ωb and ωbA, we have Tm1/A. The time to diffuse across the resonance is given by Td=dPζ2/S, with S being the slope of the diffusion curve and the resonance width dPζA. The diffusion across the nearby unperturbed KAM surfaces is given by dPζ2νt giving Sν. Thus, setting Tm = Td gives Asatν2/3. The linear dependence on A of Td and the 1/A dependence of Tm can be seen in Fig. 5.

For each value of the collision frequency, a number of determinations of Tm and Td are performed in order to find an estimate for the saturation amplitude. In Fig. 6 are shown the values of Tm and A vs collision frequency ν for the values at saturation given by this method. We see that Tmν1/3, giving ωbν1/3 and Aν2/3. Thus, dynamically, the effect of collisions is to increase the saturation level by supplying the drive given by the density gradient.

FIG. 6.

Values at saturation. Tm vs ν, showing Tmν1/3, and A vs ν, showing Aν2/3. The red lines are Tm=0.57ν1/3 and A=0.36ν2/3.

FIG. 6.

Values at saturation. Tm vs ν, showing Tmν1/3, and A vs ν, showing Aν2/3. The red lines are Tm=0.57ν1/3 and A=0.36ν2/3.

Close modal

The effect of particle collisions on resonances has been the subject of a number of theoretical investigations, and in particular, it has been suggested that the collisions effectively broaden the resonance by allowing the communication of particles within and near a resonance.21,29

To examine the effect of collisions on the resonance associated with a fixed amplitude perturbation, a particle distribution is initiated with a range of ψp, energy, and Pζ, all with the same values of μ. Particles are launched along a value of constant K, with a strong density gradient in Pζ, corresponding to a gradient in the minor radius, and in energy E. The profile is produced using standard Monte Carlo methods. We examine the modification of the resonance width in Pζ.

In Fig. 7 are shown the Poincaré plots for the two amplitudes used for our examination of profile flattening. The mode was chosen to produce an isolated resonance near the center of the distribution with wide ranges of unperturbed KAM surfaces on each side of the resonance.

FIG. 7.

Poincare plots for A = 10−4 and A = 2 × 10−4. The mode frequency is 130 kHz, with m/n = 6/5. The mode is chosen to have a large region of relatively undisturbed KAM surfaces on each side of the resonance.

FIG. 7.

Poincare plots for A = 10−4 and A = 2 × 10−4. The mode frequency is 130 kHz, with m/n = 6/5. The mode is chosen to have a large region of relatively undisturbed KAM surfaces on each side of the resonance.

Close modal

The collisions are energy conserving pitch angle scattering and directly change the value of Pζ through ρ. The gradient must be maintained in time by re-injecting lost particles using the Monte Carlo procedure. The action of the resonance produces a modification of the particle density n(Pζ,t), consisting of a local flattening in the vicinity of the resonance at Pζ = Pres as shown in Fig. 3. For small amplitude, the profile is maintained in the vicinity of the resonance hyperbolic point, and so, the flattening is only partial. For larger amplitudes, the stochastic layer about the separatrix is sufficiently large to destroy KAM surfaces around the hyperbolic point, and so, a more complete flattening can be observed. To study this flattening, we perform a time average of the density given by the particle distribution, using a coarse graining in the variable Pζ, in order to cancel out the periodic density oscillations occurring on undestroyed KAM surfaces and to improve the statistics of the data.

The modifications of the density profile are shown in Fig. 8. Six hundred thousand particles are launched using standard Monte Carlo methods to produce the initial density profile. After a delay to allow phase mixing, a time average of the density using five hundred toroidal transits is taken in each of the 50 bins in Pζ, to reduce statistical noise. This simulation is difficult because of limitations in the mode amplitude and collisionality. If the amplitude is too small, the flattening is weak and statistically difficult to observe accurately. If the amplitude is too large, the large island and associated nearby nonlinear resonances cause global chaos and large scale particle loss. If the collision rate is too large, the density profile is strongly modified, making the observation of the local flattening impossible. We choose to show δf, the difference between the distribution with and without the presence of the mode, giving a much clearer picture of the effect of the resonance. However, for each value of the collision rate, the density profile is modified in time also in the absence of the mode, and so, a determination of the collisional mode free density profile must be made using the same time interval as used in the presence of the mode and the unperturbed density profile subtracted from the perturbed one to give δf. Note that in these plots, density flattening appears as a negative slope. Thus, with no collisions, the negative slope occurs only between the edges of the Poincaré resonance. Outside, on the two sides of the resonance, the slope of δf is in fact positive. Collisionless flattening of the density profile occurs only within the island structure, due to the bounce frequency mixing. Outside the resonance island, the KAM surfaces are distorted outward, as seen in Fig. 7, leading to a local increase in the density, visible as a positive slope in the plot of δf in Fig. 8. This increase is adiabatic since there is no mixing occurring in these distorted surfaces, and the density and particle energy will return to their initial values upon decay of the mode.

FIG. 8.

Left: Density profile modification with perturbation A = 10−4. The collisionality values shown are ν=0 (black), 0.6 Hz (blue), 6 Hz (red), and 22 Hz (green). Right: For A =2 × 10−4 shown are ν=0 (black), 2 Hz (blue), 15 Hz (red), and 22 Hz (green).

FIG. 8.

Left: Density profile modification with perturbation A = 10−4. The collisionality values shown are ν=0 (black), 0.6 Hz (blue), 6 Hz (red), and 22 Hz (green). Right: For A =2 × 10−4 shown are ν=0 (black), 2 Hz (blue), 15 Hz (red), and 22 Hz (green).

Close modal

The collisionality values shown are, for A =10−4, ν=0 (black), 0.6 Hz (blue), 6 Hz (red), and 22 Hz (green). For A =2 × 10−4 shown are ν=0 (black), 2 Hz (blue), 15 Hz (red), and 22 Hz (green). The collisions are seen to initially decrease the depth of the density modification but to extend the width.

Note that these collision rates are all below the values needed for mode saturation at these amplitudes, as can be seen in Fig. 6, but the final value of 22 Hz is very near the saturation value for these amplitudes. Also note that the mean bounce time in the resonance is about 10 transit times, but the collision time is more than 1000 times this.

The collisionless resonance widths and the bounce frequencies of particles well trapped in the resonance are shown in Figs. 9 and 10. There is some variation of the width from island to island, but these figures include the separatrices, and so, a reasonably accurate determination of the width can be made.

FIG. 9.

Left—Poincare plot showing a particle trapped in the resonance for amplitude A =10−4. Right—Time evolution of the trapped bounce motion for A =10−4.

FIG. 9.

Left—Poincare plot showing a particle trapped in the resonance for amplitude A =10−4. Right—Time evolution of the trapped bounce motion for A =10−4.

Close modal
FIG. 10.

Left—Poincare plot showing a particle trapped in the resonance for amplitude A =2 × 10−4. Right—Time evolution of the trapped bounce motion for A =2 × 10−4.

FIG. 10.

Left—Poincare plot showing a particle trapped in the resonance for amplitude A =2 × 10−4. Right—Time evolution of the trapped bounce motion for A =2 × 10−4.

Close modal

In Fig. 11, it is shown the scaling of the resonance broadening with collision frequency. The calculated estimates of the broadening involve processing output data, which naturally incur errors. Numerically, the error bars for ΔPζ/ψw come from two main sources, namely, the visual inferences of the broadening from the perturbed distribution (from Fig. 8) and of ωb, which is calculated from the spacing of two consecutive peaks in Figs. 9 and 10. Therefore, the propagation of uncertainty can be estimated as

where (ΔPζ/ψw)(ωb)(ΔPζ/ψw)(ΔΩ)(ΔΩ)ωb4(Pζ/ψw)Ω, which is then calculated from Fig. 12. In both Figs. 8 and 11, note that the curves do not go to zero at ν=0 but approach the collisionless width (b). The value (ΔPζ/ψw)(Δf) is extracted from Fig. 8 and changes according to the level of collisionality.

FIG. 11.

Scaling of the collisional width with ν, for A =10−4 (squares) and A =2 × 10−4 (disks). The lines are fits to the data using ν1/3 for each amplitude. To avoid ambiguities in the definition of the effective collisional broadening, we take as a reference the location of the peaks of δf, inside which the gradient is negative, i.e., where phase mixing due to inverse Landau damping takes place.

FIG. 11.

Scaling of the collisional width with ν, for A =10−4 (squares) and A =2 × 10−4 (disks). The lines are fits to the data using ν1/3 for each amplitude. To avoid ambiguities in the definition of the effective collisional broadening, we take as a reference the location of the peaks of δf, inside which the gradient is negative, i.e., where phase mixing due to inverse Landau damping takes place.

Close modal
FIG. 12.

Plot of Ω=nωζpωθ versus Pζ for fixed K, where ωζ is the mean value of Δζt and ωθ is the mean value of Δθt for particles launched along the line K=ωPζnE, all with μB = 50 keV and p =4 which is the number of islands in the resonance, see Fig. 7.

FIG. 12.

Plot of Ω=nωζpωθ versus Pζ for fixed K, where ωζ is the mean value of Δζt and ωθ is the mean value of Δθt for particles launched along the line K=ωPζnE, all with μB = 50 keV and p =4 which is the number of islands in the resonance, see Fig. 7.

Close modal

We now fit the results of the simulations to existing theory for the collisional broadening of resonances as well as possible. The collision frequency used was ν=15 Hz, giving ν=94rad/s. From Eq. (6) of Ref. 14, we have

(4)

The sub-index that accompanies the derivative indicates that it is taken along a path that preserves the values of K [Eq. (3)] and μ at which a particle is launched. Given the constancy of K and μ for particles resonating with a mode whose frequency is much smaller than the cyclotron frequency, it is convenient to project the resonant dynamics onto the remaining free variable. The definition of the effective scattering νeff is useful in order to express the kinetic equation in terms of the relevant variable Ω(E,Pφ,μ) describing the one-dimensional resonant particle motion, for which case the scattering term can be expressed as νeff32f/Ω2 instead of in three-dimensional velocity space. We note that, in addition to collisions, a number of stochastic effects that influence resonant particle dynamics can be incorporated into νeff, such as turbulence, RF heating, and energy diffusion.15,16

From Fig. 12 

(5)

and using μB = 50 keV, R=215cm,andBpol/B=1/4, we have νeff=1.63×105rad/s for the collision frequency used. We find νeffTrun=600 for the largest value of collision frequency used and 20 for the smallest, with Trun being the length of the simulation, giving time for equilibration of the collisions with the resonance.

From Fig. 9, we have for A =10−4 the bounce period to be 0.1 ms, giving ωb = 6.28 × 104 rad/s. The resonance for A =10−4 has a width of δPζ/ψw=0.1. For the collisionless frequency width of the resonance, we have using Eq. (1), ΔΩ=aωb=2.5×105rad/s, giving a =3.9 ± 0.1.

Including collisions, using Eq. (1) and the total width of the flattening δPζ/ψw=0.36, and again using Eq. (5), we have for the frequency width including collisions ΔΩ=9×105rad/s, giving c =4 ± 0.5, with the error determined from Fig. 11.

From Fig. 10, we have for A =2 × 10−4 the bounce period to be 0.06 ms, giving ωb = 1.04 × 105 rad/s. The resonance for A =2 × 10−4 has a width of δPζ/ψw=0.15. For the collisionless frequency width, we have ΔΩ=aωb=3.7×105rad/s, and using Eq. (5), we find a =3.7 ± 0.2.

Including collisions, using Eq. (1) and the total width of the flattening δPζ/ψw=0.45, and again using Eq. (5), we have ΔΩ=9.×105rad/s giving c =3.2 ± 0.4, with the error determined from Fig. 11.

Reference 11 reported a systematic study of the value of a with the changing mode amplitude. It has been found that the maximum nominal value a =4 is only achieved when the mode amplitude is small enough in order to ensure that nonlinear resonances, island coupling effects, and mode structure non-uniformities do not play an appreciable role in modifying the resonance structure. As the amplitude is increased, a decreases. We note that the values obtained for a presented here are in agreement with those reported in Ref. 11. We employ a broad eigenstructure that introduces a minor degree of island asymmetry. As the amplitude is increased from A =10−4 to A = 2 × 10−4, a transitions from 3.9 to 3.7.

If to define the modification of the density profile, the full extent of the modification of f is used instead of the flattening, then both a and c are approximately doubled, to a =8 and c =6.

In the RBQ code,6,29 the coefficient c is calculated exploring the situation in which the system reaches a steady saturation. In this limiting case, a previously known saturation level9 can be exploited with the purpose of obtaining the value of c that would need to be used in the quasilinear system to replicate the expected saturation level. At saturation, ∂f/∂t =0 and the broadened quasilinear diffusion equation implies ΔΩ=π2ωb4νeff3γdγLγd,6 where γL is the initial linear growth rate of the mode in the absence of damping and γd is the mode background damping rate. Close to marginal stability, γLγd, the expected saturation level is ωb1.18[1γdγL]1/4νeff,9 which implies ωbνeff. Substituting the broadening recipe, Eq. (1), one finds c =1.184π/2 = 3.05. We note that in Ref. 6, the value of c was erroneously typed as 2.7 although the calculation that led to that value was correct. The results reported here for the collisional broadening coefficient c show a remarkable similarity, within the error bars, with the value of c that follows from the saturation level embedded in quasilinear asymptotic quasi-steady. We note that the contribution to the broadening due to the growth rate is neglected in RBQ1D6,29 to account for Kaufman's revision of the QL theory.31 

Consider the resonance dynamics being described by the variable Ω=ξ̇, where ξ is the angle associated with the system action. The resonance condition can be expressed as Ω=ω, and the deeply trapped particles satisfy the pendulum equation ξ̈+ωb2sin(ξωtξ0)=0.10 The resonance island is parameterized by Ω2=2ωb2[1+cos(ξωtξ0)].

We note that there is a broadening conversion factor between the full nonlinear and the quasilinear methodologies. In quasilinear theory, the angular dynamics is not resolved, which implies that the (collisionless) resonance interaction only takes place within a rectangle of width ΔΩ=aωb in action space and a width of 2π in the generalized angle. In ORBIT, the resonance platform has an elliptic parametrization. Considering δf=Ωf0Ω|Ωres, the momentum exchange in the full kinetic framework can be calculated as

(6)

In the quasilinear framework,

(7)

equating the two expressions, it follows that a quasilinear broadening of a =3 (constant for all angles) replicates the same momentum exchange with respect to an elliptic island in kinetic theory that has a maximum width of a =4 (at the elliptic point) and a minimum width of a =0 (at the hyperbolic points). One, therefore, needs to have the conversion factor of ∼0.75 in mind when extrapolating ORBIT guiding center simulations to reduce quasilinear models.

We wish to parameterize the form of the modification of the resonance due to collisions since a simple functional form is useful for further studies. The displacement of particles due to the resonance, with or without collisions, is conservative, and so, we need dPζδf=0 in either case. To this end, and with no theoretical basis other than that diffusion from a point source takes the form of a Gaussian, we fit the data for δf with the derivative of a Gaussian

(8)

with p=(PζP0)/ψw. The extrema of this function appear at p = ±W, with the local maxima (minima) being δm.

All the parameters appearing in this representation are functions of the mode amplitude A and the collision frequency ν. The width is given by W=W0(A)+W1(A)ν1/3 as seen in Fig. 11. The fit of this representation of δf is shown in Fig. 13 for an amplitude of 2 × 10−4, for no collisions and for the maximum collision frequency used, ν=22 Hz. The fit is best for these cases, and there is a minimum of noise in the data when the collisions are either absent or strong enough to smooth the distribution. Note that the center of δf, P0, shifts leftward (down the density gradient) with the increasing collision rate. Clearly, the real situation is more complicated than what is given by this simple form; undoubtedly, the actual Poincaré particle resonance with its original width is still somewhat present, although distorted by the collisions, and produces modifications of the distribution function near the resonance edges, where the collisions are unable to completely smooth out. Nevertheless, this function captures the major shape of the profile modification, in particular, when the collision rate is near that corresponding to the mode saturation level. The behavior of the maximum amplitude δm is not simple, dropping as a function of collision rate from the collisionless value to a minimum and then increasing, shown in Fig. 14. Note that the maximum value of the collision rate in this plot is very near the expected rate for saturation at this amplitude and that the maximum value of δf is more than half the original maximum without collisions.

FIG. 13.

Fit of δf for an amplitude of 2 × 10−4. Data with no collisions are in black, with black triangles. Data with collisions are in red, with red squares. The analytical fit given by Eq. (8) for each case is shown in green. The collision rate was ν=22 Hz.

FIG. 13.

Fit of δf for an amplitude of 2 × 10−4. Data with no collisions are in black, with black triangles. Data with collisions are in red, with red squares. The analytical fit given by Eq. (8) for each case is shown in green. The collision rate was ν=22 Hz.

Close modal
FIG. 14.

Profile maxima δm for mode amplitudes of 10−4 (squares) and 2 × 10−4 (triangles) and flattening as a function of collision rate.

FIG. 14.

Profile maxima δm for mode amplitudes of 10−4 (squares) and 2 × 10−4 (triangles) and flattening as a function of collision rate.

Close modal

Although the collisions produce significant broadening of the resonance, and for a rapid collision rate, the maximum value of δf is comparable to that with no collisions, the relative flattening of the profile over the extended structure is not nearly as strong as that given by the initial collisionless resonance. In Fig. 14, we see that the initial collisionless flattening, equal to δm/W, is about 0.8, but this value drops to about 0.1 as the resonance broadens with the increasing collision rate. Interestingly, this behavior appears to be only weakly dependent on the mode amplitude, at least within the limited range of values considered.

The resonance broadening of a high energy particle distribution with an ideal MHD mode is examined using guiding center analysis, including the effect of collisions. Collisional broadening of a resonance is directly observed, with the flattening scaling as ν1/3 and fit to the analytical expression ΔΩ=aωb+cνeff, with a ≃ 4 and c ≃ 3–4. Although our investigation is limited to a specific resonance, the theoretical underpinning of the broadening expression indicates that this is a universal result, applicable in general to other resonances, and useful for models of the effect of resonances on high energy particle distributions.

The visual identification of the resonance width can be challenging when more than one peak appears in δf. Numerical noise arising from the stochastic layer that forms around the island and resonance asymmetries also contribute to limiting the accuracy with which one can extract the coefficient c. The results, however, show undoubtedly that the collisional broadening can be comparable and even exceed the collisionless broadening. Figure 11 shows that the scaling ΔΩνeffν1/3 proposed in Refs. 6 and 28 is consistent with the present results.

The authors thank J. A. Krommes for enlightening discussions. This work was supported by the U.S. Department of Energy (DOE) under Contract No. DE-AC02-09CH11466.

1.
N. N.
Gorelenkov
,
S. D.
Pinches
, and
K.
Toi
,
Nucl. Fusion
54
,
125001
(
2014
).
2.
R. B.
White
,
N.
Gorelenkov
,
W. W.
Heidbrink
, and
M. A.
Van Zeeland
,
Phys. Plasmas
17
(
5
),
056107
(
2010
).
3.
H. L.
Berk
and
B. N.
Breizman
,
Phys Fluids B
2
,
2226
(
1990
).
4.
H. L.
Berk
and
B. N.
Breizman
,
Phys Fluids B
2
,
2235
(
1990
).
5.
H. L.
Berk
and
B. N.
Breizman
,
Phys Fluids B
2
,
2246
(
1990
).
6.
N.
Gorelenkov
,
V. N.
Duarte
,
M.
Podesta
, and
H. L.
Berk
,
Nucl. Fusion
58
,
082016
(
2018
).
7.
R. B.
White
and
M. S.
Chance
,
Phys. Fluids
27
,
2455
(
1984
).
8.
R. B.
White
,
N. N.
Gorelenkov
,
V. N.
Duarte
, and
H. L.
Berk
,
Phys. Plasmas
25
,
102504
(
2018
).
9.
H. L.
Berk
,
B. N.
Breizman
, and
M.
Pekker
,
Plasma Phys. Rep.
23
,
778
(
1997
).
10.
G.
Meng
,
N. N.
Gorelenkov
,
V. N.
Duarte
,
H. L.
Berk
,
R. B.
White
, and
X.
Wang
,
Nucl. Fusion
58
,
082017
(
2018
).
11.
W.
Drummond
and
D.
Pines
,
Nucl. Fusion
2
,
1049
(
1962
).
12.
A. A.
Vedenov
,
E. P.
Velikhov
, and
R. Z.
Sagdeev
,
Sov. Phys. Usp.
4
,
332
(
1961
).
13.
T. H.
Dupree
,
Phys. Fluids
9
,
1773
(
1966
).
14.
V. N.
Duarte
,
H. L.
Berk
,
N. N.
Gorelenkov
,
W. W.
Heidbrink
,
G. J.
Kramer
,
R.
Nazikian
,
D. C.
Pace
,
M.
Podestà
, and
M. A.
Van Zeeland
,
Phys. Plasmas
24
,
122508
(
2017
).
15.
J.
Lang
and
G.-Y.
Fu
,
Phys. Plasmas
18
(
5
),
055902
(
2011
).
16.
V. N.
Duarte
,
H. L.
Berk
,
N. N.
Gorelenkov
,
W. W.
Heidbrink
,
G. J.
Kramer
,
R.
Nazikian
,
D. C.
Pace
,
M.
Podestà
,
B. J.
Tobias
, and
M. A.
Van Zeeland
,
Nucl. Fusion
57
,
054001
(
2017
).
17.
H. L.
Berk
,
B. N.
Breizman
,
J.
Fitzpatrick
, and
H. V.
Wong
,
Nucl. Fusion
35
,
1661
(
1995
).
18.
A. N.
Kaufman
,
Phys. Fluids
15
(
6
),
1063
1069
(
1972
).
19.
B. V.
Chirikov
,
Phys. Rep.
52
,
263
379
(
1979
).
20.
H. L.
Berk
,
B. N.
Breizman
,
J.
Fitzpatrick
,
M. S.
Pekker
,
H. V.
Wong
, and
K. L.
Wong
,
Phys. Plasmas
3
,
1827
(
1996
).
21.
V. N.
Duarte
and
N. N.
Gorelenkov
,
Nucl. Fusion
59
,
044003
(
2019
).
22.
R. B.
White
,
The Theory of Toroidally Confined Plasmas
, 3rd ed. (
Imperial College Press
,
2014
).
23.
R. B.
White
,
Comm. Nonlinear Sci. Numer. Simul.
17
(
5
),
2200
(
2012
).
24.
R. B.
White
,
Plasma Phys. Controlled Fusion
53
,
085018
(
2011
).
25.
A. N.
Kolmogorov
, in
Proceedings International Congress Mathematicians, Amsterdam
(
1957
), Vol.
1
,
315
;
V. I.
Arnold
,
Russ. Math. Surv.
18
(
5
),
9
(
1963
);
J.
Moser
,
Nachr. Akad. Wiss. Göttingen
II
,
1
20
(
1962
).
26.
M.
Zhou
and
R. B.
White
,
Plasma Phys. Controlled Fusion
58
,
125006
(
2016
).
27.
N. N.
Gorelenkov
,
Y.
Chen
,
R. B.
White
, and
H. L.
Berk
,
Phys. Plasmas
6
,
629
(
1999
).
28.
K.
Ghantous
,
H. L.
Berk
, and
N. N.
Gorelenkov
,
Phys. Plasmas
21
,
032119
(
2014
).
29.
N. N.
Gorelenkov
,
V. N.
Duarte
,
C. S.
Collins
,
M.
Podestà
, and
R. B.
White
, “
Verification and application of Resonance Broadened Quasi-linear (RBQ) model with multiple Alfvenic instabilities
,”
Phys. Plasmas
(submitted).
30.
A. H.
Boozer
and
G.
Kuo-Petravic
,
Phys. Fluids
24
,
851
(
1981
).
31.
A. N.
Kaufman
,
J. Plasma Phys.
8
(
1
),
1
5
(
1972
).