Variational Formulation of E & M Particle Simulation Algorithms in Cylindrical Geometry using an Angular Modal Decomposition

# Variational Formulation of E & M Particle Simulation Algorithms in Cylindrical Geometry using an Angular Modal Decomposition

## Abstract

Taking advantage of the flexibility of the variational method with coordinate transformations, we derive a self-consistent set of equations of motion from a discretized Lagrangian to study kinetic plasmas using a Fourier decomposed cylindrical coordinate system. The phase-space distribution function was reduced to a collection of finite-sized macro-particles of arbitrary shape moving on a virtual Cartesian grid. However, the discretization of field quantities was performed in cylindrical coordinates and decomposed into a truncated Fourier series in angle. A straightforward finite element interpolation scheme is used to transform between the two grids. The equations of motion were then obtained by demanding the action be stationary. The primary advantage of the variational approach is preservation of Lagrangian symmetries. In the present case, this leads to exact energy conservation, thus avoiding possible difficulties with grid heating.

Plasma, Electromagnetic, Particle-In-Cell, Quasi-3D, Kinetic, Variational, Energy Conserving
###### :
52.65.-y, 52.27.Ny, 52.25.Dg, 52.38.-r, 52.65.Rr
\layoutstyle

8x11single

address=Department of Physics and Astronomy, University of Nebraska-Lincoln, Lincoln, NE 68588-0299, USA address=Department of Physics and Astronomy, University of Nebraska-Lincoln, Lincoln, NE 68588-0299, USA

## 1 Introduction

Several approaches have been used to reduce the computational costs of 3D Particle-in-Cell (PIC) simulations while maintaining predictive power. Some examples include working in the moving window reference frame, using quasi-static approximations [Huang et al.(2006), Mehrling et al.(2014)], using the ponderomotive guiding center algorithm [Gordon et al.(2000)] and the quasi-3D method [Lifschitz et al.(2009), Davidson et al.(2014)]. Lifschitz et al. [Lifschitz et al.(2009)] exploit the near cylindrical symmetry of the laser-plasma interaction by expanding the electromagnetic fields in a small number of poloidal modes while retaining the full 3D dynamics of the macro-particles. By projecting the macro-particle currents onto the grid, all macro-particles in a given annulus contribute current to the same set of grid points, greatly reducing sampling noise.

Recently a variation formulation of macro-particle model has been developed [Evstatiev and Shadwick(2013), Shadwick et al.(2014), Stamm et al.(2014)] that inherits all the benefits of a Lagrangian approach, specifically the connection between symmetries and conservation laws and the flexibility in choosing coordinates. Here we apply this technique to a incorporate a modal decomposition of the electromagnetic fields into a macro-particle model. As in previous work [Evstatiev and Shadwick(2013), Shadwick et al.(2014), Stamm et al.(2014)], we begin with the relativistic Low Lagrangian [Low(1958)] and introduce a distribution function consisting of a collection of fixed-size macro-particles of the form

 f(x,v,t)=Np∑α=1wαSx[x−ξαx(t)]Sy[y−ξαy(t)]Sz[z−ξαz(t)]δ[vx−˙ξαx(t)]δ[vy−˙ξαy(t)]δ[vz−˙ξαz(t)], (1)

where is the macro-particle weight, and are the macro-particle position and velocity respectively, , and are 1D shape functions [Evstatiev and Shadwick(2013)] and is the Dirac delta function. By taking the macro-particle shape to be the product of 1D shape functions, we are implicitly defining a regular Cartesian grid whose parameters are set by the particle sizes in , , and . As a result, it proves convenient to express the vector potential in terms of its Cartesian components, , , and , while parameterizing the functional dependence of these components by , , , and . We then use a truncated Fourier representation in for all the scalar and vector potentials, keeping only the first three modes:

 φ(r,θ,z,t)=φ(o)(r,z,t)+φ(c)(r,z,t)cosθ+φ(s)(r,z,t)sinθ (2)

and likewise for . In most circumstances, this adequately describes the wakefield (the dc term) and the laser field (the first order terms) as discussed in Ref. [Lifschitz et al.(2009)].

## 2 Discrete Lagrangian

