Modeling laser wakefield accelerators in a Lorentz boosted frame

Modeling laser wakefield accelerators in a Lorentz boosted frame J.-L. Vay1,4, C. G. R. Geddes1, E. Cormier-Michel2, D. P. Grote3,4 1Lawrence Berkeley National Laboratory, CA, USA 2Tech-X Corporation, CO, USA 3Lawrence Livermore National Laboratory, CA, USA 4Heavy Ion Fusion Science Virtual National Laboratory, USA Annapolis, MD, USA • June 13-19, 2010

Outline • Scaling, limitations, numerical Cerenkov • Novel numerical algorithms implemented in Warp • Application to: • scaled 10 GeV stage • full scale 10 GeV stage • full scale 100 GeV- 1 TeV stages • Conclusion

x/t= (L/l, T/t) FB-rest frame of “B” 0 2 20 y y 0 x x 0 0 0 Range of space and time scales spanned by two identical beams crossing each other* F0-center of mass frame space space+time •  is not invariant under the Lorentz transformation: x/t . • There exists an “optimum” frame which minimizes it. *J.-L. Vay, Phys. Rev. Lett.98, 130405 (2007)

laser plasma wake w e-beam Scaling of LPA boosted frame simulations speedup Laser plasma accelerator = plasma length/dephasing length = fraction of wake existing plasma In boosted frame at (,), the speedup S is: Example • Low  • Max S Near maximum speedup reached when ≈w Over four orders of magnitude theoretical speedup for a 10 GeV stage.

x/t= (L/l, T/t) laser plasma wake 0 w e-beam 2 20 0 0 0 0 Scaling for laser plasma accelerator is very similar to generic rigid identical beams case Two rigid identical beams Laser plasma accelerator Illustrates commonality of underlying principle to the crossing of relativistic matter and/or light.

Using conventional PIC techniques, 2-3 orders of magnitude speedup reported in 2D/3D by various groups Osiris: trapped injection* Vorpal: external injection w/ beam loading** Warp: external injection wo/ beam loading *Martins et al, UCLA/IST (talk TuPlenary) **Bruhwiler et al, Tech X (talk WG2 TuAM) • Reported speedups limited by various factors: • laser transverse size at injection, • statistics (self-trapped injection), • short wavelength instability (most severe).

Limitations from laser launching Osiris Vorpal Warp Courtesy S. Martins, IST Laser is initialized at once in simulation box. In a boosted frame at , the largest effective spot size increases as (1+)2*, i.e. the transverse box surface increases as (1+)24! Laser is emitted from all but one surface. Laser is emitted from a moving plane. Blowup of transverse box size can be avoided by using novel injection procedures for injecting the laser. *Martins et al., Comp. Phys. Comm., 2010

Instability reported for large  Warp 2D full scale simulation 10 GeV LWFA stage (ne=1017cc, =130) Longitudinal electric field laser plasma Side view laser A short wavelength instability is observed at front of plasma for large  Plasma moving near speed of light: numerical Cerenkov? plasma Top-rear view

Effect of spatial resolution x=z=13 m x=z=6.5 m x=z=3.25 m FFT FFT FFT • Amplitude of instability decreases as resolution increases, • transverse mode governed by physics, longitudinal by numerics.

Numerical Cerenkov? Numerical Cerenkov is due to macro-particles traveling > numerical speed of light in vacuum. • Numerical dispersion of Yee solver in 3D: • numerical speed of light in vacuum is 2/3 actual speed of light at Nyquist limit along main axes. x=y=z) • First identified by Boris (JCP 1973) and Haber et al (Proc. ICNSP 1973). • Studied by Godfrey (JCP 1974 & 1975) and Greenwood et al (JCP 2004). • Proposed cures: • higher shape factors for particles and filtering (commonly used), • damping of electromagnetic field (Friedman JCP 1990 identified as most by Greenwood et al), • solver with larger stencil (difficulty for preserving Gauss Law), • FFT-based solver proposed by Haber (boundary conditions difficult).

1/4 1/2 1/4 Bilinear filter 1/4 1/2 1/4 Bilinear filter with stride 2 Efficient filtering using strides Multiple pass of bilinear filter + compensation routinely used 100% absorption at Nyquist freq. Bilinear (B) Bilinear (B) +compensation (C) 4*B+C (B4C) Using a stride N shifts the 100% absorption frequency to Fnyquist/N B4C stride 1 (G1) B4C stride 2 (G2) B4C stride 3 (G3) B4C stride 4 (G4) • Combination of filters with strides allows for more efficient filtering: • G1G2  20*B+C; speedup=2 • G1G2G3  50*B+C; speedup=3.5 • G1G2G4  80*B+C; speedup=5.5

* Cole-Karkkainen solver adapted to the Particle-In-Cell methodology Cole* and Karkkainen** have applied Non-Standard Finite-Difference (NSFD) to source free Maxwell equations Warp: switched FD/NSFD to B/E. => FD on source terms, i.e. standard exact current deposition schemes still valid. NSFD: weighted average of quantities transverse to FD ().    a      FD Yee Cole-Karkkainen x=y=z) NSFD Cole-Karkkainen allows for larger time step, resulting in perfect dispersion along axes. * J. B. Cole, IEEE Trans. Microw. Theory Tech., 45 (1997). J. B. Cole, IEEE Trans. Antennas Prop., 50 (2002). ** M. Karkkainen et al., Proc. ICAP, Chamonix, France (2006).

