Quasi-linear model for the beam-plasma instability: analysis of the self-consistent evolution

TL;DR

This paper re-examines the quasi-linear model for beam-plasma instability, comparing theoretical predictions with simulations, identifying discrepancies, and proposing corrections to improve the model's accuracy during early evolution stages.

Contribution

It introduces a correction to the quasi-linear growth rate by accounting for the distribution function's time dependence, enhancing the model's predictive capability.

Findings

Diffusive asymptotic dynamics are well predicted.

Early mesoscale transport is less accurately modeled.

A revised growth rate improves spectral evolution predictions.

Abstract

We re-analyze the quasi-linear self consistent dynamics for the beam-plasma instability, by comparing the theory predictions to numerical simulations of the corresponding Hamiltonian system. While the diffusive features of the asymptotic dynamics are reliably predicted, the early temporal mesoscale transport appears less efficient in reproducing the convective feature of the self-consistent scenario. As a result, we identify the origin of the observed discrepancy in the underlying quasi-linear model assumption that the distribution function is quasi-stationary. Furthermore, we provide a correction to the instantaneous quasi-linear growth rate based on a linear expansion of the distribution function time dependence, and we successfully test this revised formulation for the spectral evolution during the temporal mesoscale.

Equations

$$⟨E_k{E}_{k^{\prime }}⟩=|E_k{|}^2\delta (k+k^{\prime }),$$

$$⟨E(t,v)⟩=0,⟨E(t,v)E(t+\tau ,v)⟩=\Delta (v)\delta (\tau ).$$

$${\partial }_{t}P={\partial }_v(\mathcal{D}(v){\partial }_{v}P),$$

$${\partial }_t{f}_0(t,v)={\partial }_v({\mathcal{D}}_{QL}{\partial }_v{f}_0),$$

$${\partial }_t|E_k(t){|}^2=2{\gamma }_k^{QL}|E_k{|}^2,$$

$$\displaystyle\mathcal{D}_{QL}(t,v)\equiv\frac{e^{2}\pi\mathcal{N}}{m^{2}v}|E_{k}|^{2}\big{|}_{k=\omega_{p}/v}\;,$$

$$\displaystyle\gamma^{QL}_{k}(t)\equiv\frac{2\pi^{2}e^{2}\omega_{p}}{mk^{2}}\partial_{v}f_{0}\big{|}_{v=\omega_{p}/k}\;,$$

$$f_B=F_B+\bar{\mathcal{N}}M{\eta }^{-1}{\partial }_u(\mathcal{J}{u}^{-5}),$$

$${\partial }_{\tau }\mathcal{J}=\pi \eta M^{-1}{u}^2\mathcal{J}{\partial }_u{F}_B+\pi \bar{\mathcal{N}}{u}^2\mathcal{J}{\partial }_u^2(\mathcal{J}{u}^{-5}),$$

$$2({\bar{\omega }}_{0\ell }+i{\bar{\gamma }}_{\ell }-1)-\frac{\eta u_{\ell }}{M}{\int }_{-\infty }^{+\infty }du\frac{{\partial }_u{F}_B}{u/u_{\ell }-{\bar{\omega }}_{0\ell }-i{\bar{\gamma }}_{\ell }}=0,$$

$$\displaystyle\bar{\gamma}_{L\ell}=\frac{1}{2}\,\pi\eta M^{-1}u_{\ell}^{2}\partial_{u}F_{B}\big{|}_{u_{\ell}}\;,$$

$$\begin{split}&\partial_{\tau}\bar{x}_{i}=u_{i}\;,\\ &\partial_{\tau}u_{i}=\sum_{\ell}\big{(}i\,\ell\;\bar{\phi}_{\ell}\;e^{i\ell\bar{x}_{i}}+c.c.\big{)}\;,\\[-15.0pt] &\partial_{\tau}\phi_{\ell}=-i\phi_{\ell}+\frac{i\eta}{2\ell^{2}N}\sum_{i=1}^{N}e^{-i\ell\bar{x}_{i}}\;,\end{split}$$

$${\partial }_t{f}_k(t)=-ikv{f}_k+\frac{e}{m}\sum _q{E}_{k-q}{\partial }_v{f}_q,$$

$${\partial }_t{E}_k(t)=-i{\omega }_p{E}_k(t)+\frac{2\pi e{\omega }_p}{k}{\int }_{-\infty }^{+\infty }{f}_k(t,v)dv.$$

$$f_{k}(t,v)=\frac{e}{m}\int_{0}^{t}dt^{\prime}E_{k}(t^{\prime})\exp\big{[}ikv(t^{\prime}-t)\big{]}\partial_{v}f_{0}(t^{\prime},v)\;.$$

$$E_{k}(t)=\bar{E}_{k}\exp\Big{[}-i\int_{0}^{t}\omega_{k}(t^{\prime})dt^{\prime}\Big{]}\;,\qquad E_{k}(t^{\prime})=E_{k}(t)\exp\Big{[}-i\int_{t}^{t^{\prime}}\omega_{k}(y)dy\Big{]}\;,$$