We introduce uniform radial and longitudinal grids with grid points and spacing and grid points and spacing , respectively and denote the coordinates of the grid points by , and , . We denote the approximation of the potentials at and by and with . The fully discretized Lagrangian takes the form . The particle contribution can be readily computed from the form of the distribution function [Evstatiev and Shadwick(2013), Stamm et al.(2014)], while the field terms in the Lagrangian can be readily evaluated by representing the spatial derivatives by finite differences:

 Lpart =−mc2Np∑α=1wα ⎷1−(˙ξα)2c2, (3) Lfield Missing or unrecognized delimiter for \Bigg +2c(c)x,lkNr∑a=1D(r)alφ(o)ak+2c(o)x,lk(Nr∑a=1D(r)alφ(c)ak+1rlφ(c)lk)+2c(s)y,lkNr∑a=1D(r)alφ(o)ak+2c(o)y,lk(Nr∑a=1D(r)alφ(s)ak+1rlφ(s)lk) Missing or unrecognized delimiter for \left +2(Nz∑b=1D(z)bkφ(o)lb+1c(o)z,lk)2+(Nz∑b=1D(z)bkφ(c)lb+1c(c)z,lk)2+(Nz∑b=1D(z)bkφ(s)lb1c(s)z,lk)2] − [2(Nz∑b=1D(z)bkA(o)y,lb)2+(Nr∑a=1D(r)alA(o)z,ak−Nz∑b=1D(z)bkA(s)y,lb)2−Nz∑b=1D(z)bkA(o)y,lb(Nr∑a=1D(r)alA(s)z,ak−1rlA(s)z,lk)+(1rlA(s)z,lk)2 +2(Nz∑b=1D(z)bkA(o)x,lb)2+(Nr∑a=1D(r)alA(o)z,ak−Nz∑b=1D(z)bkA(c)x,lb)2−Nz∑b=1D(z)bkA(o)x,lb(Nr∑a=1D(r)alA(c)z,ak−1rlA(c)z,lk)+(1rlA(c)z,lk)2 −12(Nr∑a=1D(r)alA(c)y,ak−1rlA(s)x,lk)(Nr∑a=1D(r)alA(s)x,ak−1rlA(c)y,lk)+(Nz∑b=1D(z)bkA(c)y,lb)2+(Nz∑b=1D(z)bkA(s)x,lb)2 +(Nr∑a=1D(r)alA(o)y,ak)2+(Nr∑a=1D(r)alA(o)x,ak)2+34(Nr∑a=1D(r)alA(c)y,ak−1rlA(s)x,lk)2+34(Nr∑a=1D(r)alA(s)x,ak−1rlA(c)y,lk)2 +(Nr∑a=1D(r)alA(s)z,ak)2+(Nr∑a=1D(r)alA(c)z,ak)2+14(Nr∑a=1D(r)alA(s)y,ak−Nr∑a=1D(r)alA(c)x,ak−1rlA(s)y,lk+1rlA(c)x,lk)2]}, (4)

where and are the finite difference representations of first derivatives in and respectively.

The interaction term, however, requires knowledge of the potentials at arbitrary spatial locations. The most straightforward approach is to treat the macro-particles as moving on a 3D virtual Cartesian grid. The potential on this virtual grid are interpolated between the grid points using finite elements [Evstatiev and Shadwick(2013)], allowing evaluation of the interaction Lagrangian. In terms of the potentials on the virtual grid, the interaction Lagrangian can be written as

 Lint=qNp∑α=1wα∑i,j,kρijk(ξα)[1c˙ξα⋅Aijk−φijk], (5)

where and are the values of the potential on the virtual grid, and is the finite-element projection of
[Evstatiev and Shadwick(2013)]. Assuming a static ion background, can be written as

 Lion=−q\textscI∑i,j,kn\textsc(Ion)ijkφijk, (6)

where the coefficients define the fixed ion charge distribution on the virtual Cartesian grid.

The potentials on the virtual grid are obtained from the Fourier decomposition along with interpolation in (see Fig.2):

 φijk=∑lΛijl(φ(o)lk+φ(c)lkcosθ+φ(s)lksinθ)=∑lΛijl(φ(o)lk+xirijφ(c)lk+yjrijφ(s)lk), (7)

where are the interpolation coefficients, and the index is simply carried through since the -axis is unaffected by the transformation. An analogous expression relates to .

This approach of adopting a virtual grid has general applicability as it decouples the computational grid used to advance the potentials (or fields) from the grid used to interpolate forces. For example, in the case of a non-uniform grid, this approach eliminates the need to have the macro-particle size conform to the local grid spacing.

## 3 Discrete Equations of Motion

The equations of motion are obtained by requiring the action to be stationary under variations of the particle position and potentials. For the particles, the Euler–Lagrange equations give

 dπαdt=qc∑i,j,k[(˙ξα⋅Aijk−cφijk)∂ρijk(ξα)∂ξα−˙Aijkρijk(ξα)−Aijk˙ξα⋅∂ρijk(ξα)∂ξα], (8)