More generally, NSFD-based solver offers tunability of numerical dispersion Yee Cole-Karkkainen (CK) CK2 intermediate CK3 CK4 CK5 more compact (=0) perfect dispersion 2D diag. most isotropic Solver can be tuned to better fit particular needs.

Application of Friedman damping straightforward B push modified to with where is damping parameter. Yee-Friedman (YF) Cole-Karkkainen-Friedman (CKF) Dispersion degrades with higher values of  Damping more potent on axis and more isotropic for CKF than YF.

PML implemented with tunable solver Same high efficiency as with Yee.

Scaled 10 GeV stage Physical and numerical parameters of scaled 10 GeV stage* (relevance to BELLA**) Full scale 1017 cm-3 1.5 m Wake frame:≈13 *Similar to E. Cormier-Michel et al, AAC 08; C. Geddes et al, PAC 09 ** W. Leemans et al, talk MoPlenary

Percent level agreement demonstrated on scaled simulation of 10 GeV stage Beam energy Momentum spread at peak energy Timing 2D 3D Yee solver; no damping; standard filtering; Max speedup ≈200 in 2D; 130 in 3D

Different frames = different views of same physics 2D simulation of 100 MeV stage Laboratory frame E (laser) E// Wake frame E (laser) E// Amount of short wavelength content much reduced in wake frame Fixed station diagnostics confirm that physics is the same 3D

Calculating in boosted frame allows for more damping, filtering and for less constraint on cell aspect ratio Friedman damping Wideband filtering Yee vs CK (with cubic cells) Beam energy For a given level of damping or filtering, the accuracy is improved for higher frame boost, best in wake frame. In wake frame, Yee solver as accurate as CK solver with cubic cells. Timing

Simulations in 2D of full scale 10 GeV stage Conclusion regarding damping and filtering holds • Yee solver • In wake frame (≈130): • square cells • ct=z/√2 =130 Laser/2 • No instability • % level agreement in =30-130 range

Simulations in 3D of full scale 10 GeV stage in wake frame CK solver; in wake frame (≈130): cubic cells, ct=z. =0. =0.1 =0.5 =130 =130 =130 • Strong instability (worse with Yee solver) • Filtering more potent than damping for control of instability

Simulations in 3D of full scale 10 GeV stage Conclusion regarding damping and filtering unchanged • CK2 solver • In wake frame (≈130): • cubic cells • ct=z/√2 =130 • No instability • Good agreement in =30-130 range

Simulations in 2D and 3D with same time step • 2D • Yee solver • In wake frame (≈130): • square cells • ct=z/√3 • 3D • Yee solver • In wake frame (≈130): • cubic cells • ct=z/√3 • Strong instability at similar level • Similar spectrum in 2D and 3D • Is the time step a key parameter?

Time step scan 2D; CK solver; in wake frame (≈130); square cells; 0.5 ≤ ct/z ≤ 1. Sharp decrease of instability level at ct=z/√2!

Other dependencies? 2D; CK solver; in wake frame (≈130); square cells; 0.5 ≤ ct/z ≤ 1. Node-centered field gathering (“momentum conserving”) Yee-mesh field gathering (“energy conserving”) A clue to the mystery: existence of “magical” time step depends on field gathering procedure.

Main findings • can apply high levels of filtering in wake frame • instability growth reduced by orders of magnitude when • ct=z/√2 “magical” time step • (provided that one uses Yee mesh centered field gathering)* • *This result was recently confirmed by D. Bruhwiler using Vorpal.

Modeling of 100 Gev-1 TeV stages demonstrated using magical time step + filtering • 2D • Yee solver • In wake frame: • square cells • ct=z/√2 • 3D • CK2 solver • In wake frame: • cubic cells • ct=z/√2 • 0.1-10 GeV: filter S(1) • 100 GeV: filter S(1:2) • 1 TeV: filter S(1:2:3) • 0.1-1 GeV: filter S(1) • 10-100 GeV: filter S(1:2)

Over 5-6 orders of magnitude speedup for 100 GeV-1 TeV stages Osiris: trapped injection Vorpal: external injection w/ beam loading Warp: external injection wo/ beam loading 100 GeV stage in 3D: only 4 hours using 2016 CPU of Cray Franklin at NERSC

Conclusion • Speedup of Lorentz boosted frame simulations was limited by: • transverse expansion of injected laser: Speedup≤102 in 3D • short wavelength instability: Speedup≤104 in 2D; ≤103 in 3D • Limitation from laser injection removed with special injection procedures • Algorithms implemented in Warp for instability control • efficient filtering using strides • novel field solver with tunable numerical dispersion (including PML) • Friedman damping • Analysis of boosted frame simulations of 10 GeV stages reveals that • damping and filtering can be used more aggressively in wake frame • there is a magical time step for which the instability is minimized • field gathering procedure matters

Conclusion (2) • Instability does not seem to be Numerical Cerenkov • interest in CK field solver resides in allowing CFL limit ≥ magical time step • First time verification of scaling of deeply depleted stages up-to 1 TeV • Applied to external injection wo/ beam loading so far • Application to case w/ beam loading underway and trapped injection pending looking promising

Warp: a parallel framework combining features of plasma (Particle-In-Cell) and accelerator codes Parallel scaling of Warp 3D PIC-EM solver on Franklin supercomputer (NERSC) • Geometry:3D (x,y,z), 2D-1/2 (x,y), (x,z) or axisym. (r,z) • Python and Fortran:“steerable,” input decks are programs • Field solvers:Electrostatic - FFT, multigrid; AMR; implicit • Magnetostatic - FFT, multigrid; AMR; implicit • Electromagnetic - Yee, Cole-Kark.; PML; AMR • Parallel:MPI (1, 2 and 3D domain decomposition) • Boundaries:“cut-cell” --- no restriction to “Legos” • Lattice:general; non-paraxial; can read MAD files • solenoids, dipoles, quads, sextupoles, linear maps, arbitrary fields, acceleration • Bends:“warped” coordinates; no “reference orbit” • Particle movers:Boris, large time step “drift-Lorentz”, novel relativistic Leapfrog • Reference frame:lab, moving-window, Lorentz boosted • Surface/volume physics: secondary e-/photo-e- emission, gas emission/tracking/ionization • Diagnostics:extensive snapshots and histories • Misc.:trajectory tracing; quasistatic & steady-flow modes; space charge emitted emission; “equilibrium-like” beam loads in linear focusing channels; maintained using CVS repository. 32,768 cores

Extras

More challenging: large intense beam; plasma taper Early Ramp on plasma density provides phase locking Later (Lab frame) - beam not matched transversely and experienced losses, => significant population of “halo” particles, - filtered out of the diagnostics using cutoffs in position and energy. Very good agreement between runs in frames at =1-10

Widening the band of the digital filter improves significantly filter applied to current density and gathered field, NOT on Maxwell EM field Smoothing 1 Smoothing 2 Smoothing 3 Smoothing 4 E// (GV/m) E// (GV/m) E// (GV/m) E// (GV/m) Z Z Z Z laser Result is only slightly affected, even with most aggressive filtering tested. wideband filter applied along longitudinal direction only

Other possible complication: inputs/outputs • Often, initial conditions known and output desired in laboratory frame • relativity of simultaneity => inject/collect at plane(s)  to direction of boost. • Injection through a moving plane in boosted frame (fix in lab frame) • fields include frozen particles, • same for laser in EM calculations. • Diagnostics: collect data at a collection of planes • fixed in lab fr., moving in boosted fr., • interpolation in space and/or time, • already done routinely with Warp for comparison with experimental data, often known at given stations in lab. frozen active z’,t’=LT(z,t)

H H E E E E E E Is the instability due to Yee solver numerical dispersion errors?Implementation of a low-dispersion solver in Warp Enlarged stencil* => no disp. in x,y,z Implementation in Warp • E and H switched, • => E push same as Yee, • exact charge conservation preserved in 2D & 3D with unmodified Esirkepov current deposition and implied enlarged stencil on div E. Stencil on H unchanged 1-2-1 filter needed at t=x/c 1-D simulation of current Heaviside step • Issues for PIC: • source terms are not given, • odd-even oscillation when t=x/c. full amplitude odd-even oscillations for t=x/c 1-D simulation of LWFA at t=x/c 121 filter No filter *M. Karkkainen, et al., Proceedings of ICAP’06, Chamonix, France

n+1Vn+1 + nVn qt m 2 n+1/2 Vn+1 + Vn qt m 2 (with , , , , , , ). Seems simple but ! . Algorithms which work in one frame may break in another. Example: the Boris particle pusher. • Boris pusher ubiquitous • In first attempt of e-cloud calculation using the Boris pusher, the beam was lost in a few betatron periods! • Position push: Xn+1/2 = Xn-1/2 + Vn t -- no issue • Velocity push: n+1Vn+1 = nVn + (En+1/2 + Bn+1/2) issue: E+vB=0 implies E=B=0=> large errors when E+vB0 (e.g. relativistic beams). • Solution • Velocity push: n+1Vn+1 = nVn + (En+1/2 + Bn+1/2) • Not used before because of implicitness. We solved it analytically* *J.-L. Vay, Phys. Plasmas15, 056701 (2008)

level 0 level 1 level 2 level 3 Mesh Refinement expected to offer further savings Example of a 2D Warp simulation of LPA stage with up-to 3 levels of mesh refinement*. E// (V/m) PIC MR-PIC Higher resolution Lower resolution • Successful initial benchmarking: beam emittance history from calculation at various resolutions with standard PIC (one grid) recovered with MR-PIC runs, • resolution varied transversely only; further tests pending. *Vay et al, IPAC 2009.

Modeling laser wakefield accelerators in a Lorentz boosted frame