Generation of controllable plasma wakefield noise in particle-in-cell simulations

TL;DR

This paper introduces a new method to generate controllable plasma wakefield noise in particle-in-cell simulations, enabling more accurate modeling of beam-plasma instabilities by reducing unrealistic initial noise levels.

Contribution

A novel technique using randomly placed charged rods is proposed to produce adjustable wakefield noise without parasitic effects in simulations.

Findings

Method effectively controls noise amplitude in simulations.

No parasitic Cherenkov radiation observed.

Improves accuracy of beam-plasma instability modeling.

Abstract

Numerical simulations of beam-plasma instabilities may produce quantitatively incorrect results because of unrealistically high initial noise from which the instabilities develop. Of particular importance is the wakefield noise, the potential perturbations that have a phase velocity which is equal to the beam velocity. Controlling the noise level in simulations may offer the possibility of extrapolating simulation results to the more realistic low-noise case. We propose a novel method for generating wakefield noise with a controllable amplitude by randomly located charged rods propagating ahead of the beam. We also illustrate the method with particle-in-cell simulations. The generation of this noise is not accompanied by parasitic Cherenkov radiation waves.

Equations

$$\lambda (\xi )=k_{p}Q/\pi ,-{\lambda }_p/2<\xi -{\xi }_0<0,$$

$$\lambda(\xi)=\frac{k_{p}Q}{\pi}\Bigl{(}1-\cos\bigl{(}2k_{p}(\xi-\xi_{0})\bigr{)}\Bigr{)},\quad-\lambda_{p}/2<\xi-\xi_{0}<0.$$

$$F(k_z)=\left|{\int }_{-\infty }^{\infty }\lambda (\xi )e^{i{k}_z\xi }d\xi \right|$$

$$\rho ({\vec{r}}_{\perp },\xi )=\delta ({\vec{r}}_{\perp })\lambda (\xi ),{\vec{r}}_{\perp }=(x,y),$$

$$E_{z}(\vec{r}_{\perp},\xi)=2{k_{p}}^{2}\int_{\xi}^{\infty}d\xi^{\prime}\int d\vec{r}_{\perp}{\!\!\!}^{\prime}\;\rho(\vec{r}_{\perp}{\!\!\!}^{\prime}\>,\xi^{\prime})\\ \times K_{0}\left(k_{p}|\vec{r}_{\perp}-\vec{r}_{\perp}{\!\!\!}^{\prime}\>|\right)\cos\bigl{(}k_{p}(\xi-\xi^{\prime})\bigr{)}\\ =\frac{2k_{p}^{2}Q}{\pi}K_{0}(k_{p}r)G(\xi-\xi_{0}),$$

$$G(\xi)=k_{p}\int_{\max(-\pi k_{p}^{-1},\xi)}^{0}d\xi^{\prime}\bigl{(}1-\cos(2k_{p}\xi^{\prime})\bigr{)}\cos\bigl{(}k_{p}(\xi-\xi^{\prime})\bigr{)}\\ =\begin{cases}-\frac{8}{3}\sin(\xi),&k_{p}\xi<-\pi,\\ \frac{2}{3}\bigl{(}\sin(2\xi)-2\sin(\xi)\bigr{)},&-\pi<k_{p}\xi<0,\\ 0,&\xi>0.\end{cases}$$

$$E_{\text{rms}}^2=N⟨{E_z}^2⟩,$$

$$N=2{\pi }^2{k}_p^{-1}{R}^{2}n$$

$$E_{\text{rms}}^{2}(\xi)=N\left\langle\left(\frac{2k_{p}^{2}Q}{\pi}K_{0}(k_{p}r)G(\xi-\xi_{0})\right)^{2}\right\rangle\\ =2n\frac{4k_{p}^{4}Q^{2}}{\pi^{2}}\int_{0}^{R}2\pi rK_{0}^{2}(k_{p}r)dr\,\int_{\xi}^{0}G^{2}(\xi-\xi_{0})\,d\xi_{0}\\ =\frac{8k_{p}^{2}Q^{2}n}{\pi}\Bigl{(}1+k_{p}^{2}R^{2}\bigl{(}K_{0}^{2}(k_{p}R)-K_{1}^{2}(k_{p}R)\bigr{)}\Bigr{)}\\ \times\int_{\xi}^{0}G^{2}(\xi-\xi_{0})\,d\xi_{0}.$$