$$\displaystyle\partial_{t}f_{0}(t,v)=\frac{e^{2}}{m^{2}}\sum_{k>0}|E_{k}(t)|^{2}\partial_{v}\Big{(}$$

$$\displaystyle\partial_{t}|E_{k}|^{2}=\frac{2\pi e^{2}\omega_{p}}{mk}|E_{k}|^{2}\Big{(}\int_{-\infty}^{\infty}dv\int_{0}^{t}dt^{\prime}\exp\Big{[}ikv(t^{\prime}-t)-i\int^{t^{\prime}}_{t}dy\,\omega_{k}(y)\Big{]}\partial_{v}f_{0}(t^{\prime},v)+c.c.\Big{)}\;.$$

$$\int_{0}^{t}dt^{\prime}\frac{1}{i(kv-\omega_{k}(t^{\prime}))}\partial_{v}f_{0}(t^{\prime},v)\partial_{t^{\prime}}\Big{(}\exp\Big{[}ikv(t^{\prime}-t)-i\int_{t}^{t^{\prime}}\omega_{k}(y)dy\Big{]}\Big{)}\;.$$

$$\displaystyle\int_{0}^{t}dt^{\prime}\frac{1}{i(kv-\omega_{k})}\partial_{v}f_{0}(t,v)\partial_{t^{\prime}}\Big{(}\exp\Big{[}ikv(t^{\prime}-t)-i\int_{t}^{t^{\prime}}\omega_{k}(u)du\Big{]}\Big{)}+\qquad\quad$$

$$\displaystyle-\int_{0}^{t}dt^{\prime}\frac{1}{i(kv-\omega_{k})}\partial_{t}\partial_{v}f_{0}(t,v)\exp\Big{[}ikv(t^{\prime}-t)-i\int_{t}^{t^{\prime}}\omega_{k}(u)du\Big{]}\,,$$

$$\frac{1}{i(kv-{\omega }_k(t))}{\partial }_v{f}_0(t,v)+\frac{1}{(kv-{\omega }_k(t){)}^2}{\partial }_t{\partial }_v{f}_0(t,v).$$

$${\partial }_t{f}_0={\partial }_v({\mathcal{D}}_{QL}(t,v){\partial }_v{f}_0+{\mathcal{D}}_1(t,v){\partial }_t{\partial }_v{f}_0),$$

$$\displaystyle\mathcal{D}_{QL}(t,v)\equiv\frac{e^{2}}{m^{2}}\sum_{k>0}|E_{k}(t)|^{2}\Big{[}\frac{1}{i(kv-\omega_{k})}-\frac{1}{i(kv-\omega^{*}_{k})}\Big{]}\;.$$

$$\displaystyle\mathcal{D}_{1}(t,v)\equiv\frac{e^{2}}{m^{2}}\sum_{k>0}\mid E_{k}(t)\mid^{2}\Big{[}\frac{1}{\left(kv-\omega_{k}\right)^{2}}+\frac{1}{\left(kv-\omega^{*}_{k}\right)^{2}}\Big{]}\;.$$

$${\partial }_t|E_k(t){|}^2=2{\Gamma }_k|E_k{|}^2,$$

$${\Gamma }_k(t)={\gamma }_k^{QL}+\delta {\gamma }_k,$$

$$\displaystyle\delta\gamma_{k}(t)=\frac{\pi e^{2}\omega_{p}}{mk}\int_{-\infty}^{+\infty}dv\Big{[}\frac{1}{\left(kv-\omega_{k}\right)^{2}}+\frac{1}{\left(kv-\omega_{k}^{*}\right)^{2}}\Big{]}\partial_{t}\partial_{v}f_{0}\;,$$

$$2{\int }_{-\infty }^{+\infty }dv\frac{(kv-{\omega }_p{)}^2-{\Gamma }_k^2}{{\left[(kv-{\omega }_p{)}^2+{\Gamma }_k^2\right]}^2}{\partial }_t{\partial }_v{f}_0=-2{\int }_{-\infty }^{+\infty }dvG(kv-{\omega }_p){\partial }_t{f}_0,$$

$$G(y)=k\,\partial_{y}\Big{[}\frac{y^{2}-\Gamma_{k}^{2}}{\left[y^{2}+\Gamma_{k}^{2}\right]^{2}}\Big{]}=2k\,y\,\frac{3\Gamma_{k}^{2}-y^{2}}{\left[y^{2}+\Gamma_{k}^{2}\right]^{3}}\;.$$

$$-2{\int }_{-\infty }^{+\infty }dvG(kv-{\omega }_p){\partial }_t{f}_0\sim -2{\partial }_t{f}_0(kv={\omega }_p-\bar{y}){\int }_{-\infty }^{{\omega }_p/k}dvG(kv-{\omega }_p)+$$

$$-2{\partial }_t{f}_0(kv={\omega }_p+\bar{y}){\int }_{{\omega }_p/k}^{+\infty }dvG(kv-{\omega }_p)\sim -\frac{4(\sqrt{2}-1)}{k{\Gamma }_k}{\partial }_v{\partial }_t{f}_0{|}_{v={\omega }_p/k}.$$

$$\delta {\gamma }_k=-\frac{4(\sqrt{2}-1)\pi e^2{\omega }_p}{m{k}^2{\Gamma }_k}{\partial }_t{\partial }_v{f}_0{|}_{v={\omega }_p/k},$$

$${\Gamma }_k={\gamma }_k^{QL}-\frac{4(\sqrt{2}-1)\pi e^2{\omega }_p}{m{k}^2{\Gamma }_k}{\partial }_t{\partial }_v{f}_0{|}_{v={\omega }_p/k}.$$

$$\bar{\Gamma }=\frac{{\bar{\gamma }}^{QL}}{2}+\frac{{\bar{\gamma }}^{QL}}{2}\sqrt{1-4{\bar{\delta }}^2/({\bar{\gamma }}^{QL}{)}^2},\text{with}{\bar{\gamma }}^{QL}(u)=\frac{\pi \eta u^2}{2M}{\partial }_u{f}_B,{\bar{\delta }}^2(u)=\frac{(\sqrt{2}-1)\eta u^2}{M}{\partial }_{\tau }{\partial }_u{f}_B,$$

Full text

I Introduction

One of the most useful paradigm of the theoretical analysis of Plasma Physics is the so-call bump-on-tail (BoT) scheme, describing the instability induced in a thermal plasma due to a suprathermal population of charged particles ZCrmp . The relevance the BoT paradigm relies on its capability to capture the main features of the fast particle radial transport in Tokamak devices, due to a one-to-one correspondence between the gradients in the velocity distribution finction and the pressure gradients in the corresponding radial direction BB90a ; BB90b ; BB90c . As shown in Refs.EEbook ; AEE98 ; EE08 and in Refs.CFMZJPP ; ncentropy , the BoT scenario can be suitably described via a beam-plasma system (BPS) model by a cold background dielectric on which fast electrons interact with the longitudinal Langmuir waves, as first postulated in Refs.OM68 ; OWM71 , see also Refs.MK78 ; TMM94 ; KV14 and Ref.L72 for the case of a warm beam.

We will use a $N$-body Hamiltonian BPS formulation for testing the predictivity of the quasi-liner (QL) dynamics. When the Langmuir spectrum with which the tenuous fast electron beam interacts is sufficiently continuous and large (in comparison to the non-linear velocity spread), it is well-known, since from the original treatment Pines (see also the detailed monographs Refs.MT64 ; KT73 ; D82 ), that the transport in the velocity space acquires a marked diffusive character and the particle distribution function obey a Fokker-Planck equation, whose diffusion coefficient is proportional to the electric field spectrum intensity. This behavior accounts for a random walk of the electron velocity, due to the continuous scatter with the superposition of many different wave components, having a random phase distribution pesme94 . Since the beam is tenuous, the electric field amplitude is to be regarded as small even after its saturation and it constitutes a small control parameter to ensure the consistency of the QL approximation, especially about the possibility to neglect cross coupling terms in the wave amplitudes. For a satisfactory discussion about the validity of the QL model and its limit of application, see Ref.Laval99 and refs. therein. While an implementation of the QL diffusive scenario to describe radial transport in Tokamak devices is given in Ref.GG12 , where the so-called $1.5$ model is addressed. Morever, in Ref.GBG14 it is provided a re-analysis of the Line Broadened QL model BBFW95 and a comparison with respect to a fully Vlasov solver for the one dimensional reduced scheme in the presence of collisional effects and focusing on the case of isolated modes.

Clearly, the pure diffusive QL model properly accounts for the BoT evolution only in sufficiently late phases when each particle is, on average, scattering with all the field components simultaneously and no significant trapping processes take place. As it has been outlined in Ref.ncentropy , such a situation is not appropriate to describe the temporal mesoscale evolution, where convective transport can be important. In fact, in this intermediate region of the evolution, a significant fraction of particles can coherently interact with a finite number of wave components, giving rise to ballistic motion in the velocity space and avalanche phenomena. For the relevance of such features in the radial transport of fast ions in ITER relevant simulations, see Ref.spb16 . In Ref.ncentropy , it has been emphasized that trace of such coherent evolution of the particle-field interaction can be found even in the QL dynamics, as far as the self-consistent evolution is accounted to some extent.

In this paper, we focus our attention just to the emergence of such non-diffusive transport in the early non-linear phases of the BPS, in order to evaluate the predictivity of the QL paradigm in the temporal mesoscale, when structures in the phase space (like clumps and holes) can still survive. To this end, after a qualitative introduction of the problem, we write down the QL model equations, in particular the single partial differential equation (PDE) for the spectrum evolution. and, in order to characterize the prediction of the QL theory, we compare the profiles of the calculated distribution function with those ones obtained by the Hamiltonian $N$-body BPS CFMZJPP ; ncentropy . In this work, we adopt such a description which directly follows from the original paper Ref.OWM71 . It is worth mention that also particle-in-cell (PIC) numerical simulations of the BoT instability is widely addressed in plasma physics (see, for example, the seminal papers Refs.DS68 ; MN69 ; D90 ). As well known, PIC codes are based on a Lagrangian frame and sample the plasma density as a series of macro-particles, which represent parcels of density. Particles themselves are assigned a weighting, which represents the amount of phase-space density represented by the particle, and a position in phase space (details of this approach can be found in Refs.H55 ; HE88 ). For more recent applications, see, e.g., Refs.D94 ; EB15 ; BDE17 and refs therein, and Ref.DF14 . The most significant discrepancy, emerging in the above discussed comparison, concerns the system time “reactivity”, *i.e., *the QL dynamics appears as slower with respect the exact Hamiltonian one (see also Ref.LP83 ). While the non-linear velocity spread is qualitatively reproduced, say with enough precision in view of the radial Tokamak transport (see the analysis in Ref.nceps16 ), the achievement of the plateau configuration and, on average, the whole mesoscale evolution appear a bit delayed in the QL scenario.

The main merit of the present analysis consists of giving a clear explanation for such a slower response of the QL model in comparison to a real BPS, in terms of the quasi-stationary assumption at the ground of the derivation of a pure diffusive dynamics. In the early non-liner stages of the BPS, bouncing-like processes, as discussed in Ref.OWM71 , take place and the dynamics of the distribution function is associated to the same time scale of the electric field: in such a regime, the quasi-stationary assumption must be relaxed, or at least improved, to get a more realistic picture of the system, as it appears in the Hamiltonian model. This information is of a significant impact when, like in Ref.nceps16 , we want to predict the radial transport of fast ions in Tokamak configuration via a one-to-one mapping from the velocity to the physical space. While, as shown in Ref.GG12 the prediction of the diffusive QL evolution are of interest for the radial transport, here we suggest that the mesoscale evolution can not be reliably mapped before a suitable upgrade of the QL self-consistent picture is reached.

We make a first step to upgrade the QL model by considering a more refined approximation scheme for the quasi-stationary distribution function, in which a linear in time contribution is retained. As shown by using the numerical distribution function, the modified growth rate is enhanced with respect to the QL case, implying that the system reactivity increases. Hence, a better agreement with simulations, thus a better characterization of the mesoscale dynamics, is achieved.

The paper is organized as follows. In Sec.II, we provide some basic concepts introduced via a qualitative picture of the problem. Furthermore, using the fundamental hypotheses of the QL model, we arrive to the basic equation for the electric field spectral evolution, which constitutes the tool to test the self-consistency of the QL paradigm. In Sec.III, we introduce a schematic description of the Hamiltonian BPS and a discussion of the comparison with the QL prediction for the distribution function morphology is provided. In Sec.IV, the growth rate is explicitly evaluated by relaxing the quasi-stationary assumption for the distribution function, and tests for the reliability of the predicted dynamics are numerically provided. In Sec.V, concluding remarks follows.

II Quasi-linear dynamics

We focus on a warm beam of electrons, having an initial homogeneous spatial distribution with number density $n_{B}$, interacting with a cold homogeneous plasma of density $n_{p}\gg n_{B}$. The dynamics of a beam in the sea of Langmuir waves present in the plasma is faced on a 1D level. The addressed electric field $E(x,t)$ is expressed using the Fourier decomposition $E(t,x)=\sum E_{k}(t)e^{ikx}$, where $k$ denote the considered wave numbers for the Langmuir spectrum.

Those electrons which have a velocity $v$ very close to the phase velocity of an electric wave component $v_{Ph}=\omega_{p}/k$ ($\omega_{p}$ being the plasma frequency, $\omega_{p}^{2}=4\pi n_{p}e^{2}/m$, where $e$ and $m$ denote the electron charge and mass, respectively) can efficiently exchange power with such a mode, resulting trapped in the local electric well. In this sense, let us introduce the trapping (bounce) frequency $\omega_{B}$ defined by $\omega_{B}^{2}\equiv 2ek|E_{k}^{S}|/m$, where $|E_{k}^{S}|$ denotes the saturated values of $E_{k}$. In fact, after a linear growth due to the inverse Landau damping process oneil65 , the electric field reaches a saturated value corresponding to the maximum energy extraction from the resonant electrons.

The time scale $\tau_{B}=2\pi/\omega_{B}$ characterizes the typical trapping time of an electron within the potential. We did not label $\omega_{B}$ by the suffix $k$, because, in what follows, we will intend to deal with an average value, estimated with an average electric field saturation level or with the most relevant Fourier component. Another important time scale is the autocorrelation time $\tau_{ac}$, *i.e., *the typical time an electron takes to cross the spatial size $\Delta x$ of the wave packet, namely $\tau_{ac}=\Delta x/v=1/(v\Delta k)$, $\Delta k$ being the packet spread in wavenumbers. Also in this estimate, the velocity must be intended as the typical resonant one of the beam. Clearly, the trapping of a large number of electrons in the electric field, associated to the formation of structures in the phase space, can take place only if $\tau_{B}\ll\tau_{ac}$. In such a situation, resonant particles can get trapped by a small set of modes before fell the whole fluctuation spectrum. Otherwise, for $\tau_{B}\gg\tau_{ac}$, the particles feel simultaneously all the field components and a diffusion-like process takes place.

The random character of the Langmuir electrostatic field leads us to require that its mean value, taken over a sufficiently large time, larger than $\tau_{ac}$, essentially vanishes, *i.e., *$\langle E(t,x)\rangle\simeq 0$. Furthermore, because of the average spatial homogeneity of the considered system, the variance of the electric field $\langle E^{2}\rangle$ must depend on time only: the weak turbulence arising in the system does not produce structures on spatial scales comparable to the fundamental one $L$ (here $L$ deotes the system periodicity spatial size, *i.e., *$k=2\pi\ell/L$, $\ell$ being an integer number). If we deal with a Langmuir continuum ($L$ is sufficiently large), we can evaluate the variance obtaining the basic condition

[TABLE]

where $\delta$ denotes the Dirac function. This issue clearly suggests that, as far as macroscopic homogeneity is preserved, no significant coupling can take place, on average, between Fourier components associated to different wavenumbers. Such a consideration is at the ground of the QL model.

Furthermore, in the limit $\tau_{B}\gg\tau_{ac}$, the superposition of different harmonics interacting with the same electron allows us to model the electric field as a white noise. Thus, the force equation of the system can be written in the form of a random walk in the velocity space, namely described by the Langevin equation $dv/dt=-eE(t,v)/m$, where the dependence on the particle velocity is inferred by means of the resonance condition, which selects certain groups of harmonics, while the electric field satisfies the condition

[TABLE]

where $\Delta(v)$ is a given function of the velocity. By other words, when the electron interacts with a large number of Fourier harmonics having random distributed phases, we deal with a white noise in the electron velocity space. It is well-known that the probability $d\mu$ of finding an electron in the range ($v$, $v+dv$) is $d\mu=P(t,v)dv$, where the probability density $P$ satisfies a Fokker-Planck equation of the form

[TABLE]

where $\mathcal{D}=e^{2}\Delta/2m^{2}$. This picture closely reproduces some hypotheses and the main predictions of the QL model in the limit of sufficiently large times with respect to $\tau_{ac}$. However, we stress how the condition (1) (spatial homogeneity of the field) is much general than the more restrictive assumption (2) (Gaussian field) responsible for a diffusive behavior. In this paper, we are interested to study the self-consistent dynamics allowed by the preservation of the system homogeneity, in order to characterize the evolution from the initial linear instability phase, up to the diffusive regime.

II.1 QL equations

The QL model is derived from the Vlasov-Poisson coupled system specified for the the interaction of the suprathermal electron beam with a cold background plasma and it is expressed, in Fourier space, in terms of $E_{k}(t)$ and the distribution function components $f_{k}(t,v)$ of the beam particles. The attention is focused on behavior of $f_{k=0}\equiv f_{0}$, which is the only contribution having a non-zero initial condition since the BPS is initially spatially homogeneous. The dynamics of $f_{0}$ receives contributions from all $k$-values and, thus, it accounts for the main transport properties of the system. Moreover, according to the homogeneity condition Eq.(1), the component $f_{k}$ is assumed to receive mainly contribution only by the correspondent harmonics. Introducing the assumption that the function $f_{0}$ is not fast varying in time (quasi-stationarity hypothesis associated to the plateau formation), the QL model equations can be written as follows Pines ; LP81 :

[TABLE]

where $\mathcal{D}_{QL}(t,v)$ and $\gamma^{QL}_{k}(t)$ are the QL diffusion coefficient and instantaneous growth rate, respectively, and read

[TABLE]

here we have used the resonance condition $kv=\omega_{p}$ and since the model requires a sufficiently large spectral density, we have approximated the discrete $k$-space by a continuum. In fact, we have introduced $\mathcal{N}(k)$ as the state density and it can be assumed constant $\mathcal{N}\simeq K/\Delta k$, where $K$ is the number of the considered modes while $\Delta k=k_{max}-k_{min}$ indicates the spectral width (the interval where the spectral intensity is significantly different from zero). Out of $\Delta k$, both the spectral density and intensity vanish.

The electric field and the distribution function dynamics can be combined (see also Ref.ncentropy ) to get a partial differential equation in the spectral intensity only. First of all, by means of the resonance condition and the consideration above, we move to the velocity space introducing the electric potential modes as $E_{k}(t)\to-i\omega_{p}\varphi(t,v)/v$, where $\varphi(t,v)$ has to be intended as the continuous mode spectrum derived from the discrete one, specified for $v=\omega_{p}/k$. In order to compare the QL model with respect to numerical simulations, let us now introduce the following normalization. The 1D cold plasma equilibrium is taken as a periodic slab of length $L$ while time is taken in $\omega_{p}$ units $\tau=t\omega_{p}$. We also introduce the beam to plasma density ratio as $\eta=n_{B}/n_{p}$. Other quantities are defined as as: $\ell=k(2\pi/L)^{-1}$ (integer numbers), $u=v(2\pi/L)/\omega_{p}$, $\phi=(2\pi/L)^{2}e\varphi/m\omega_{p}^{2}$. In this scheme, we analyze the spectral intensity $\mathcal{J}(\tau,u)=|\phi|^{2}$ and we use the beam distribution function $f_{B}$ defined by transforming $f_{0}(t,v)/n_{B}\to f_{B}(\tau,u)$. The QL equations write now as

[TABLE]

where $F_{B}(u)$ indicates the initial beam profile in the $u$ space, $\bar{\mathcal{N}}=m/\Delta\ell$ (here $\Delta\ell=\ell_{max}-\ell_{min}$) and $M=\int_{-\infty}^{+\infty}duF_{B}$.

III Comparison with numerical simulations

We now compare the QL distribution function and spectrum evolutions with those ones emerging from numerical simulations of an $N$-body Hamiltonian system describing the BPS CFMZJPP ; ncentropy .

From the linear theory, the dispersion relation for a single mode having corresponding (normalized) resonant velocity $u_{\ell}=1/\ell$ reads as OM68 ; LP81

[TABLE]

(for an analysis on the accuracy of this equation, see Ref.nceps18 ). Here, the dielectric function is expanded near $\omega_{\ell}\simeq\omega_{p}$ (according to motion equations, see below) and we used $\bar{\omega}_{\ell}=\bar{\omega}_{0\ell}+i\bar{\gamma}_{\ell}$, where $\bar{\omega}_{0\ell}$ denotes the real part of the frequency and $\bar{\gamma}_{\ell}$ ($\ll 1$) the growth rate of the considered mode (the bar indicates that these quantities are in $\omega_{p}$ units). The linearized form of Eq.(10) (by assuming $\bar{\omega}_{0\ell}=1$) reduces to the well know expression

[TABLE]

which of course corresponds to the steady form Eq.(7) and, by construction, it is the discretization of the spectral growth rate of Eq.(9) (first term of the right-hand side). It is important to notice that a break-down of the perturbative Landau damping expression (Eq.(11)) is observed (compared with respect to numerical integration of Eq.(10)) for specific setup of the beam plasma system, enlightening non-perturbative and non-local character of the dispersion relation, especially in the flat region of the initial distribution function.

In what follows, we careful choose the initial conditions of the system so that the numerically evolved modes have growth rate values well predicted by the Landau expression. This way, we ensure that discrepancies in the system evolution between the numerical simulation and QL model do not rely on different drive values. In fact, the QL theory naturally defines an instantaneous Landau growth rate as far as the field amplitude is small enough.

III.1 N-body beam plasma system description

We adopt the Hamiltonian formulation of the BPS CFMZJPP ; ncentropy , where the broad fast particle beam self-consistently evolves in the presence of several modes at the plasma frequency. This ensures a nearly vanishing dielectric function of the cold background plasma ($\epsilon=1-\omega_{p}^{2}/\omega^{2}$) and allows casting the Poisson equation for each mode into the form of an evolution equation (while particle trajectories are solved from the force equation). Each beam particle position along the $x$ direction is labeled by $x_{i}$ ($N$ indicates the total particle number) and we use the normalization introduced in the previous Section (positions are scaled as $\bar{x}_{i}=x_{i}(2\pi/L)$). The BPS is governed by the following set of equations:

[TABLE]

and we recall that resonance conditions rewrite as $\ell u_{\ell}=1$. In the non-linear simulations, we initialize the beam particle numbers distributed according to a given $F_{B}(u)$ (for a total $N=10^{6}$ particles). We solve Eqs.(12) using a Runge-Kutta (4th order) algorithm and the initial conditions for the particle positions $\bar{x}_{i}$ are given uniformly in $[0,\,2\pi]$, while modes are initialized at amplitudes of order $\mathcal{O}(10^{-14})$ to ensure the starting linear regime. The resulting evolutive velocity distribution function is simply evaluated by a smoothed $x$-dimension averaged procedure of the correspondent phase-space.

III.2 Numerical results

We assume an initial Gaussian warm beam distribution function in velocity space $F_{B}(u)=\text{Exp}[a(b+u)^{2}]$ (with $a\simeq-2.2\times 10^{-6}$ and $b=-0.003$). Having set $\eta=0.0002$, we select and simulate modes for which the relative error between $\bar{\gamma}$ (numerical integration of Eq.(10)) and $\bar{\gamma}_{L}$ (Eq.(11)) is less than $5\%$ (see left-hand panel of FIG.1, where we plot the corresponding values and we also represent (off-scale) the intial distribution function to underline the mode positions). In the non-linear simulations of Eqs.(12), we run 45 modes with equispaced resonant velocities and $\ell_{min}=342$ and $\ell_{max}=489$, thus getting $\bar{\mathcal{N}}\simeq 0.3$.

The mode temporal evolution is plotted in the right-hand side of FIG.1, where the linear and non-linear phase of the dynamics can be easily recognized. We underline how the mode behavior depicted in FIG.1 is consistent with the emergent QL approximation in the sense described in Sec.II (see also Ref.ncentropy ). In fact, the present setup corresponds to an averaged estimate $\tau_{ac}/\tau_{B}\simeq 0.02$.

Let us now compare the evolution of the distribution function as derived from simulations (left-hand panel of FIG.2) with respect to the numerical integration of Eqs.(8)-(9) (right-hand panel of FIG.2). The single PDE for the spectral evolution (Eq.(9)) is integrated using a (linearized) Crank-Nicolson algorithm. The initial condition for $\mathcal{J}(0,u)$ is set Gaussian with features (amplitude, mean and variance) chosen to mimic the discrete simulated mode spectrum.

The comparison of the two behaviors clearly shows a time delay in the evolution of the QL model with respect to $N$-body simulations. In fact, from FIG.2 we see that the formation of the plateau (green curve) is slower in the QL model with respect to simulations. This result is consistent with the assumptions under which we derived the QL equations. By other words, the ansatz that the distribution function evolves slowly enough to be considered as constant in Eq.(19) affects the whole evolution of the QL model, making it less reactive in time than the real system.

Analogously, we can analyze the behavior of the spectrum in the two schemes. In FIG.3, we over-plot (for different times) the function $\mathcal{J}$ integrated from Eq.(9) and the values of $|\phi_{\ell}|^{2}$ extracted from simulations. We see again that the QL spectrum evolves more slowly since $\tau\simeq 1000$ and it is just this behavior which is responsible for the discrepancy in the distribution function according to Eq.(8). Nonetheless, at sufficiently late times, the QL spectrum satisfactory envelops the simulation results 111We would like to underline that a certain degree of scattering in the numerical spectrum determined via the $N$-body simulations is not depending on the number of simulated particles, but it is associated to the limit our mode number to reproduce the real continuum of the Langmuir spectrum. However, the addressed modes are adequate to reproduce all the feature of a QL profile and the resulting spectral morphology does not prevent a satisfactory comparison with the PDE numerical plots.. It is worth noting that Eq.(9) naturally contains the linear regime of the instability. In fact, when the field intensity is very small (according to the level of the Langmuir fluctuations), we naturally recover the linear growth of the intensity according to the standard inverse Landau damping rate. However, as soon as the non-linear contribution grows enough, it accounts for a diffusive-like behavior in the velocity space, which has no counterpart in the real simulations where such effects are absent.

Since, as discussed above, we set up our initial conditions in order to avoid non-perturbative effects in the dispersion relation, we can claim that, in the initial linear phases, the numerical and QL spectrum have to coincide, as clearly emerges from the first panel of FIG.3 ($\tau=1000$). The crucial point is that the diffusive contributions, naturally associated with the Fokker-Planck equation (4), induces a too early termination of the linear regime with respect to the exact simulations. As a consequence, the spectrum remains at a slower level with respect to the real one and the whole temporal mesoscale evolution turns out to be delayed. These considerations suggest the necessity of a revised approach able to account for the coherent particle-field interactions in the mesoscale temporal evolution of the BPS.

Finally, it is important to stress that the different width of the plateau region in the two approaches (in particular, cf. the red curves of FIG.2) has to be mainly attribute to the fact that the QL spectrum is larger than the numerical one according to the presence of a continuum of linearly stable modes, against the fixed character of the mode array used in the simulations.

IV Corrections to the instantaneous QL growth rate

The comparison with simulations outlines how the system reactivity at mesoscales is underestimated in the QL model. We are now going to show how a more reactive system is predicted by relaxing the quasi-stationary assumption for the distribution function. In Fourier space, the Vlasov-Poisson coupled system writes as

[TABLE]

We now implement the assumption that the $k$ component of the distribution function receives mainly contribution just by the correspondent harmonics. Setting $q=0$ in Eq.(13), we get the solution

[TABLE]

Without loss of generality, we can set the electric field as

[TABLE]

where $\bar{E}_{k}=const.$, $\omega_{k}=\omega_{p}+\delta\omega_{k}$, while $\delta\omega_{k}$ is a small complex function. Using the last expression in Eq.(15), the dynamics of $f_{0}$ assumes the following Dyson-like form

[TABLE]

where the reality of the electric field allowed to use the relations $E_{-k}=E^{*}_{k}$ and $\delta\omega_{-k}=-\delta\omega^{*}_{k}$. For a derivation of an equivalent Dyson-like equation in the $\omega$-space, see Refs.AK66 ; ZCrmp .

Similarly, Eq.(14) can be manipulated leading to the following evolutive equation for the electric harmonics intensity

[TABLE]

The system of Eqs.(17) and (18) constitutes the basic framework to self-consistently evolve the beam-plasma interaction under the hypothesis of average homogeneity only.

If in the region where the distribution function $f_{0}$ remains significantly different from zero (namely $v_{min}<v<v_{max}$), we assume that the Langmuir spectrum is sufficiently dense in order to have $kv\simeq\omega_{p}$ for $\forall v$, the integral taken in $dt^{\prime}$ of the two expressions above has a small exponent (the contribution $(t^{\prime}-t)\omega_{p}$ cancels out between the two parts of the exponent). In this case, such an integral can be rewritten as

[TABLE]

If $\omega_{k}$ is not a fast varying function in time and we Taylor expand $f_{0}$, retaining only the next-to-the-leading order term, i.e. $\partial_{v}f_{0}(t^{\prime},v)=\partial_{v}f_{0}(t,v)+(t^{\prime}-t)\,\partial_{t}\partial_{v}f_{0}(t,v)$, the expression (19) rewrites

[TABLE]

and the integrals can be explicitly computed (see appendix A in Ref.OWM71 ) finding

[TABLE]

A similar Taylor expansion is performed in Ref.Klimontovich , where it is applied to the calculation of the collisional integral, from which the spectral densities of the particle numbers and of the electric field are evaluated for non-ideal stationary and homogeneous plasmas, where the fluctuations are non-stationary and the delayed contribution could be relevant.

Finally, Eq.(17) can be restated as as

[TABLE]

As already discussed, the QL diffusion coefficient $\mathcal{D}_{QL}$ in Eq.(23) is found by approximating the discrete $k$-space by a continuum, *i.e., *$\sum_{k}\left(...\right)\to\int_{k_{min}}^{k_{\max}}dk\;\mathcal{N}(k)(...)$ and assuming a constant state density $\mathcal{N}$, thus getting the formal expression in Eq.(6). It is worth stressing how, in the QL model, the assumption that the function $f_{0}$ is quasi-stationary relies on the idea that the velocity beam profile reaches a stationary regime in the resonance region.

Applying similar considerations to the time integral in Eq.(18), we get the following generalized expression of Eq.(5):

[TABLE]

where we mention that the QL growth rate $\gamma^{QL}_{k}$ (see Eq.(7)) is obtained for $\mathrm{Im}(\omega_{k})\ll\omega_{p}$. The additional contribution $\delta\gamma_{k}$ has already been written in Ref.D82 (see Eq.(49) in chapter 14), where it is called the non-resonant contribution to the growth rate, but it is not evaluated and simply neglected with respect to the non-linear wave-particle and wave-vave interactions. Its impact on the BPS is here addressed.

We can rewrite the integral in Eq.(27) via a partial integration as follows

[TABLE]

where the function $G(y)$ within the integral reads

[TABLE]

This function has absolute maximum and minimum at $y=\pm(\sqrt{2}-1)\,\Gamma_{k}=\pm\bar{y}$ and it rapidly tends to vanish for $y\gg\Gamma_{k}$. Hence, we split the integral into the sum of the two contributions for negative and positive $y$ and we approximate the function $\partial_{v}f_{0}$ by its values at $\pm\bar{y}$, so getting222Here, we approximate the variation of the function $\partial_{t}f_{0}$ by its derivative with respect to velocity, valid for $\Gamma_{k}\ll\omega_{p}$.

[TABLE]

Therefore, $\delta\gamma_{k}$ reads

[TABLE]

and the following equation can be written for $\Gamma_{k}$

[TABLE]

Using the dimensionless variables introduced for Eqs.(8) and (9), the solution of the equation above writes as

[TABLE]

where we recall that barred growth rates are expressed in $\omega_{p}$ units and it is worth noting how a small deviation from the QL growth rate is found for $4|\bar{\delta}^{2}|\ll(\bar{\gamma}^{QL})^{2}$. The predicted correction to the QL growth rate is expected to be positive approaching the plateau region, where the distribution function flattens and $\partial_{u}f_{B}$ is a decreasing function of time, giving a negative value for $\bar{\delta}$. As a consequence, the value of the modified growth rate is larger than the QL instantaneous one.

The modified expression for the linear instantaneous growth rate could allow us to restate also the equation for the distribution function into an upgraded form able to account for the mesoscale dynamics. Up to first approximation, we could infer that this non-steady correction affects the diffusion coefficient only, leaving unchanged the morphology of Eq.(4) accordingly with the considerations made in Ref.Laval99 about the proper renormalization of the theory.

It is worth stressing that our correction to the QL theory is valid as far as $(t^{\prime}-t)$ remains a small quantity in the Taylor expansion of the distribution function, and this can take place for sufficiently early times only (surely it is not true in the late phases). In this respect, we can get an analytical form of the ratio between the revised expression for the spectral intensity and the QL one, at least when, during the early mesoscale evolution, they do not differ too much, *i.e., *when $\bar{\delta}^{2}/(\bar{\gamma}^{QL})^{2}\ll 1$. From Eq.(25), we naturally obtain the relation

[TABLE]

and, in the considered limit, we can perform the following expansion of Eq.(33):

[TABLE]

Substituting this expression in Eq.(34) and taking the primitive of the exact differentiation, it naturally provides the relevant form of the ratio $\mathcal{J}/\mathcal{J}_{QL}$, i.e.,

[TABLE]

where, henceforth, $\mathcal{J}_{QL}$ denotes the QL spectrum explicitly plotted in FIG.3.

Since the derivative of the distribution function clearly decreases with respect to the Gaussian initial condition as the system evolves, it naturally comes out that the expression above is, on average, describing an enhancing of the predicted spectrum with respect to the QL one.

As a check, we use, for the early temporal mesoscales of the system (says $\tau\simeq 1600$) the numerically calculated distribution function to estimate the ratio above. By other words, we use the distribution function extracted from the $N$-body dynamics to calculate $\partial_{u}f_{B}$ at the considered times. As shown in FIG.4, the amended spectrum (orange solid line) is significantly close to the numerical one (we also plot the original QL spectrum for comparison). This good agreement between the predicted and the simulated spectra indicates the validity of our derivation for the growth rate and that the addressed temporal stages are consistent with a small deviation from the QL model, when Eq.(36) reliably holds.

In order better understand the dynamical implications of the amended dynamics, we now evaluate the ratio $\bar{\Gamma}/\bar{\gamma}^{QL}$ from Eq.(33) in the full mesoscale temporal region, *i.e., *when Eq.(36) is no longer valid.

In particular, we use the morphology of the distribution function at $\tau=2000$, *i.e., *when most of the spectrum is saturating, again by means of the numerical data. This ratio (together with the individual growth rates) is plotted in FIG.5. We see that $\bar{\Gamma}$ is significantly enhanced with respect to the instantaneous QL growth rate clarifying how, as soon as the system leaves the linear phase of the instability, our calculation provides the requested increase of the system reactivity at the ground of the temporal mesoscale inadequacy of the QL model.

The anaysis above clearly otulines how the present reformulation of the growth rate for the BPS is able to account for deviations at the mesoscale evolution. However, it is worth noting that our approach still relies, like the QL one, on the negligibility of cross-coupling terms among spectral components. Actually, new physical features could be enclosed in a revised picture by addressing such non-diagonal effects as, for instance, consireded in K69 where the so-called trapped-particle instability is investigated. This scenario relies on the presence of a large amplitude electrostatic mode which is able to trigger sidebands in the spectrum. Such kind of physics is surely relevant for the possible deviation from the QL model but, as far as the beam is sufficiently tenuous, the homogeneity of the problem is on average preserved and the diagonal hypothesis remains reliable during the dynamics.

V Conclusions

In this work, we compared the dynamical evolution of the QL model with respect to the results obtained from a $N$-body simulation, outlining a qualitative but definitely distinguishable delay of the former with respect to the latter. The nature of our analysis identifies the source of such a discrepancy in the too restrictive character of the quasi-stationary assumption for the function. Thus, we suggest that, in order to significantly upgrade of the 1D QL dynamics, it is relevant to make precise the evolution of the distribution function. In this respect, we amended the QL instantaneous growth rate by accounting for the distribution function evolution via a linear temporal expansion, and we provided the implications on the self-consistent dynamics.

Although not all the original derivations do make explicit account for the quasi-stationary behavior of the distribution function, their statistical characterization of the field and particle fluctuations are implicitly underestimating the real time response of the BPS. Thus, it clearly emerges how the QL model is really predictive only for the linear stages and for the late evolution, while the derivation must be refined for the temporal mesoscales, for instance, as in Sec.IV, taking into account corrections due to the distribution function time dependence.

We validated such corrections to the growth rate, obtained by retaining the next-to-leading order in a Taylor expansion, by using the behavior of the distribution function extracted from numerical simulations. We first evaluated an analytical expression for the ratio between the predicted and QL spectra in those early mesoscale region where they are sufficiently close: we showed how a very good agreement exists between the restated spectrum and the numerical one. We also evaluated the ratio between the modified and the QL growth rates when this approximation is no longer valid and the revised expression was obtained two time larger than the QL one. The obtained result suggests to look for a more appropriate partially differential representation of the system evolution at the temporal mesoscale. This requires the possibility to restate the Vlasov equation in a suitable predictive and differential form. The determination of such a reduced model for the Dyson-Poisson system would constitute a significant improvement in the description of BPS.

Summarizing, it is worth noting how the main goal of our analysis consists of having outlined the necessity to revise the expression of the QL instantaneous growth rate in order the model be able to better fulfill $N$-body simulations. Actually, we also traced the proper way how to get a revise formula, by refining the time evolution of the distribution function contribution to the spectrum of the Langmuir waves.

The QL model satisfactory accounts for the initial linear phase of the beam-plasma instability and for the late time diffusive evolution in the velocity space, but it significantly fails when applied to the temporal mesoscale evolution. This has relevant implications also for the Tokamak physics, when a one-to-one mapping between the velocity space of the BPS and the radial coordinate of a Tokamak equilibrium is performed (for details of the mapping, seeRef.nceps18 and Refs.BB90a ; BB90b ; BB90c for the classical bump-on-tail paradigm). By other words, the present study, when mapped into the radial profile of fast ions in a Tokamak, shows that the temporal mesoscale radial transport is not satisfactory addressed by a QL code. This means that all phenomena which could be relevant during this phase of the evolution must be described out of the QL paradigm, at least by a correction of the instantaneous growth rate, due to the radial gradients of the thermodynamic parameter of the plasma. When using a QL code, the instantaneous growth of the spectrum is significantly underestimated in the field saturation regime and then relevant features of the fast ions transport could not be adequately reproduced.

Acknowledgements.

This work has been partly carried out within the framework of the EUROfusion Consortium [Enabling Research Projects: NAT (AWP17-ENR-MFE-MPG-01), MET (CfP-AWP19-ENR-01-ENEA-05)] and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.