Contributions to the linear and nonlinear theory of the beam–plasma interaction

Abstract

We focus our attention on some relevant aspects of the beam–plasma instability in order to refine some features of the linear and nonlinear dynamics. After a re-analysis of the Poisson equation and of the assumption dealing with the background plasma in the form of a linear dielectric, we study the non-perturbative properties of the linear dispersion relation, showing the necessity for a better characterization of the mode growth rate in those flat regions of the distribution function where the Landau formula is no longer predictive. We then upgrade the original $N$-body approach in O'Neil et al. (Phys. Fluids, vol. 14, 1971, pp. 1204–1212), in order to include a return current in the background plasma. This correction term is responsible for smaller saturation levels and growth rates of the Langmuir modes, as result of the energy density transferred to the plasma via the return current. Finally, we include friction effects, as those due to the collective influence of all the plasma charges on the motion of the beam particles. The resulting force induces a progressive resonance detuning, because particles are losing energy and decreasing their velocity. This friction phenomenon gives rise to a deformation of the distribution function, associated with a significant growth of the less energetic particle population. The merit of this work is to show how a fine analysis of the beam–plasma instability outlines a number of subtleties about the linear, intermediate and late dynamics which can be of relevance when such a system is addressed as a paradigm to describe relevant nonlinear wave–particle phenomena (Chen & Zonca, Rev. Mod. Phys., vol. 88, 2016, 015008).

1. Introduction

A fundamental model of plasma physics consists of the beam–plasma instability, originally studied in O'Neil & Malmberg (1968), O'Neil, Winfrey & Malmberg (1971), Mynick & Kaufman (1978) and Tennyson, Meiss & Morrison (1994) (see also Elskens & Escande 2003), mainly having in mind real experiments, while it is currently a scenario to mimic some of the most relevant features of the interaction between fast ions and the Alfvén spectrum present in tokamak devices (Berk & Breizman 1990a,b,c; Breizman & Sharapov 2011).