where is the usual relativistic momentum, with . Variation with respect to yields Poisson’s equation, leading to

 Nr∑a=1raD(r)la(1c(c)x,ak+1c(s)y,ak+2Nr∑b=1D(r)baφ(o)bk)+2rlNz∑a=1D(z)ka(1c(o)z,la+Nz∑b=1D(z)baφ(o)lb) Missing or unrecognized delimiter for \right (9) 1rlφ(c)lk+1c(o)x,lk+Nr∑a=1raD(z)la(1c(o)x,ak+Nr∑b=1D(r)baφ(c)bk)+rlNz∑a=1D(z)ka(1c(c)z,la+Nz∑b=1D(z)baφ(c)lb) Missing or unrecognized delimiter for \right (10) 1rlφ(s)lk+1c(o)y,lk+Nr∑a=1raD(z)la(1c(o)y,ak+Nr∑b=1D(r)baφ(s)bk)+rlNz∑a=1D(z)ka(1c(s)z,la+Nz∑b=1D(z)baφ(s)lb) Missing or unrecognized delimiter for \right (11)

Similarly, variation with respect to leads to the following modal wave equations:

## 4 Energy

Since the Lagrangian has no explicit time dependence there is a conserved energy of the form

 W=Np∑α=1˙ξα⋅∂L∂˙ξα+∑l,k⎛⎝˙A(o)lk⋅∂L∂˙A(o)lk+˙A(c)lk⋅∂L∂˙A(c)lk+˙A(s)lk∂L∂˙A(s)lk⎞⎠−L. (21)

Using the equations of motions, it is straightforward to show that this energy is a conserved quantity. After computing the time derivative of the energy, one can use the Lorentz force equation (8) to evaluate the kinetic energy term through giving . The remaining terms in the particle summation can be evaluated by using the modal Poisson and wave equations. Specifically, to evaluate the first term in the brackets, we take the time derivative of the modal Poisson equations (9)–(11) before summing with the corresponding to make use of (7). We apply the same approach to obtain the second term summing the modal wave equation (12)–(20) with . This leaves only field-like terms which cancel to zero, proving that .

## 5 Conclusion

Building on previous work [Evstatiev and Shadwick(2013), Shadwick et al.(2014), Stamm et al.(2014)], we have derived a macro-particle model that exploits the near cylindrical symmetry of the laser-plasma interaction by introducing a poloidal mode expansion of the potentials. We have obtained a set of continuous-time equations of motion for the macro-particles and potentials that exactly conserves energy. As this method obtains the force on the macro-particles directly from the potential, we expect much lower sampling noise than when employing fields in the usual PIC algorithm. Overall, we expect this method to exhibit the same order of computational savings as found by Lifschitz et al. [Lifschitz et al.(2009)]. Future work will involve a non-uniform radial grid, which would be useful for resolving the sheath in the blowout regime and is well adapted to the approach discussed. Additionally, it is desired for computational savings to work in the moving window, which involves a trivial transformation as discussed in Ref. [Stamm et al.(2014)]. Lastly, it would be interesting to develop a symplectic integrator for the time integration of these equations of motion.

###### Acknowledgements.
This work was supported by the U. S. Department of Energy under Contract No. DE-SC0008382 and by the National Science Foundation under Contract No. PHY-1104683.

### References

1. C. Huang, V. K. Decyk, C. Ren, M. Zhou, W. Lu, W. B. Mori, J. H. Cooley, T. M. Antonsen Jr., and T. Katsouleas, J. Comput. Phys. 217, 658–679 (2006).
2. T. Mehrling, C. Benedetti, C. B. Schroeder, and J. Osterhoff, Plasma Phys. Control. Fusion 56, 084012 (2014).
3. D. F. Gordon, W. B. Mori, and T. M. Antonsen, Jr., IEEE Trans. Plasma Sci. PS-28, 1224–1232 (2000).
4. A. Lifschitz, X. Davoine, E. Lefebvre, J. Faure, C. Rechatin, and V. Malka, J. Comput. Phys. 228, 1803 – 1814 (2009).
5. A. Davidson, A. Tableman, W. An, F. S. Tsung, W. Lu, J. Vieira, R. A. Fonseca, L. O. Silva, and W. B. Mori, arXiv:1403.6890 (2014).
6. E. G. Evstatiev, and B. A. Shadwick, J. Comput. Phys. 245, 376–398 (2013).
7. B. A. Shadwick, A. B. Stamm, and E. G. Evstatiev, Phys. Plasmas 21, 055708 (2014).
8. A. B. Stamm, B. A. Shadwick, and E. G. Evstatiev, IEEE Trans. Plasma Sci. 42, 1747–1758 (2014).
9. F. E. Low, Proc. R. Soc. London Ser. A. Math. Phys. Sci 248, 282–287 (1958).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters