High order integration schemes to solve the  Maxwell's equations for
the fields and Newton's equations for the particles are presented

A. Sgattoni

C. Bendetti

G. Turchetti

P. Londrillo

High order integration schemes to solve the  Maxwell's equations for
the fields and Newton's equations for the particles are presented

SGATTONI, ANDREA

TURCHETTI, GIORGIO

Archivio istituzionale della ricerca - Alma Mater Studiorum Università di Bologna

High order integration schemes  for particle in cell (PIC)

The detailed numerical investigation of the highly nonlinear physics involved in the interaction of a laser-pulse with a plasma and/or an externally injected beam requires suitable simulation tools which are able to retain the basic features of the process without increasing too much the computational needs. In our group at the University of Bologna we developed an electromagnetic PIC code using high-order schemes (explicit/compact finite differences derivatives and a 4th-order Runge-Kutta) instead of the most used second-order accurate algorithm, this allows to reach higher accuracy with similar number of grid-points or to use coarser grid and relax the computer load without compromising the accuracy. Here we present some validation tests and comparisons for the different schemes

Sgattoni, A.

Londrillo, P.

Benedetti, C.

Turchetti, G.

Scientific Open-access Literature Archive and Repository

DOI 10.1393/ncc/i2009-10418-4Colloquia: CSFI 2008IL NUOVO CIMENTO Vol. 32 C, N. 2 Marzo-Aprile 2009High-order integration schemes for Particle In Cell (PIC) methodA. Sgattoni(1), P. Londrillo(2), C. Benedetti(1) and G. Turchetti(1)(1) Dipartimento di Fisica, Universita` di Bologna and INFN, Sezione di BolognaVia Irnerio 46, 40126 Bologna, Italy(2) INFN, INAF Bologna Osservatorio Astronomico - Bologna, Italy(ricevuto il 4 Agosto 2009; pubblicato online il 9 Ottobre 2009)Summary. — The detailed numerical investigation of the highly nonlinear physicsinvolved in the interaction of a laser-pulse with a plasma and/or an externally in-jected beam requires suitable simulation tools which are able to retain the basicfeatures of the process without increasing too much the computational needs. In ourgroup at the University of Bologna we developed an electromagnetic PIC code usinghigh-order schemes (explicit/compact finite differences derivatives and a 4th-orderRunge-Kutta) instead of the most used second-order accurate algorithm, this allowsto reach higher accuracy with similar number of grid-points or to use coarser gridand relax the computer load without compromising the accuracy. Here we presentsome validation tests and comparisons for the different schemes.PACS 52.38.kd – Laser-plasma acceleration of electrons and ions.PACS 52.65.-y – Plasma simulation.PACS 02.60.Jh – Numerical differentiation and integration.1. – IntroductionMost existing PIC codes to study the interaction of ultra intense and short laser pulsewith plasma are implemented on the basis of classical schemes [1, 2], where: i) time in-tegration is discretized by second-order leap-frog schemes and space derivatives by thesecond-order finite difference on staggered mesh (Yee’s module) [1,2]; ii) shape functionsfor charge weighting and field interpolation are linear or quadratic B-splines, and iii)charge conservation is usually implemented by the Esirkepov procedure [3]. To improvethis basic numerical framework on accuracy and efficiency, as needed by the challengingnew physical phenomena arising in high-intensity laser-plasma interaction, the Bolognagroup started a new project aimed at developing a PIC code based on high-order integra-tion schemes (HOPIC). Consolidated results of this numerical activity have been imple-mented and tested on the new code ALaDyn [4,5]. Main features of the ALaDyn code are i)a fourth-order Runge-Kutta (RK4) algorithms for time integration, ii) a high-order (upto eighth) compact/explicit finite differences scheme for space integration; iii) a stretchedc© Societa` Italiana di Fisica 261262 A. SGATTONI, P. LONDRILLO, C. BENEDETTI and G. TURCHETTIcomputational grid in the transverse direction in order to ensure a high resolution whereneeded (usually around the propagation axis) together with substantially far boundaries;iv) a hierarchical numerical particles sampling that allow to put more particles in the“dynamically interesting” zones and, at the same time, decrease the particles in the out-lying regions. In the PIC framework computational parameters controlling accuracy aregiven by: i) space-time resolution measured by the number of grid points which sam-ple the smallest relevant spatial scale to be resolved (nλ = λ/Δx = σλ/Δt, where λ isthe relevant spatial scale, Δx is the grid size, Δt is the related time step and σ is theCourant number); ii) the number of representative numerical particles per cell nc and iii)the order of the shape functions. Numerical experiments using classical (second-order)PIC schemes show a very slow convergence to high-accuracy results by increasing thecontrol parameters [6]. HOPIC are expected to allow: i) lower numerical errors, for fixedspace-time resolution and the number of particle per cell (ppc); ii) faster convergence ratewhen increasing space-time resolution. In the following, we document on the improvedaccuracy of a HOPIC approach by considering 1D non-trivial test problems where somereference analytical solution is available.2. – Numerical integration of Maxwell equationsCompact finite differences schemes for space differentiation [7], coupled to a RK4scheme for time integration assure efficient and highly accurate numerical solution forthe Maxwell equations. In a one-dimensional grid with uniform spacing h and node pointsxj = jh, j = 0, 1, N − 1, the compact first derivative u′(xj) of a u(xj+1/2) function withstaggered collocation is expressed by(1) u′ ≡ Dˆ[u] =[Pˆ−1CˆΔu],where Pˆ = [α, 1, α] and Cˆ = [b, (a+b), b] are tridiagonal matrix with coefficients a = 3(3−2α)/8, and b = (22α−1)/24 and [Δu]j = [uj+1/2−uj−1/2]/h is the two-point second-orderexplicit derivative. This one-parameter (α) family of compact schemes gives for α = 9/62,a sixth-order scheme (SC6) and for α = 1/22, b = 0 a fourth-order scheme (SC4). Usingthe same formula, the α = 0 case reproduces the fourth-order explicit scheme (SE4). Theclassical second-order derivative (SE2) is here expressed by the limiting case α = b = 0and a = 1. By expressing u(x) in terms of Fourier modes eikx, where k is the discretizedwave number, the derivative Dˆ matrix has purely imaginary eigenvalues iZ(w)/h wherew = kh, and Z(w) is the numerical or modified wave number, replacing the spectral(exact) ik numerical derivative. The one-dimensional wave equation in semi-discreteform reads ∂tu(x, t) + Dˆu(x, t) = 0 gives a dispersion relation ω(k) = Z(w)/h and anumerical wave speed vph = ω/k = Z(w)/w. The difference R = Z(w)/w − 1 providesthen a measure of the resolution property of a derivative scheme. There are significantdifferences between the SE2 and high-order compact derivative (SC): to keep the phaseerror below some tolerance value, say R ≤ 10−4, a value of points per wavelength (PPW)larger than 40 is necessary for SE2, while a number of PPW in the range 10–12 is alreadysufficient for SC4 and SC6. In a fully discretized wave equation, using a fourth-orderRunge-Kutta (RK4) scheme for time integration with time step Δt = σh, the resultingdispersion relation for the complex frequency Ω = ω + iγ, takes the formHIGH-ORDER INTEGRATION SCHEMES FOR PARTICLE IN CELL (PIC) METHOD 263Fig. 1. – The error in phase speed in the 1D wave equation as a function of PPW parameter,using composite schemes: (a) (RK4 + SC6), (b) (RK4 + SC4), (c) (LP2 + SE2).(2) eγΔt = |G| =√G2R + G2I , ω(k) =kσwsin−1GI|G| ,where G(w) = [GR, GI ] is the complex amplification factor, with(3) GR(w) = 1− 12(σZ(w))2 +124(σZ(w))4, GI(w) = σZ(w)− 16(σZ(w))3.The stability condition |G| ≤ 1 entails an upper limit for the Courant number σ  1for compact schemes SC4 and SC6. The corresponding dispersion relation for a Leap-Frog scheme (LPF2) is given by ω(k) = 2kσw sin−1 σZ(w)/2 with local stability conditionsσZ(w) < 2, where now Z(w) is the modified wave number of the SE2 derivative scheme.In fig. 1 we plot the phase speed of a 1D wave freely propagating versus the number ofpoints per wavelength.2.1. Multidimensional issues. – In the multidimensional case, high-order compactschemes having spectral-like behaviors along each coordinate direction Zs(w) ≈ ws s =x, y, z give accurate isotropy properties of the numerical Laplacian operator(4) Z2 = Z2x + Z2y + Z2z ,Fig. 2. – The error in phase speed in the 2D wave equation as a function of PPW, for differentpolar angle θ. Left: (RK4 + SC6) scheme; right: (LPF2 + SE2) scheme.264 A. SGATTONI, P. LONDRILLO, C. BENEDETTI and G. TURCHETTIFig. 3. – The maximum error in the numerical wave solution of the (LPF2+SE2) and (RK4+SC4) schemes, for different resolution PPW = λ/Δz parameter. The corresponding referencescaling laws in the Δt = σΔz grid size are also reported.while the reference (LPF2+SE2) scheme has dispersion relation strongly dependent onthe polar angle with respect to the propagation direction. This can be appreciated infig. 2 where we have plotted the dispersion relations at different angles.2.2. Testing 1D wave propagation. – Together with the analytical consideration, wetested the different schemes in propagating a 1D wave packet. To measure the error inphase speed for different numerical schemes, we consider a linear polarized wave packetpropagating along the z-axis, defined by the vector potential Ax(z, t) with analyticalform(5) Ax(z, t) = a0 cos(ω0(z − ct)) cos4(β(z − ct)), Ex(z, t) = By(z, t) = −1c∂tAx(z, t),with ω0 = 2πλ , β =π2τ , where λ and τ are, respectively, the wavelength and the durationof the pulse. We compare then the cumulative error, by considering maxz |Enumx (z, t)−Ex(z, t)| after a fixed evolution time T = 50λ/c. Figure 3 shows the maximum error fordifferent integration schemes, when increasing the resolution parameter PPW.3. – Integration of particles motionTo test accuracy in the particle motion, we consider test particles in the freely prop-agating (Ex, By) fields given in eq. (5), with cold px = pz = 0 initial conditions.In the reference analytical model, each particle preserves the two integrals of motion:h1(t) = p2x/2 − pz h2(t) = Px = Ax − px giving hk(t) = hk(0) = 0, k = 1, 2. In thediscretized problem, the total numerical error results as a combination of: i) the errorinduced by the time integration scheme of the particle motion; ii) the numerical errorinduced by the space-time integration scheme wave equation; iii) the interpolation errorto evaluate (Ax, Ex, By) fields at the particle position zi(t).We compare the LPF2 − SE2 and the RK4 − SC4 space-time integrators. Thecomparison, see fig. 4, shows significantly different errors and very different scaling lawswith respect to the resolution.3.1. Fully self-consistent problem. – When considering the fully self-consistent prob-lem, the errors in the particle motions manifest themselves as an anomalous numericalHIGH-ORDER INTEGRATION SCHEMES FOR PARTICLE IN CELL (PIC) METHOD 265Fig. 4. – Scaling of the maximum error |h(t) − h(0)| of the particle invariants of motion atincreasing space-time resolution Δt = σΔz for PPW = 16, 24, 32, 48, σ = 0.6.heating causing spurious injections of electrons in the plasma wave (numerical wavebreaking), another consequence is the incorrect value of the plasma wavelength. Figure 5shows the nonlinear generation of 1D plasma wave by an intense laser pulse, as can beseen, at low resolution both, high- and low-order algorithms show spurious wave break-ing (the high spike in the phase space) that disappears increasing the resolution. WhileRunge-Kutta results show a very fast conversion rate, the Leap-Frog results show a muchslower one and the need of a high number of grid points. This results in a reward for theincreased number of operation per time step needed for the Runge-Kutta: the possibilityto relax the resolution by a factor of roughly 2 in each dimension, i.e. a factor 8 in 2Dand a factor 16 in 3D less operation (the time step is related to the spacial grid-step)gives a significant advance to higher-order scheme.Fig. 5. – Longitudinal phase space of the electrons in the wake of a 1D wave packet. Left:leap-frog with different spatial/temporal resolution. Right: RK4+SE6. The figure shows theplasma wave generated in the wake of the pulse. For LPF at least 40 PPW are needed to get acorrect phase space, while for RK4 20 PPW give already a good accuracy.266 A. SGATTONI, P. LONDRILLO, C. BENEDETTI and G. TURCHETTIREFERENCES[1] Birdsall C. K. and Langdom A. B., Plasma Physics via Computer Simulation (McGraw-Hill, New York) 1985.[2] Hockney R. W. and Eastwood J. W., Computer Simulation Using Particles (IOP) 1988.[3] Esirkepov T. Zh., Comput. Phys. Comm., 135 (2001) 144.[4] Benedetti C., Sgattoni A., Turchetti G. and Londrillo P., IEEE Trans. PlasmaSci., 36/4 (1985) 1790.[5] Benedetti C., Londrillo P., Petrillo V., Serafini L., Sgattoni A., Tomassini P.and Turchetti G., Nucl. Instrum. Methods Phys. Res. A, 608 (2009) S94.[6] Cormier-Michel E., Shadwick B. A., Geddes C. G. R., Esarey E., Schroeder C. B.and Leemas W. P., Phys. Rev. E, 78 (2008) 016404.[7] Lele S. K., J. Comput. Phys., 103 (1992) 16.

English

High-order integration schemes for Particle In Cell (PIC) method

http://eprints.bice.rm.cnr.it/16576/1/ncc9422.pdf

High order integration schemes for particle in cell (PIC)

Abstract

Similar works

Full text

Available Versions

Archivio istituzionale della ricerca - Alma Mater Studiorum Università di Bologna

Scientific Open-access Literature Archive and Repository