The basic mechanism underlying the beam–plasma instability is the inverse Landau damping mechanism (O'Neil 1965; Shapiro & Shevchenko 1971; Lifshitz & Pitaevskii 1976), according to which fast particles are able to pump the Langmuir spectrum, up to a saturation limit, resulting in trapping into the instantaneous potential well. The original treatment of the problem in O'Neil et al. (1971) is particularly successful because it reduces the dynamics to a Hamiltonian $N$-body system, using dimensionless universal units. Such a model relies on the hypothesis that the beam propagates in a cold background plasma, feeling its response in the form of a real dielectric and treating the system as collisionless, de facto neglecting any other reciprocal back-reaction among these two components. Furthermore, the particle beam is assumed to be tenuous, i.e. the ratio between the beam and the plasma densities is a small control parameter of the interaction. However, as discussed in § 3.1, when the beam–plasma system (BPS) (say also the bump-on-tail instability) is implemented to describe the outgoing radial transport of fast ions interacting with Alfvén waves (for a precise mapping between the velocity to the radial space, see Carlevaro et al. (2016b, 2019)), some of the approximations underlying the original model could result in being inadequate.

In this paper, we start by providing, in § 2, a detailed derivation of the Poisson equation at the ground of the dielectric approximation for the background plasma, based on a careful investigation of the instantaneous linear response of the bulk when it is crossed by fast particles. Although the main ideas of the present derivation were implicitly contained in O'Neil et al. (1971), nonetheless the present contribution explicitly clarifies which approximations are required to summarize the plasma response in terms of a dielectric function. Such an analysis is of interest in order to better clarify the approach developed in § 5, where the thermal plasma back-reaction is described via fluid dynamics, and therefore the emergence of a return current can be described.

Then, after providing in § 3 the standard $N$-body description of the BPS, in § 4 we investigate the details of the linear dispersion relation in the case of an initial positive slope of the beam distribution in the velocity space. In particular, we address non-perturbative effects and the breakdown of the inverse Landau damping expression where the ratio between the growth rate and the real frequency shift with respect to the plasma frequency becomes finite.

In § 5, we analyse the problem concerning the induction of a return current within the background plasma, by coupling the $N$-body dynamics with the electrostatic fluctuations in the plasma, described via a linear fluid representation. We arrive at a revised model, in which the Poisson equation has memory of the presence of the Langmuir oscillations. The numerical implementations of this scheme, both for a single mode and for the so-called quasi-linear model (when many modes are included in a broad spectrum), are developed. As a result, we demonstrate that the spectrum and the beam distribution function morphology are similar to those predicted by the original treatment. However, now it is possible to calculate a return current in the background plasma, whose presence depresses the saturated spectrum energy and the mode linear growth rates. This is of clear importance in view of the realistic prediction for relaxed fast-particle profiles.

Finally, in § 6 we investigate the effect of friction on the motion of fast particles, in order to clarify how their distribution function and the spectrum are affected in the quasi-linear scenario. The friction we are addressing consists of the interaction of the fast particles with the bulk of the plasma charges, as integrated at a distance greater than the Debye scale. In other words, we neglect the effects of collisions along the fast-particle motion, because of their extreme dilution, but we account for collective interaction, acting as a force contrasting the particle dynamics, in analogy to what has already been studied in the literature for heavy ions (Neufeld & Ritchie 1955). We find the interesting feature of a resonance detuning, able to significantly deform the particle distribution function and the mode spectrum on a given time of the nonlinear dynamics. This effect is due to the loss of energy that the particles undergo under the collective friction of the plasma charges. Actually, a rather efficient displacement of particles takes place in the phase space towards the low-velocity region, which depresses those resonant modes at smaller wavelength. Clearly, after a sufficiently long time, the resonant detuning would result in a depletion of the spectrum energy, since all particles would essentially be no longer resonant.

2. Derivation of the Poisson equation for the BPS

The BPS (O'Neil & Malmberg 1968; O'Neil et al. 1971) deals with a fast electron beam injected into a one-dimensional plasma, which is treated as a cold linear dielectric medium supporting longitudinal electrostatic Langmuir waves. This section is devoted to setting up the proper form of the Poisson equation for the electric field. Such a construction follows and clarifies the general ideas contained in O'Neil et al. (1971). In particular, we explicitly outline how the shift in frequency appearing in the electric susceptivity coincides with the modulation of the electric field amplitude induced by the interaction with the beam particles.

Let us now consider the distribution function in the velocity space $F^p(v)$ of the equilibrium one-dimensional thermal plasma, the associated perturbed plasma electron distribution function $\delta f^p(t,x,v)$ and the fast electron tenuous beam distribution function $f(t,x,v)$ (here $x$ is the coordinate of the one-dimensional periodic slab of plasma). The interaction between the plasma and the beam electrons is mediated by the electric field $E(x,t)$, according to the Poisson equation

(2.1)\begin{equation} \partial_x E = - 4{\rm \pi} e \int_{-\infty}^{\infty}\,\mathrm{d}v (\,f+ \delta f^p), \end{equation}

where $e$ denotes the electron charge. The equation above determines the electric field via the charge density of the external beam and the charge displacement induced in the thermal plasma. It is worth noting that the charge distribution associated with $F^p$ is cancelled by the ion neutralizing background, whose dynamics are neglected completely in this treatment, due to the slow response of ion motion. Expanding in Fourier integral both $E$ and all the space-dependent distribution functions, the Poisson equation is written, for each wavenumber $k$, as

(2.2)\begin{equation} \mathrm{i}kE_k(t)=-4{\rm \pi} e \int_{-\infty}^{\infty}\,\mathrm{d}v(\,f_k(t,v)+\delta f_k^p(t,v)). \end{equation}

Neglecting the coupling between different (non-zero) $k$, $\delta f_k^p$ obeys the following reduced Vlasov equation:

(2.3)\begin{equation} \partial_t \delta f_k^p(t,v) + \mathrm{i}kv\delta f_k^p(t,v) = \frac{e}{m}E_k(t)\partial_v ( F^p(v) + \delta f_0^p(t,v)). \end{equation}

Here, we considered only the contribution coming from the zeroth wavenumber because such homogeneous ($x$-independent) contribution is the only one containing the equilibrium distribution function derivative in the particle velocity. All of the other components $\delta f_k^p$ would provide an intrinsically nonlinear contribution to the Vlasov equation. If we additionally consider $|\delta f_0^p|\ll F^p$, this equation describes the linear response of the plasma as a dielectric structure. The solution of (2.3) reads

(2.4)\begin{equation} \delta f_k^p(t,v) = \frac{e}{m} \int_0^t \,\mathrm{d}t^{\prime} E_k(t^{\prime}) \exp({\mathrm{i}kv(t^{\prime} - t)})\partial_v (F^p(v)+\delta f_0^p(t',v)). \end{equation}

Without loss of generality, we can use the following functional form for the electric field:

(2.5)\begin{equation} E_k(t)=\bar{E}_k \exp \left[-\mathrm{i}\int_0^t \,\mathrm{d}t^{\prime}\omega_k(t^{\prime}) \right], \end{equation}

where $\bar {E}_k=\textrm {const.}$ As a natural consequence, we can also write

(2.6)\begin{equation} E_k(t^{\prime})=E_k(t)\exp \left[-\mathrm{i}\int_t^{t^{\prime}}\,\mathrm{d}y\,\omega_k(y)\right]. \end{equation}

Using the general expression of the electric field above, (2.4) can be rewritten as

(2.7)\begin{equation} \delta f_k^p(v,t) = \frac{e}{m} E_k(t) \int_0^t\,\mathrm{d}t^{\prime} \exp \left[\mathrm{i}kv(t^{\prime}-t)-\mathrm{i}\int_t^{t^{\prime}}\,\mathrm{d}y\, \omega_k(y)\right] \partial_v (F^p(v)+\delta f_0^p(t',v)). \end{equation}

In this expression, in order to focus on resonant-driven Langmuir modes, we can set $kv\simeq \omega _k\simeq \omega _p$ (O'Neil & Malmberg 1968), where $\omega _p=\sqrt {4{\rm \pi} n_p e^2/m}$ denotes the plasma frequency ($m$ is the electron mass). Hence we explicitly set $\omega _k=\omega _p+\delta \omega _k$, where $\delta \omega _k$ is a small complex function. We observe that the integral in $t^{\prime }$ contains only one rapidly oscillating term $\exp ({-\mathrm {i}kv(t^{\prime }-t)})\simeq \exp ({-\mathrm {i}\omega _p(t^{\prime }-t)})$ multiplied by slow varying terms: $F^p$ does not depend on time and, for a cold plasma, it is peaked around $v\simeq 0$, while $\delta f_0^p$ is expected to have a typical time scale much larger than $1/\omega _p$. Thus, we can bring out of the integral in $t^{\prime }$ the term in the distribution function derivative. Moreover, the time integral in $t^{\prime }$ is dominated by the values taken in its extremes (where the rapidly oscillating term does not phase mix), especially at $t=t^{\prime }$ where one has no damping due to the imaginary part of $\delta \omega _k$: there, the integral in $y$ gives almost $\omega _k(t)(t^{\prime }-t)$ (see the appendix in O'Neil et al. (1971)). We finally get

(2.8)\begin{equation} \delta f_k^p(v,t) = \frac{e}{m}E_k(t) \frac{\partial_v (F^p(v)+\delta f_0^p(t,v))}{\mathrm{i}(kv-\omega_k(t))}. \end{equation}

Substituting this expression in the Poisson equation (2.2), we obtain

(2.9)\begin{align} &\mathrm{i}kE_k(t)\left[1-\frac{4{\rm \pi} e^2}{mk}\int_{-\infty}^{\infty}\,\mathrm{d}v \frac{\partial_v F^p(v)}{kv-\omega_k(t)} \right]\nonumber\\ &\qquad\equiv \mathrm{i}kE_k(t)\epsilon(\omega_p+\delta \omega _k(t), k) \nonumber\\ &\qquad=-4{\rm \pi} e\left[ \int_{-\infty}^{\infty}\,\mathrm{d}v f_k(t,v)\right] - \mathrm{i}kE_k(t) \delta \epsilon (t,k), \end{align}

where $\epsilon (\omega _p+\delta \omega _k(t), k)$ is the instantaneous (linear) dielectric function, while the correction $\delta \epsilon (t,k)$ is defined as

(2.10)\begin{equation} \delta \epsilon \equiv - \frac{4{\rm \pi} e^2}{mk} \int_{-\infty}^{\infty}\,\mathrm{d}v \frac{\partial_v\delta f_0^p(t,v)}{kv - \omega_k(t)}. \end{equation}

In the case of a cold plasma, i.e. $\epsilon =1-(\omega _p/\omega )^2$, for a Langmuir wave ($\omega _k\simeq \omega _p$), we immediately get the relation

(2.11)\begin{align} & \mathrm{i}k\epsilon E_k(t) \simeq \mathrm{i}k\frac{2}{\omega_p}\delta \omega _k(t) E_k(t) \nonumber\\ &\qquad = -\frac{2k}{\omega _p}\partial_t E_k(t) - 2\mathrm{i}kE_k(t), \end{align}

where we made use of the representation (2.5) for the electric field. Finally, if we neglect the small nonlinear contribution $\delta \varepsilon$, the Poisson equation takes the well-known form

(2.12)\begin{equation} \partial_t E_k(t) = - \mathrm{i} \omega _p E_k(t) + \frac{2{\rm \pi} e\omega _p}{k}\int_{-\infty}^{+\infty}\,\mathrm{d}v f_k(t,v). \end{equation}

In this derivation, we have outlined the relation existing between the function $\delta \omega _k$ appearing in (2.5) and the effective frequency shift of the dielectric function.

3. Hamiltonian description of the beam–plasma interaction

The BPS is here addressed by considering the basic assumption that the beam density $n_B$ is much smaller than that of the background plasma $n_p$. We adopt the Hamiltonian $N$-body formulation of the problem described in Carlevaro et al. (2015) and Carlevaro et al. (2016a) and references therein, where the broad supra-thermal particle beam self-consistently evolves in the presence of $M$ modes at the plasma frequency, i.e. $\omega _j\simeq \omega _p$ for $j=1, \ldots , M$. This ensures that the dielectric function of the cold background plasma, i.e. $\epsilon =1-\omega _p^2/\omega _j^2$, is nearly vanishing (O'Neil & Malmberg 1968) and allows casting the Poisson equation for each plasma oscillation into the form of a simple evolution equation. In this scheme, particle trajectories are solved from Newton's law (O'Neil et al. 1971), while one single mode in the fluctuation spectrum is initially set in order to be resonantly excited (linearly unstable) with a wavenumber $k=\omega /v_r$, where $v_r$ is a given initial resonance velocity.

As already discussed, the one-dimensional cold plasma equilibrium is taken as a periodic slab of length $L$, and the position along the $x$ direction for each particle is now labelled by $x_i$ while $N$ denotes the total particle number ($i=1, \ldots , N$). Beam particle positions are scaled as $\bar {x}_i=x_i(2{\rm \pi} /L)$, and the Langmuir wave scalar potential $\varphi (x,t)$ is expressed in terms of the Fourier components $\varphi _j(k_j,t)$ (where $k_j$ is the wavenumber). Introducing the parameter $\eta =n_B/n_p$, we use the following dimensionless variables: $\tau =t\omega _p$, $u_i=\bar {x}_i'=v_i(2{\rm \pi} /L)/\omega _p$, $\ell _j=k_j(2{\rm \pi} /L)^{-1}$ (integers), $\phi _j=(2{\rm \pi} /L)^2 e\varphi _j/m\omega _p^2$ and $\bar {\phi }_j=\phi _j\,\mathrm {e}^{-\mathrm {i}\tau }$. The prime denotes $\tau$ derivative and barred frequencies and growth rates associated with the beam–plasma instability are normalized as $\bar {\omega }=\omega /\omega _p$ and $\bar {\gamma }=\gamma /\omega _p$, respectively.

The BPS is governed by the following system:

(3.1a)\begin{gather} \bar{x}_i' =u_i \;, \end{gather}

(3.1b)\begin{gather}u_i' =\sum_{j=1}^{M} \left(\mathrm{i}\ell_j\bar{\phi}_j\exp({\mathrm{i}\ell_j\bar{x}_{i}})+\mathrm{c.c.}\right), \end{gather}

(3.1c)\begin{gather}\bar{\phi}_j' =-\mathrm{i}\bar{\phi}_j+\frac{\mathrm{i}\eta}{2\ell_j^2 N}\sum_{i=1}^{N} \exp({-\mathrm{i}\ell_j\bar{x}_{i}}). \end{gather}

For a single mode, we remark that resonance conditions are rewritten $\ell u_{r}=\bar {\omega }$. Moreover we assume that the warm beam is initialized in the velocity space with an assigned distribution function $F(u)$ (the time-dependent beam profile is denoted by $f_B(\tau ,u)=f_0(\tau ,u)/n_B$). In the simulations, we analyse a total of $N=10^{6}$ particles and we solve (3.1) using a Runge–Kutta (fourth-order) algorithm. The initialization in the velocity space is formal and, thus, particles are free to spread in this coordinate, while we set uniform initial conditions for particle positions in $[0, 2{\rm \pi}$] (we also apply periodic boundary conditions for $\bar {x}_i$ in this range). More specifically, we discretize the adopted initial beam profile in 500 delta-like beams each of them having a single fixed velocity. Each cold beam is initialized also with random generated positions and with a number of particles provided by the discretization of the considered initial profile (normalized to $N$). This provides the initial $2\times N$ array for $(\bar {x}_i,u_i)$, which is dynamically evolved in time self-consistently with the modes. Finally, the modes are initialized at $\mathcal {O}(10^{-12})$, in order to guarantee the initial linear stage of evolution. For the considered time scales and for an integration step $\Delta \tau =0.1$, as we discuss in detail in § 5, both the total energy and momentum are safely conserved.

3.1. Relevance of the BPS to the tokamak configuration

It is well known that the energetic particle (EP) excitations of global instabilities in burning plasmas and the corresponding nonlinear transport are multi-scale processes (Chen & Zonca 2016) and, thus, reduced models for the profile relaxation are often considered in view of computational savings. When addressing multiple Alfvén eigenmodes (AEs), the quasi-linear model (Berk et al. 1995, 1996; Ghantous et al. 2012) corresponds to the typical reduced description (a detailed comparison between reduced and fully nonlinear schemes for turbulence-driven EP transport can be found in Bourdelle et al. (2016)). As already stated, in the group of reduced models, the BPS can be considered to face the analysis of the nonlinear interactions between EPs and Alfvénic fluctuations (Berk & Breizman 1990a,b,c; Chen & Zonca 2013; Zonca et al. 2015). In particular, a comprehensive study of the applications and specific restrictions is given in Breizman & Sharapov (2011) and Chen & Zonca (2016).

In the following, with the aim of illustrating the direct implementation of the BPS as a paradigm for describing the burning plasma scenario, we introduce a mapping technique between the reduced radial profile of real systems and the velocity coordinate of the BPS. A validation of this approach can be found in Carlevaro et al. (2019) (see also Carlevaro et al. 2016b), where the ITER 15 MA baseline scenario is addressed in the presence of the least damped 27 toroidal AEs, and also in Biancalani et al. (2018) for an application to the case of EP-driven geodesic acoustic modes.

The relation between the independent variables is a linear correspondence which preserves the nonlinear EP profile redistribution. In particular, the AE/EP scheme must be dimensionally reduced (using standard averaged distribution functions) in order to describe the dynamics as a one-dimensional non-autonomous system (Briguglio et al. 2014). The map can be defined for a fixed single resonance, starting from the corresponding resonance condition:

(3.2)\begin{equation} \bar{\omega}^{ \rm AE}(s)-\bar{\omega}^{ \rm AE}(s_r)=\varOmega\, \ell(u-u_r), \end{equation}

where $\varOmega$ is a normalization constant to be determined during the procedure, and we use standard nomenclature for the EP/AE scenario.¹ This mapping procedure is defined as local, by using the expansion $\bar {\omega }^\textrm {AE}(s)=\bar {\omega }^\textrm {AE}(s_r)+(s-s_r)\partial _s\bar {\omega }^\textrm {AE}$, which is always valid when the frequency results in a sufficiently smooth function of the radius over the resonance region. The procedure can be finalized by implementing the suitable boundary conditions $s=0 \mapsto u=u_{\max }$ and $s=1 \mapsto u=0$, and by introducing a parameter $\ell _1$ which substantially fixes the arbitrary periodicity length $L$ for the BPS one-dimensional slab (cf. the normalization introduced above) and it corresponds to the minimum mode number that can be analysed. The closed linear link between the radial profile and the velocity coordinate reads now

(3.3)\begin{equation} u=(1-s)/\ell_1. \end{equation}

In this scheme, multiple modes can be introduced in the dynamics simply by considering the corresponding resonance conditions written in the radial space: $\ell _{j}=\ell _1/(1-s_{rj})$. A real predictive comparison of the two systems can be implemented once the density parameter of the BPS, i.e. $\eta$, is determined. Without entering the details of this step, $\eta$ can be evaluated by assuming a direct proportionality between the linear growth rates (normalized to the mode frequency) of the BPS and the AE/EP scenario. The proportionality constant can be fixed requiring that the nonlinear transport is preserved (Biancalani et al. 2018; Carlevaro et al. 2019). Then, the integration of the linear dispersion relation for the BPS, which is analysed in the next section, directly provides the value of the density parameter.

This analysis underlines how the BPS can be addressed as a paradigm for describing specific burning plasma scenarios. The powerful features of this approach can be seen in the linear link between $s$ and $u$, and in the very reduced computational demand of the BPS. In this respect, in the following sections, we analyse some specific aspects of the beam–plasma instability having in mind their direct outcomes when implemented in the radial profile of real systems.

4. Linear dispersion relation: non-perturbative effects

In this section, we review and go into more depth of some specific features of the dispersion relation. The following analysis is valid on a purely linear regime, having the aim of clarifying the non-local character of the distribution function in the velocity space.

The dispersion relation, in correspondence to electric field perturbations $\propto \exp ({\mathrm {i}(kx-\omega t)})$, is formally written (O'Neil & Malmberg 1968; Lifshitz & Pitaevskii 1976)

(4.1)\begin{equation} 1-\frac{\omega_p^2}{\omega^2}=\frac{\eta \omega_p^2}{k^2} \int_{-\infty}^{+\infty}\,\mathrm{d}v\,\frac{k\partial_v \hat{F}(v)}{k v-\omega}, \end{equation}

where the left-hand side corresponds to the cold plasma dielectric function $\epsilon (\omega )$ and $\hat {F}(v)$ is the normalized (to unity) initial beam distribution function in the velocity space.

Since the dielectric is expanded near $\omega \simeq \omega _p$ in (3.1) (O'Neil et al. 1971; Carlevaro et al. 2015), the left-hand side of the dispersion relation must be consistently rewritten as $\epsilon \simeq 2(\bar {\omega }-1)$. In order to study the instability features of one single resonant mode having mode number $\ell$, we explicitly write $\bar {\omega }=\bar {\omega }_0+\mathrm {i}\bar {\gamma }_L$, where $\bar {\omega }_0$ contains a small real frequency shift with respect to the Langmuir mode frequency $\omega _p$ (the resonance velocity is, thus, $u_r=\bar {\omega }_0/\ell \simeq 1/\ell$) and $\bar {\gamma }_L$ denotes the normalized linear growth rate. Using the proper normalization, we can finally rewrite the dispersion relation as

(4.2)\begin{equation} 2(\bar{\omega}_0+\mathrm{i}\bar{\gamma}_L-1)-\frac{\eta}{\ell M} \int_{-\infty}^{+\infty}\,\mathrm{d}u\,\frac{\partial_u F(u)}{\ell u-\bar{\omega}_0-\mathrm{i}\bar{\gamma}_L}=0, \end{equation}

where we transformed $\hat {F}(v)\to F(u)$ with $M=\int \,\mathrm {d}u F(u)$.

Equation (4.2) is numerically integrated for different mode number by fixing the parameter $\eta$ and, as initial distribution function, we consider a half-Gaussian positive slope $F(u)=0.5\, \text {Erfc}[a-b\,u]$ (with $a\simeq 6.8$ and $b\simeq 4537$).The corresponding resonant velocities occur in different regions of the distribution function varying from $u_r\simeq 0.0011$ to $u_r\simeq 0.0019$. From the numerical integration plotted in figure 1 and specified for a reference case $\eta =0.0007$, it emerges how growth rate values increase as the associated resonant velocity approaches the inflection point of $F(u)$, where, of course, the drive of the inverse Landau damping $\partial _u F$ is maximum. Furthermore, the real part of the frequency $\bar {\omega }_0$ starts to deviate from unity on the right-hand flat profile of the distribution function. Actually, in that region the local value of $\partial _u F$ is essentially negligible and a finite growth rate (necessarily associated with a real frequency shift with respect to the plasma frequency) marks the breakdown of the perturbative Landau damping expression.

Figure 1. Plot of the linear growth rate $\bar {\gamma }_L$ $(a)$ and of the real part of the frequency $\bar {\omega }_0$ $(b)$ as a function of the mode-number, obtained numerically by integrating (4.2) for the reference case $\eta =0.0007$. Green dashed line qualitatively depicts (off-scale) the initial distribution function represented assuming $F(u=1/\ell )$.

Let us now compare the obtained results with respect to a standard analytical and a semi-analytical integration of the dispersion relation. By linearizing (4.2) fixing $\bar {\omega }_0=1$, one can get the well-known Landau formula for the linear growth rate (here denoted by the subscript lin) (Lifshitz & Pitaevskii 1976):

(4.3)\begin{equation} \bar{\gamma}_{\textrm{lin}}=\left. \frac{{\rm \pi}\eta}{2\ell^2 M}\partial_u F\right|_{u=1/\ell}. \end{equation}

Meanwhile, the dispersion relation can be also analytically integrated in terms of the Faddeeva function $w(x)$ and residue contributions (Fried & Conte 1961). We easily get $\partial _u F=A\exp [-B^2(C-Du)^{2}]$, and (4.2) provides the following expression (specified using the subscript $F$) (Idouakass et al. 2016):

(4.4)\begin{equation} 2\left(\bar{\omega}_{0F}+\mathrm{i}\bar{\gamma}_F-1\right)-\frac{\eta A}{\ell M}\,\mathrm{i}{\rm \pi}(2 \mathrm{e}^{-X}-w(X))=0, \end{equation}

with $X=BC-BD(\bar {\omega }_{0F}+\mathrm {i}\bar {\gamma }_F)/\ell$ (for the sake of completeness, in our specific case, we have: $A=2559.88$, $B=5.89$, $C=1.15$ and $D=770.0$), whose roots can be numerically determined.

In figure 2, we plot the ratios between the growth rates $\bar {\gamma }$ fitted from the time behaviour of single-mode simulations of (3.1) and the obtained values of $\bar {\gamma }_L$, $\bar {\gamma }_{\textrm {lin}}$ and $\bar {\gamma }_F$ from (4.2), (4.3) and (4.4), respectively, as a function of the analysed mode number for $\eta =0.0007$. While the perturbative inverse Landau damping expression can deviate up to $100\,\%$, the predicted $\bar {\gamma }_F$ values from (4.4) appear more precise; nonetheless an underestimation of the simulated growth rate up to $20\,\%$ emerges in correspondence with the associated $u_r$ larger than the inflection point. The $\bar {\gamma }_L$ values calculated from (4.2) almost perfectly match the numerical simulations.

Figure 2. Ratios between the growth rate $\bar {\gamma }$ (fitted from simulations) and $\bar {\gamma }_L$, $\bar {\gamma }_{\textrm {lin}}$ and $\bar {\gamma }_F$ evaluated from (4.2), (4.3) and (4.4), respectively, for different mode numbers ($\eta =0.0007$).

It is worth nothing that the considered linear problem actually contains a subtle issue concerning the non-perturbative character of the dispersion relation. The integral over the whole distribution function and the nonlinear dependence upon the frequency make the estimate of the growth rate and the real frequency shift with respect to the plasma frequency intrinsically non-local. In particular, focusing on a region close to the inflection point of the distribution, the two methods giving $\bar {\gamma }_L$ and $\bar {\gamma }_F$, respectively, almost coincide, while the perturbative inverse Landau damping growth rate $\bar {\gamma }_{\textrm {lin}}$ overestimates the simulation results by approximately $10\,\%$. This is due to the fact that the residue and the principal value contribution are of the same order in this region.

We conclude by observing that, when addressing the Landau growth rate expression (4.3), we evaluate the quantity $\partial _u F$ assuming the resonance at $u=1/\ell$, i.e. we impose $\bar {\omega }_0=1$ ($\omega = \omega _p$) without considering the real frequency shift due to the presence of the beam. Actually, when the beam is tenuous (as in the considered case), we are near the marginal stability corresponding to $\bar {\gamma }_L\ll 1$ and $\bar {\omega }_0\simeq 1$. Thus, neglecting the frequency shift when evaluating (4.3) is expected to be a good approximation up to the considered order in $\bar {\gamma }$. This is certainly true when we consider modes resonant with velocities for which $\partial _u F$ is near a maximum value, i.e. near the inflection point of the distribution function. In fact, if we calculate the correction due to the real frequency shift, it is of order $\partial _u^2 F|_{u=1/\ell } (\bar {\omega }_0-1)$, being clearly of higher order in the addressed velocity region. This allows us to consider the assumption $\bar {\omega }_0=1$ as a closed and simple prescription to describe the inverse Landau damping associated with the presence of a tenuous beam.

However, since we are investigating the dispersion relation in the flat regions of the distribution function (where $\partial _u F$ is very small but $\partial _u^2 F$ is finite), it becomes important to evaluate the modification of the Landau growth rate when the real frequency shift is taken into account as directly evaluated from the exact dispersion relation (4.2). In figure 3, we plot the ratio between $\bar {\gamma }_{\textrm {lin}}$ calculated using $\bar {\omega }_0=1$ and that (dubbed $\bar {\gamma }_{\textrm {lin}}^*$) determined using the numerical frequency shift, i.e. by considering $\partial _u F|_{u=\bar {\omega }_0/\ell }$ in (4.3). It is immediately evident that this ratio is $\lesssim 1$ and the real frequency shift contribution is, thus, unable to eliminate the discrepancy with respect to the fitted growth rates from the simulation. This fact is significant because it guarantees that such observed discrepancy is really due to the non-local contributions coming to the whole profile of the distribution function, which can never be mimicked by a locally calculated growth rate. It is just the relevant effect of the non-local nature of the dispersion relation, the main issue of the present analysis, deepening the pioneering studies in O'Neil & Malmberg (1968); see also Fried & Gould (1961) for a study of the dispersion relation in a hot plasma (in this latter work, the existence was outlined of a strongly damped electronic acoustic mode).

Figure 3. Ratio between $\bar {\gamma }_{\textrm {lin}}$ using $\bar {\omega }_0=1$ and $\bar {\gamma }_{\textrm {lin}}^*$ evaluated considering $\partial _u F|_{u=\bar {\omega }_0/\ell }$ in (4.3), where $\bar {\omega }_0$ are plotted in figure 1 ($\eta =0.0007$).

5. Return current in the beam–plasma instability

As already discussed, in the original treatment the plasma is considered as a linear dielectric. The validity of this assumption relies on the condition $\eta ^{1/3}\ll 1$, inasmuch as nonlinear terms of the expansion of the background charge density are shown to give higher-order contributions to the Poisson equation (see appendix in O'Neil & Malmberg 1968).

In this section, we provide a derivation of the BPS dynamics based on the description of the background plasma as a fully ionized fluid. In the magnetohydrodynamic scheme, its density and velocity are governed by the continuity and Euler equations, respectively. In the latter, we neglect any dissipative effects, such as friction, which we will separately examine in the next section. According to the assumption of a cold plasma, we also neglect the pressure gradient by considering only the electrostatic interaction expressed by the Poisson equation.

Although we consider only linear perturbations to the background plasma, as in the linear dielectric approach, the fluid representation allows us to give a quantitative estimate of the emerging non-zero current in the plasma, consisting of a combination of two effects: a fast oscillation of the electrons at the plasma frequency $\omega _p$, having a small amplitude and essentially non-relevant impact on the long-time-scale dynamics; and a return current in the opposite direction to that associated with the beam dynamics. The return current evolution closely resembles that of the electric field: after an initial exponential growing phase, its amplitude saturates and starts oscillating around a non-zero value at the bounce frequency of the plasma electrons in the potential well. Its inclusion affects the growth rate of the electric field instability, even more markedly in the multi-mode case, where also the saturation levels are systematically lower.

In the following, $\bar {n}_e$ and $\bar {n}_B$ denote the plasma electron and the beam densities, respectively, and $u_e$ characterizes the plasma electron fluid velocity. Barred densities are normalized in units of the uniform density $n_p$ of the neutralizing ion background, assumed still.

In this setting, the equations describing the BPS take the form

(5.1a)\begin{gather} \bar{n}_e' + \partial_{\bar{x}}\left(\bar{n}_e\,u_e\right) = 0, \end{gather}

(5.1b)\begin{gather}u_e' + u_e\partial_{\bar{x}}u_e - \partial_{\bar{x}}\phi=0 , \end{gather}

(5.1c)\begin{gather}-\partial^2_{\bar{x}}\phi-1+\bar{n}_e+\bar{n}_B=0. \end{gather}

We now look for small perturbations in the plasma electron density and velocity, indicated by subscript ‘1’, i.e. $\bar {n}_e=1+\bar {n}_1$ and $u_e=0+u_1$. Substituting these expressions in the system above, and retaining only first-order terms, we find the following relation between $\bar {n}_1$ and $\bar {n}_B$:

(5.2)\begin{equation} (\partial^2_{\tau} + 1)\,\bar{n}_1=-\bar{n}_B, \end{equation}

whose formal solution is given by

(5.3)\begin{equation} \bar{n}_1=\tilde{n}_{1} - \int_0^{\tau}\,\mathrm{d}\tau' (\sin(\tau-\tau') \,\bar{n}_B(\bar{x},\tau')), \end{equation}

where $\tilde {n}_{1}=A\cos (\tau +\zeta )$ is the solution of the associated homogeneous equation, corresponding to fast Langmuir oscillations at frequency $\omega _p$. In the numerical simulations, we will assume the amplitude $A$ to be spatially constant and random phase $\zeta \in [0,2{\rm \pi} )$. The integral term in (5.3) can be rewritten as

(5.4)\begin{equation} -\bar{n}_B(\bar{x},\tau)+\eta\cos(\tau)+\int_0^{\tau} \,\mathrm{d}\tau' ( \cos(\tau-\tau') \,\partial_{\tau^{\prime}}\bar{n}_B(\bar{x},\tau')). \end{equation}

Let us analyse the physical meaning of the expression above. The first term is an instantaneous neutralization of the net negative charge introduced by the beam in the system. The second term corresponds to a fast oscillation at the plasma frequency, with amplitude directly proportional to the density ratio $\eta$, which is a small parameter in our analysis. Finally, the integral term contains the non-trivial interaction between the two densities (its behaviour will be characterized using numerical simulations at the end of the section).

Considering now (5.1a), still retaining only linear terms, we find the following expression for the velocity perturbation:

(5.5)\begin{equation} u_1= \int\,\mathrm{d}\bar{x} \left[ A\sin(\tau+\zeta)+\int_0^{\tau}\,\mathrm{d}\tau^{\prime}(\cos(\tau-\tau^{\prime})\,\bar{n}_B(\bar{x},\tau^{\prime}))\right], \end{equation}

which also provides, with opposite sign, the first-order plasma current defined as $\bar {J}_1=J_1/e=-\bar {n}_e u_e = -u_1$. As stressed above, also this velocity is a combination of a rapidly oscillating term, coming from Langmuir waves and giving no net effect, and an integral term akin to that encountered in (5.4). The latter gives a net contribution that can be identified as a return current which is directed in the opposite direction to the stream of the beam electrons.

Let us now perform a numerical study implementing the same approach of O'Neil et al. (1971), as in § 3. More specifically, we consider the beam density as a sum of particles at coordinates $\bar {x}_i$. By Fourier transforming the space variable, we obtain

(5.6)\begin{equation} \bar{n}_B(\bar{x},\tau)=\frac{\eta}{N}\sum_{j=0}^{M}\sum_{i=1}^N \exp(\mathrm{i}\ell_j\bar{x}_i), \end{equation}

where we recall that $N$ is the number of particles and $M$ the number of Fourier modes. The final system reads

(5.7a)\begin{gather} \bar{x}_i'=u_i, \end{gather}

(5.7b)\begin{gather}u_i'=\sum_{j=1}^{M}(\mathrm{i}\,\ell_j\,\phi_j\exp({\mathrm{i}\ell_j\bar{x}_{i}})+\textrm{c.c.}), \end{gather}

(5.7c)\begin{gather}\phi_j =\frac{1}{\ell_j^2}\left( \tilde{n}_{1j} + \frac{\eta}{N}\sum_{i=1}^N \exp({-\mathrm{i}\ell_j\bar{x}_i}) + \frac{\mathrm{i}\eta}{2N}(g^{+}_{j}-g^{-}_{j}) \right), \end{gather}

(5.7d)\begin{gather}\bar{J}_{1j} = -\frac{\mathrm{i}}{\ell_j} \left[ A\sin(\tau+\zeta)-\frac{\eta}{2N}(g^+_{j}+g^-_{j})\right], \end{gather}

(5.7e)\begin{gather}\left(g^{\pm}_{j}\right)' = \pm \mathrm{i} g^{\pm}_{j} + \frac{\eta}{N}\sum_{i=1}^N \exp({-\mathrm{i}\ell_j\bar{x}_i}), \end{gather}

where we underline that the final expression of the Poisson equation, (5.7c), is algebraic, and we have introduced the two auxiliary variables $g^{\pm }$ simply to obtain a fully differential system.

As a first step, we study the single-mode case setting $\ell _1=666$ (associated with a resonant velocity $u_r\simeq 1.5\times 10^{-3}$) as the only mode included in the dynamics and by initializing the particles with the positive slope $F(u)$ introduced in the previous section (also the density parameter $\eta$ is fixed as 0.0007). The temporal evolution with and without back-reaction (i.e. (3.1)) of the electric potential is shown in figure 4, while the current perturbation is plotted in figure 5.

Figure 4. Single-mode simulation: temporal evolution of the electric potential modulus in the case without $(a)$ and with $(b)$ back-reaction in logarithmic scale. Dashed lines correspond to the numerical fit for the linear phase.

Figure 5. Single-mode simulation: temporal evolution of the first-order current perturbation (equation (5.7d)) induced in the plasma ($\eta =0.0007$).

As it emerges from figure 4, after a longer pre-linear stage, the standard linear regime is clearly recovered and followed by mode saturation (nonlinear phase) and oscillation. The value of the field at saturation results in being slightly lower if the back-reaction is on (for the analysed set-up, it is lower by a factor of about $0.97$), indicating how part of the total energy is spent on the excitation of the plasma velocity perturbation. More specifically, the total energy of the system in the presence of back-reaction can be written as follows:

(5.8)\begin{equation} E=\frac{1}{2N}\sum_{i=1}^Nu_i^2+\sum_{j=1}^M\frac{\ell_j^2}{\eta}\left|\phi_j\right|^2+\sum_{j=1}^M\frac{1}{\eta}\left|\bar{J}_{1j}\right|^2 = \textrm{const.}, \end{equation}

where we stress that the first two terms of this expression correspond to the total energy of the system without back-reaction (O'Neil et al. 1971). In that case, the energy could only be exchanged between beam electrons and the electric potential, while now the third term proportional to $\bar {J}_1$ enters the dynamics. This corresponds to the energy of the excited plasma back-reaction current, whose incidence on the total energy conservation is shown in figure 6.

Figure 6. Single-mode simulation: evolution of the total energy $E$ from (5.8), normalized to its initial value $E_0$, with (solid line) and without (dashed line) the contribution of the background plasma current.

To further analyse the energy exchange between the plasma and beam electrons, in figure 7 we plot the fitted linear growth rates in the single-mode dynamics as a function of the density parameter $\eta$. Of course, in the case with back-reaction, the obtained growth rates are always below the standard ones leading to relative errors (for the analysed cases) of the order of 15–20 %. Since the linear growth rate is proportional to the nonlinear velocity spread of beam particles (Carlevaro, Montani & Zonca 2018), the outlined behaviour underlines the importance of the return current for the transport properties of the system and the resulting relaxed beam profiles.

Figure 7. Single-mode simulation: growth rate with and without back-reaction (specified in the plot) as a function of the density parameter $\eta$.

Let us now study the case of a broad beam initialized with a full Gaussian $F(u)=\exp (-(u-\bar {u})^2/2\sigma )/\sqrt {2{\rm \pi} \sigma }$, with $\bar {u}=1/600$ and $\sigma =5.67\times 10^{-5}$. We address $M=21$ modes chosen to resonate in the interval $0.00156\leqslant u \leqslant 0.00178$, i.e. $\ell _{\min }=560$ and $\ell _{\max }=640$, having equispaced resonant velocities ($\eta =0.0007$). With the chosen set-up, the particle velocity spreads associated with the mode saturation allow the beam particles to move from one resonance to the other. Each mode, due to the limited width of the velocity spread, overlaps with its neighbours allowing a cascade process of diffusive kind, i.e. the so-called quasi-linear dynamics (cf. Carlevaro et al. (2016a) and references therein for a Hamiltonian treatment of this process, in which the back-reaction is neglected). The evolution of the electric potentials is shown in figure 8, where we also report the case without back-reaction as reference.The obtained results confirm the previous analysis outlining the larger energy exchange between particles and modes than in the standard case (here, it is evident also in the saturation levels). In view of an extension of this approach towards a non-tenuous beam interacting with a plasma, it is clear that the return current would play an important physical role, not only because its intensity would be large, but also because it would be a source (in a two-dimensional configuration) of a no-longer-negligible magnetic field gradient.

Figure 8. Multi-mode simulation: mode evolution in the case with (solid lines) and without (dashed lines) back-reaction. We follow a standard rainbow colour scheme from $\ell _{\min }=560$ (red) to $\ell _{\max }=640$ (blue).

6. Source of friction and resonance detuning

Up to now, we have considered the BPS as unaffected by dissipation phenomena. In what follows, we address a dissipative effect on the beam particles, corresponding to the friction force exerted by the plasma electrons at sufficiently large impact parameter. This phenomenon has the merit of being suitably modelled in the considered scenario, and it can be explicitly evaluated assuming that a cylindrical portion of the plasma, with its axis aligned with the beam trajectory, is removed. The radius of this cylinder is taken larger than the Debye length of the plasma, so that we are excluding short impact parameter collisions in the present study. Such a simplification of the friction term clearly disregards a significant contribution due to local collisions which, however, would unavoidably yield scattering, thus violating the one-dimensionality of the BPS and requiring an extension of the model to more spatial dimensions. Furthermore, we are interested in investigating how friction affects processes like wave–particle resonance and frequency detuning more than providing a quantitative estimate of the whole friction effect in the plasma.

We observe that, when implementing the BPS in the characterization of the interaction of fast ions with Alfvénic modes, the presence of sources and sinks is often included (see Berk & Breizman 1990a,b,c). Furthermore, such interaction is characterized by a drive, due to the radial gradient of the pressure profile, and a damping rate, due to the background plasma properties. It is just this damping effect that the following analysis is able to mimic, since it also deals with the beam–plasma interaction due to the thermal plasma properties. Thus, in addition to the direct physical interest for the interaction of the beam particles with large-impact-factor thermal electrons, this scenario offers also an improved dynamical picture to be used for predicting features of fast ions in real tokamak devices.

In this respect, to take into account dissipative effects on beam particle dynamics, a new force term is added to the motion equations. In particular, in Neufeld & Ritchie (1955), it was suggested to deal with contributions having a short range and long range separately, and then to consider them in one single force term (there, heavy particles have been considered). Here, for the sake of consistency with the BPS model, we neglect collisions with short range and we focus on the collective effects of the far region of the plasma on the motion of beam particles. As stressed above, such effects are obtained by subtracting a region of the plasma contained in a cylinder of radius $\rho _{\min }$ around the trajectory of the electron beam, and then by calculating the potential induced by the injected particles in the far region of plasma, wherein polarization is produced. Then, the field generated by the polarized plasma region on the trajectories is the cause of the friction force undergone by each particle. In this way, in accordance with Neufeld & Ritchie (1955), the force term (specified for the beam electrons) can be expressed as follows:

(6.1)\begin{equation} F_{\textrm{fric}}=\frac{e^{2}\rho_{\min}}{\omega_p {x'}^{3}}K_{0} \left( \frac{\rho_{\min}}{{x'}} \right) K_{1} \left( \frac{\rho_{\min}}{{x'}} \right), \end{equation}

where $K_{0}$ and $K_{1}$ are the second kind modified Bessel functions of order 0 and $1$, respectively. The critical parameter $\rho _{\min }$ has to obey the following two restrictions:

(6.2)\begin{equation} p_{\bot} = \frac{2 e^{2}}{\rho_{\min}\omega_p {x'}} \ll m \langle v_e^{2}\rangle^{1/2}, \end{equation}

ensuring that the motion of incident beam particles slightly perturbs the plasma unperturbed electrons (having mean square velocity $\langle v_e^{2}\rangle$) at a certain distance $\rho _{\min }$ (here $p_{\bot }$ denotes the momentum gain of the plasma electrons by the passage of the incident particle), and

(6.3)\begin{equation} \rho_{\min} \gg \lambda_{D}, \end{equation}

where $\lambda _D$ denotes the Debye length of the plasma.

Considering a cold plasma, i.e. $\omega _p {x'}\gg \langle v_e^{2}\rangle ^{1/2}$, we can assume as reference case $\omega _p {x'}=100\langle v_e^{2}\rangle ^{1/2}$ (thus, the minimum value of $\rho _{\min }$ is of order $10\lambda _{D}$). Finally, the set of equations in (3.1) can be generalized to include the considered friction effect as follows:

(6.4a)\begin{gather} \bar{x}_i' =u_i , \end{gather}

(6.4b)\begin{gather}u'_{i} = \sum_{j=1}^{M}\left( \mathrm{i}\ell_{j} {\bar{\phi_{j}}} \exp\left({\mathrm{i}\ell_{j}\bar{x}_{i}}\right) + \mathrm{c.c.}\right) - A \frac{\delta^{4}}{u_{i}^{3} N_{D}} K_{0}\left(B \frac{\delta}{u_{i}}\right) K_{1}\left(B \frac{\delta}{u_{i}}\right), \end{gather}

(6.4c)\begin{gather}\bar{\phi}'_{j} = -\mathrm{i} \bar{\phi}_{j} + \frac{\mathrm{i} \eta }{2 \ell_{j}^{2} N} \sum_{i=1}^{N} \exp({-\mathrm{i}\ell_{j}\bar{x}_{i} }), \end{gather}

where $A=(2 {\rm \pi} )^{4}/4 {\rm \pi} \times 10$, $B=10 \times 2 {\rm \pi}$ , $\delta =\lambda _{D}/L$ and $N_D$ is the Debye number of the plasma. The parameter $\delta$ weights the relevance of the friction term with respect to the electric force. The range of the values of such a parameter is fixed by the upper bound $\delta \leqslant 10^{-5}$, a consequence of the condition $u \gg \langle v_e^{2}\rangle ^{1/2}$, and by a lower bound due the computational demand of the simulation.

In figure 9, we plot the spectral evolution for the case of 21 self-interacting modes, with and without the friction term, assuming an initial Gaussian distribution function in the velocity space (as represented in figure 11). Other parameters characterizing the simulation are $\eta =7 \times 10^{-6}$, $N_{D}=10^{3}$ and again $\ell _{\min }\leqslant \ell _j\leqslant \ell _{\max }$, having equispaced resonant velocities. It is worth stressing that several modes are suppressed at large time scales, but this is not a general feature since some modes remain unaffected. In order to better understand this feature, we also analyse the behaviour of the spectrum, plotted in figure 10. First of all, it is notable that, at late times, while the ideal BPS shows a steady state of the spectrum, the considered dissipative model does not reach this limit. In fact, Coulomb friction produces a progressive loss of energy in the beam particles, with kinetic energy transferred to the thermal plasma in the form of a very weak re-heating effect (we recall that the beam is here a tenuous one). As a result, the resonance between particles and field is therefore shifted towards lower velocity values, corresponding to higher wavenumbers in the field spectrum (here, the frequency of the plasma is assigned as a fixed parameter) and consequently, we observe the deformation of the field spectrum amplitude, as described in figure 10.

Figure 9. Warm beam, multi-mode simulation: mode evolution without dissipation, (3.1), and in the case of a friction coefficient $\delta =10^{-5}$ in (6.4), as indicated in the plots. Colour scheme as in figure 8.

Figure 10. Spectral evolution at different times, as indicated above the plots, with and without friction (same as figure 9). We note that, as in the previous section, we fixed $\ell _{\min }=560$ and $\ell _{\max }=640$, resonating in the interval $0.00156\leqslant u \leqslant 0.00178$ (see also figure 11).

Figure 11. Plots of the temporal evolution of the velocity distribution function at different instants, as indicated above the graphs, with and without friction (same as figure 9). The green dashed line denotes the initial beam velocity profile.

It is easy to realize that the basic mechanism entering the spectral evolution is a progressive detuning of each resonance due to the particle loss of energy. Clearly, the modes which are resonant with high-velocity particles are the first to be affected by this mechanism along the system evolution. Differently, those modes which exchange power with low-velocity particles are still supported by that amount of them which was interacting with the previous set of modes. The validity of this conjecture is clearly confirmed by figure 11, where the distribution function behaviour is shown and the dragging of particles in the velocity space is outlined. In fact, the distribution function evolution without friction (upper panels) is characterized by the typical flattening, while, in the presence of the friction term, a new morphology emerges in which a bump corresponding to small velocities is present. If the number of available modes at large wavenumbers were infinite, we would see the particle distribution function progressively flatten at low velocity and the kinetic energy of the beam unavoidably dissipated. In this sense, a real steady state here cannot exist, but the described effect must be considered as a very slow phenomenon, so that, on the time scale on which the $N$-body approach remains valid (see the discussion in Carlevaro et al. (2014) on the so-called quasi-stationary states), its effect appears negligible, while still playing a significant role in the resonance phenomenon.

This study is relevant when the BPS model is implemented to mimic the radial transport in a tokamak device. In fact, if EPs lose energy for specific phenomena taking place in the tokamak configuration, this can affect the AE spectrum. For instance, in Schneller, Lauber & Briguglio (2016), it is shown that the low toroidal number modes can be destabilized via a cascade process. These modes would be, using the mapping procedure addressed in § 3.1, those ones surviving in figure 10. Thus, including friction term in the BPS model can be important in reproducing specific tokamak features due, for instance, to collisional effects.

7. Concluding remarks

We analysed different aspects of the beam–plasma instability clarifying which features of the original treatment have to be improved to adapt this scheme to the interaction between fast ions and the Alfvénic spectrum in a tokamak device.

As a first step, we gave a reformulation of the Poisson equation for the BPS to outline the basic hypotheses at the ground of the dielectric assumption which, in the original analysis, was made for the background plasma.

Then, we studied the so-called non-perturbative effects of the beam–plasma dispersion relation, showing how the inverse Landau damping formula fails approaching a flat region of the fast particle distribution function. Actually, in such a condition, a significant frequency shift is naturally predicted by the distribution function and a non-vanishing growth rate emerges from the integral nature of the dispersion relation.

After that, we analysed the determination of a return current into the bulk of the plasma as an effect of the fast particle crossing. In this respect, we treated the interaction of the beam with the background plasma fluctuations, describing the former still as an $N$-body Hamiltonian system, while the latter in terms of linear fluid dynamics. We then arrived at a restated $N$-body dynamics, whose Poisson equation is now able to account for the interaction of fast particles with the Langmuir waves present in the plasma. We showed how the emergence of a return current implies a suppression of the saturation level of the spectrum and the obtained mode growth rates. This fact suggests that, in the case of a non-tenuous beam, the dielectric approximation, if still adopted, must be amended for effects of the plasma back-reaction that, in turn, influence the beam dynamics and the resulting beam profile.

Finally, we studied the correction to the BPS when the interaction between fast particles and plasma particles far from the beam trajectory is taken into account. A friction term arises due to the plasma charges at scales much greater than the Debye length. This effect determines a significant displacement of the beam particles in the distribution function from high velocities to smaller ones. The resulting field spectrum is initially damped in correspondence to larger wavelengths, and eventually the whole profile is affected.

We stress that, in both these last analyses, the ions of the thermal plasma are considered as essentially at rest, so that they behave here simply as a ‘neutralizing background’ for the system. This hypothesis is well grounded for a sufficiently cold plasma, since the ion response takes place on a time scale much slower than both the inverse of electron plasma and bounce frequencies. This scenario is natural for the BPS and easily amendable if the physical conditions require one to account for ion motion. In other words, we are applying to plasma physics the so-called Born–Oppenheimer approximation typical of the physics of matter, when molecular spectra are concerned. In plasma physics, such an assumption lives on a pure classical sector and it is justified by the mass separation existing between ions and electrons. However, it is important to distinguish here between the single-electron velocity, like that of the beam particles, and the electron fluid velocity, characterizing the background response via a return current. The single-electron velocity is always much larger that the single-ion one, but the electron fluid, the velocity of which is an average of the single-particle motion, remains in principle small, as happens in our study, due to the validity of quasi-neutrality in the unperturbed plasma at the physical scale of interest.

We conclude that, when implementing the results of this study into real tokamak configurations, the considered effects are not directly relevant: return current and friction are associated with the motion of fast ion beams, but their description requires a different setting with respect to the present one, especially because of the radial non-uniformity of the tokamak equilibrium.

The relevance of the considered deviations from a standard BPS picture lies in the possibility of mimicking corresponding features of the mapped scheme from the velocity space to the radial profile of a tokamak. In particular, while the drive of the instability of fast ions interacting with Alfvénic modes is essentially governed by the radial profile of the pressure of the hot particle population, i.e. the radial slope of its distribution function, the damping undergone by the ion beam depends on intrinsic features of the background equilibrium plasma. The relevance to dealing with a BPS accounting for damping features relies on modelling this effect when the mapping procedure is implemented.