$$E_{\text{rms}}^2(\xi )=\frac{8k_p^2{Q}^{2}n}{\pi }{\int }_{\xi }^0{G}^2(\xi -{\xi }_0)d{\xi }_0,$$

$$E_{\text{rms}}^2(-\infty )=\frac{256k_p{Q}^{2}n}{9}.$$

$${\bar{E}}_z\left(x,y,\xi \right)=\frac{1}{1000}\sum _{t=1201\Delta t}^{2200\Delta t}{E}_z\left(x,y,\xi ,t\right).$$

$$\mu (\xi )=\frac{1}{n_s}\sum _{i=1}^{n_s}{\bar{E}}_z^i\left(x_m,y_m,\xi \right),$$

$$\sigma (\xi )=\sqrt{\frac{1}{n_s-1}\sum _{i=1}^{n_s}{\left({\bar{E}}_z^i\left(x_m,y_m,\xi \right)-\mu \left(\xi \right)\right)}^2},$$

$$\begin{gathered} {\int }_0^{R}2\pi r{K}_0^2(k_{p}r)dr\to \\ {\int }_{-X/2}^{X/2}{\int }_{-Y/2}^{Y/2}dxdy\sum _{i,j}{K}_0^2(k_p{r}_{ij}), \end{gathered}$$

$$r_{ij}=\left|{\vec{r}}_{\perp }+iX{\vec{e}}_x+jY{\vec{e}}_y\right|,$$

Full text

I Introduction

Understanding a collective interaction of relativistic charged particle beams with plasmas is important for a wide variety of physical problems. Among them are space plasmasPoP17-120501 , fast ignition schemes for inertial fusionPoP14-055502 ; PoP12-057305 , turbulent plasma heating for magnetic fusionFST55(2T)-63 ; PoP13-062312 ; PoP17-083111 , positron bunch instabilities driven by an electron cloud in collidersPRST-AB7-124801 , plasma wakefield acceleration,PRL104-255003 ; PoP22-103110 ; PoP22-123107 and many others. In most cases, the interaction has the form of an instability that develops starting from some low-amplitude shot noise.

A realistic noise levelPRST-AB16-041301 is difficult to reproduce in numerical simulations, since the number of simulated quasi-particles (or macro-particles) is usually much smaller than the number of real particles in the system. Fewer quasi-particles, each carrying a larger charge, produce random fields with amplitudes which are orders of magnitude too high in comparison with the experiment. This, in turn, may result in faster growth of unstable perturbations and quantitatively wrong simulation results.

Since one-to-one simulations of typical beam-plasma systems of interest fall far beyond state-of-the-art computational capabilities, the only available choices are ordered initial distributions of particlesPRST-AB16-041301 and extrapolating high-noise simulation results to low-noise physical systems. The second approach is schematically illustrated in Fig. 1. Obviously it would benefit from having several simulation points with different noise levels. This makes it desirable to develop a simple method with the capability to produce a noise field with a controllable amplitude. For reliable extrapolation, the noise level must be controlled independently of key simulation parameters, like the grid resolution or the number of quasi-particles.

For many instabilities of relativistic beams, the noise harmonics of interest are potential (Langmuir) plasma waves with a phase velocity which is close to the speed of light $c$ and the wavevector $\vec{k}$ directed along the axis of beam propagation. We will call these waves wakefield noise, which is an analogy to regularly excited plasma wakefields, which have the same properties. At first glance, a controllable wakefield noise might be easily excited by an ensemble of randomly located point-like quasi-particles propagating ahead of the beam (Fig. 2). By changing the charge of the quasi-particles and their number, it is possible to control the noise level. However, relativistic quasi-particles usually emit numerical Cherenkov radiation, if simulated using particle-in-cell codes.JCP15-504 ; PRST-AB16-021301 ; EPJD68-177 ; NIMA-829-353 This radiation critically affects a clean study of the beam instability and must therefore be avoided.

In this paper we propose a novel method for generating wakefield noise with a controllable amplitude by randomly located charged rods (Sec. II). The noise is free from parasitic Cherenkov radiation. We also give expressions for the relations between the amplitude of the noise field and the rod parameters (Sec. III) and illustrate the method using particle-in-cell (PIC) simulations (Sec. IV). The main findings are summarized in Sec. V. Studies of particular beam instabilities fall beyond the scope of this paper.

For all simulations detailed in this paper we use the PIC code PSC.PSC The coordinates are either Cartesian $(x,y,z)$ or cylindrical $(r,\phi,z)$ with the $z$-axis being the direction of beam propagation. The co-moving coordinate $\xi=z-ct$ is used wherever convenient.

II The idea of charged rods

The problem of numerical Cherenkov radiation depends on the details of the numerical solver. In the present paper we use of the finite-difference time-domain (FDTD) scheme. The discretization of the FDTD scheme leads to a phase velocity of the propagated radiation, which is strictly smaller than the velocity of light. A quasi-particle with sufficient energy can therefore travel faster than the radiation on the grid. This leads to numerical Cherenkov radiation, as illustrated in Fig. 3(a). A point-like charge $Q$ moving with the speed of light emits short-wavelength electromagnetic radiation, the amplitude of which is much higher than the amplitude of the useful plasma wave. The radiation wavelength is of the order of the grid size and is much shorter than the plasma wavelength $\lambda_{p}=2\pi k_{p}^{-1}$. The plasma wavelength is determined by the plasma density $n_{0}$ through the plasma wavenumber $k_{p}=\sqrt{4\pi n_{0}e^{2}/(mc^{2})}$, where $e$ is the elementary charge and $m$ is the electron mass. The amplitude of the emitted Cherenkov electromagnetic waves is proportional to the amplitude of short-wavelength harmonics in the Fourier spectrum of the emitting charge. More specifically, the $z$-component of the wavevector is of importance. This already hints at the idea of how the numerical Cherenkov radiation can be reduced in comparison to the longer-wavelength plasma waves: The radiating source must be long and smooth. In other words, it should be shaped like a charged rod.

We consider two variants of rods: i) Rectangular rods with the linear charge density

[TABLE]

and ii) Cosine-shaped rods with

[TABLE]

For efficient excitation of Langmuir waves, the rod length must be $\lessapprox\lambda_{p}$. We set it to $\lambda_{p}/2$. This value does not maximize the wakefield of the rod and is chosen only for the simplification of the subsequent analytical calculations it provides.

The Fourier spectra

[TABLE]

of these rods and for the point charge $Q$ are shown in Fig. 4. For the wave number of plasma waves holds that $k_{z}\approx k_{p}$ (vertical dotted line in Fig. 4) and they are excited equally well by all charge distributions. Numerical Cherenkov waves for typical resolutions of PIC code simulations have wave numbers for which $k_{z}/k_{p}\sim 10^{2}-10^{3}$ (shaded area in Fig. 4) holds, and their excitation is suppressed by some orders of magnitude in the case of cosine-like rods. The smoothness degree of the charge distribution $\lambda(\xi)$ determines the decrease rate of the spectrum $F(k_{z})$ for high $k_{z}$. In the case of the cosine-like distribution this decrease rate is given by $F(k_{z})\propto k_{z}^{-3}$.

A “beam” of equally charged rods would produce not only the noise field, but also a regular wakefield which is excited by the total rod charge. This is undesirable. To get rid of the regular field, the sign of the charge of the different rods must be chosen randomly.

Another necessary condition for achieving a correct noise is the equivalence of all wakefield phases. The contributions of different rods must, on average, uniformly cover the whole period of the plasma wave. Using alternating rod charges the minimum length of the rod area is $\lambda_{p}/2$, as rods of opposite charge contribute to opposite half-periods of the wave.

Note that an exact nulling of the total charge of the rod ensemble, for example by choosing an equal number of positively and negatively charged rods, is erroneous. This additional constraint on the rod ensemble would break the equivalence of wakefield phases. This can be easily seen in the extreme case of two oppositely charged rods.

III Amplitude of the wakefield noise

In this section we will analytically calculate the root mean square (rms) of the longitudinal electric field excited by an ensemble of cosine-shaped rods (2). Assume that the rod heads (characterized by random coordinates $\vec{r}_{\perp 0}$ and $\xi_{0}$) are uniformly distributed in $z$-direction over an infinitely wide layer of thickness $\pi k_{p}^{-1}$. The average density in this region is given by $2n$, with $n$ being the average number density of rods of each charge sign ($Q$ or $-Q$).

We will first find an expression for the wakefield of a single rod. Given the charge density of the rod

[TABLE]

the longitudinal component of the wakefield isPAcc22-81

[TABLE]

where $K_{0}$ is the modified Bessel function of the second kind, and

[TABLE]

In Ref. PRST-AB16-041301, , the rms noise field of $N$ charges is calculated as

[TABLE]

where the angle brackets denote the averaging of the field of a single charge over possible locations of this charge with respect to the observation point. We have to modify our approach in this case, as we take the number of rods to be infinite, which makes the average field vanish. Assume the rods occupy a cylindrical area with a large radius $R$. Then

[TABLE]

and the average field

[TABLE]

In the limit $R\to\infty$ it holds that

[TABLE]

and for $\xi<-\pi k_{p}^{-1}$ (trailing the rods),

[TABLE]

IV Wakefield noise in PIC simulations

We now describe our approach to the production of wakefield noise in three-dimensional PIC simulations, how the noise looks like, and which additional actions are required to observe a close quantitative agreement between simulation results and the developed theory.

We configure the code to use widespread algorithms: a standard Boris pusherBoris with second order particles to perform particle pushing and a standard FDTD-schemeFDTD to perform field pushing. The simulation window moves with the speed of light in a patch-based manner, that is, by appending a simulation grid and quasi-particles on one end of the box and detaching the same amount of volume on the opposite side. This process is iterated continuously. The window size is given by $X\times Y\times Z$ with $X=Y=0.3\lambda_{p}$, $Z=2.02\lambda_{p}$. The boundary conditions are periodic in transverse dimensions ($x$ and $y$) and reflecting in the $z$ dimension. The simulated propagation length, moving in the $z$ direction, is $10\lambda_{p}$. Initially the plasma is cold and uniform. It is composed using 3 quasi-particles per cell for the electrons, while the ions are not simulated and treated as immobile charges. The spatial resolutions are $\Delta x=\Delta y=\Delta z=\lambda_{p}/130\approx 0.05k_{p}^{-1}$ and the time step is $\Delta t=0.99\Delta x/(c\sqrt{3})\approx 0.28\omega_{p}^{-1}$. The simulated time for the whole simulation is given by $T=10\lambda_{p}/c$, which results in a number of time steps of about $2260$.

At the start of the simulation, $100$ random positions are computed in the cuboid given by $0\leq x<X$, $0\leq y<Y$ and $1.5\lambda_{p}\leq z\leq 2\lambda_{p}$. These positions $\left(x_{i},y_{i},z_{i}\right)$ are used as head positions for $100$ rods of length $\lambda_{p}/2$. The rods are constructed as strings of quasi-particles with a distance of $\Delta z$ between each adjacent pair. The strings start at $\left(x_{i},y_{i},z_{i}-0.5\Delta z\right)$ and continue with decreasing z-coordinates. These quasi-particles have the same shape as the plasma quasi-particles, a relativistic factor of $\gamma=10^{10}$, and a charge that varies according to the cosine-like distribution (2). The total charge of a single rod (the sum of all charges in the constituting string of quasi-particles) is $|Q|\approx 1.8\times 10^{-3}en_{0}k_{p}^{-3}$, where the sign of the charge is chosen randomly. This charge is sufficiently high for a clear observation of the wakefield and sufficiently low to stay within the regime of a linear plasma response.

To enhance the wakefield noise against other noise harmonics (which are always present in PIC simulations), we average the fields over many time steps according to the formula

[TABLE]

This averaging strongly suppresses all perturbations except those propagating with the speed of light and being stationary in the co-moving frame.

The simulation is then run for a large number of samples $n_{s}$ with different random number generator seeds, which leads to different rod positions. Examples of the produced field distributions are shown in Fig. 5. Since the frequency $\omega_{p}$ and the phase velocity $c$ of the perturbation are fixed, the wavelength $\lambda_{p}$ is fixed as well, which makes the wakefield noise periodic in $\xi$. It therefore does not look like usual noise. The transverse distance for field correlations is about $k_{p}^{-1}$, which makes random field changes in transverse directions only visible in the case of wider simulation areas. The amplitude and the phase of the generated fields are quite random even for narrow simulation areas (Fig. 6), and their properties can be characterized by the mean value $\mu(\xi)$ and the standard deviation $\sigma(\xi)$:

[TABLE]

[TABLE]

where the superscript $i$ denotes the sample number, and $x_{m}$ and $y_{m}$ are the coordinates of the observation line. For the following we put $x_{m}=y_{m}=0.15\lambda_{p}$.

The resulting values for $\mu(\xi)$ and $\sigma(\xi)$ are shown in Fig. 7. The dotted line in Fig. 7 gives the value of $E_{\text{rms}}$ calculated for these rods according to the formula (10). It is substantially lower than the established average field in simulations. The difference results from the periodical boundary conditions. Each rod in the simulation domain has an infinite number of replicas, the positions of which are correlated with the rod position (Fig. 8). Consequently, for correct comparison of the theory and the simulation results, the radial averaging in (9) has to be modified:

[TABLE]

where

[TABLE]

$\vec{e}_{x}$ and $\vec{e}_{y}$ are unit vectors, and indices $i,j\in\mathbb{Z}$, where $i=j=0$ corresponds to the position of one specific rod (situated in the pink area in Fig. 8) and all other values correspond to the positions of its replicas (situated in gray boxes in Fig. 8). The dashed line in Fig. 7 shows the modified values for the average field.

V Summary

We have shown a reliable method of producing controllable noise levels in PIC simulations. This controllable wakefield noise does not suffer from numerical Cherenkov radiation. Analytical expressions for the rms amplitude of this noise have been given. These expressions include formulas for two different external boundaries, namely for periodic boundaries as well as for open or absorbing boundaries. These expressions and the method in general were illustrated using very common PIC simulation techniques, making them very comparable and applicable for the research community. A very good agreement has been found. This method can be used to perform noise level scans for a multitude of noise seeded physical processes. This paper also provides an understanding of the to-be-expected noise structure. Using these noise level scans, the correct behavior can be predicted by extrapolation of the generated data points. Further research should focus on the application of this novel method on specific beam-plasma non-linearities.

Acknowledgements.

Contribution of R.S. and K.L. to this work is supported by The Russian Science Foundation, grant No. 14-50-00080. N. M. and H. R. acknowledge the hospitality of the Arnold Sommerfeld Center for Theoretical Physics at the Ludwig Maximilians University. N. M. acknowledges the very useful advice of K. U. Bamberg during the development of the computational simulation. This work was supported by the Cluster-of-Excellence Munich Centre for Advanced Photonics (MAP) and the Gauss Centre for Supercomputing (GCS), project PLASMA SIMULATION CODE, LRZ-ID pr84me.