Tutorial Computation for Advanced Accelerators Thomas M Antonsen

Tutorial: Computation for Advanced Accelerators Thomas M. Antonsen Jr, IREAP, University of Maryland College Park MD, 20742 AAC 2012 Work supported by NSF, and DOE

Computation has become integral to Accelerator Science Why is Computation Important ? - Understanding of basic physical processes - Understanding and diagnosing particular experiments - Designing improved experiments - Optimizing designs

Areas of Progress Improved Physics Models Find better descriptions of what’s going on inside a an accelerator More Efficient Algorithms Solve the governing equations with fewer operations and/or less memory usage. New Computing Architectures Restructure your code to take advantage of new hardware.

Examples of Areas where Computation is Important Particle trajectories in analytic fields* Space charge fields* Magnetic fields due to currents in coils RF fields in cavities and waveguides Radiation generated by energetic particles Multipactor Plasma-Based Acceleration and more *See R. Ryne http: //uspas. fnal. gov/materials/09 UNM/Computational. Methods. pdf

Meshes and Grids Calculations are carried out on digital computers Time is discretized Unstructured mesh Space may also be discretized - meshes or grids If field is a dynamical variable (Maxwell’s Equations) a mesh is needed If field describes pair-wise interaction forces (Coulombs law) it can be calculated with or without a grid. Structured mesh

Tree Code Calculation of Coulomb interaction without grids J. Barnes and P. Hut, “A hierarchical force-calculation algorithm, ” Nature, vol. 324, pp. 446– 449, 1986 Direct application of Coulomb’s law to N particles requires O(N 2) operations. Not practical for large N. Tree code: Divide domain into cells Each cell has fewer than NL particles Particles within same cell interact pairwise Particles in different cells interact through multipole moments O(Nlog. N) operations possible

Mesh Types Structured Mesh Unstructured Mesh Complex geometries may be represented using a multi-block representation: conformal meshes: boundary of mesh conforms with physical boundary

Existing Discretization Techniques u Finite Integration Technique Applicable to orthogonal Cartesian mesh Degrees of freedom are integrals of fields over appropriate mesh elements u Finite Element Method Applicable to curvilinear mesh (structured or unstructured) Field described by interpolating basis functions in each cell Basis & test functions 1 -½ Suffers due to non-conformal representation of geometry 0 ½ Solves integral (weak) form of 2 nd order differential equations Evaluates matrix elements using overlap integrals of basis & test functions Large number of coupling terms

Fields and Structured Mesh for FDTD solution of Maxwell’s Equations Yee Mesh (1966) • Widely used • Simple (explicit) • Second order accurate • Stable if Dt <cdx Electric fields defined on edges Magnetic fields defined on faces Electric fields defined at integer time steps Magnetic fields defined at half integer time steps

Leap Frog Explicit FDTD Solutions to Maxwell’s Equation time differencing Alternately apply Faraday’s law to faces Explicit: knowing E at time t allows for direct calculation of H at time t+Dt/2

Finite Differencing Plane Waves Assume all fields vary as Dispersion relation for Plane waves Continuous FDTD

Stability - Courant Freidricks Lewy Condition (CFL) Dispersion Relation Left side < Stable if CFL Condition Maximum of Right Side

Dispersion and Accuracy Dispersion in discrete system is periodic. Departures from continuous system are quadratic in Dt, D. Waves have phase velocity <c Numerical Cerenkov Radiation - Requires filter

More complicated evaluation of curl

Conducting Boundaries on Structured Grids When boundary does not coincide with mesh FDTD must be modified. Boundary Without modification errors of order D are introduced. S. Dey, R. Mittra, A locally conformal finite-difference time-domain (FDTD) algorithm modeling three-dimensional perfectly conducting objects, IEEE Microwave Guided Wave Lett. 7 (1997) 273 Problem: If cells are small CFL condition is violated Separate modification to compute surface losses.

ADI (Alternating Dimension Implicit) Solutions to Maxwell’s Equation) Scalar: Douglas, Peaceman, and Ratchford (1955) ME: Namiki T (1999), Zheng, Chen, Zhang (2000) time differencing • Split time steps • Recently developed • Complex (implicit) • Not subject to Courant condition, Dt > cdx Second half of time step has x, y, and z permuted

Leapfrog ADI-FDTD equations S. Cooke et al (2009) Tri-diagonal operators u u u “curl” terms Electric and magnetic field terms are offset in time Electric fields only appear on the half time steps Magnetic fields only appear on the full time steps Still need to solve six tri-diagonal equations per time step 3 for the E-components, 3 for the H-components ADI-FDTD uses 6 for the E-components (3 components, 2 half steps) No additional explicit equations are necessary

PIC Algorithm - From Plasma Physics via Computer Simulation C. K. Birdsall and A. B. Langdon FDTD ME in here

Integration of Equations of Motion Basic Leap Frog: Modified for Lorenz Force: “Boris Push” p and x known on staggered time grid Half step with E Full step with B Half step with E

Interpolation of Fields and Accumulation of Currents and Charges E and J are known on the same grid, but are staggered in time Forces interpolated from grid to location of particles Particle contribution to current density accumulated on grid. Algorithm should conserve charge. J. Villasenor and O. Buneman. Computer Physics Communications, 69. 306 -316, (1992)

21 IVEC 2012 FDTD-PIC Time-stepping loop Cooke NRL H E J P X S. Charge is not computed, only current is used. Space charge fields occur only through integration of Maxwell’s equations. Advance H Assign external current sources Add new particles Advance P, then X (accumulating J) Remove dead particles Charge-conserving methods do this exactly – no need for divergence-cleaning. Advance E

Plasma Based Accelerators • Direct Acceleration in Modulated Channels

Physical Processes Important • Relativistic Motion • Plasma Wave Generation / Cavitation • Laser Self Focusing / Scattering • Ionization Not Important • Turbulence - most plasma particles interact for a short time, several plasma periods, then leave Our job is “relatively” easy. Simulation has played an important role in the development of the field. Guiding our understanding Designing new experiments

Relevant Time Scales (LWFA) • Laser Period =800 nm • Driver Duration ~ Plasma Period TL = 2. 7 fs TD = 50 fs • Propagation Time (driver evolution, L= 1 cm) TP = 3. 3 104 fs (driver evolution, L= 1 m) TP = 3. 3 106 fs Propagation Time >> Driver Duration >> Laser Period TP >> TD >> TL • Disparity in time scales leads to both complications and simplifications

MODELS APPROACHES / APPROXIMATIONS • Laser Full EM - Laser Envelope Particles - Fluid • Plasma Full Lorenz force - Ponderomotive Dynamic response - Quasi-static

Hierarchy of Descriptions 1. Full Format Particle in Cell (PIC) - Relativistic equations of motion for macro-particles - Maxwell’s Equations on a grid - Most accurate but most computationally intensive • Laser Full EM - Laser Envelope • Plasma Particles - Fluid Full Lorenz force - Ponderomotive Dynamic response - Quasi-static

Windows and Frames Three versions: - Lab Frame - Moving Window - Boosted Frame t=L/c Lab Frame - # Grids NLF = (L/dz)2 ~ (L/ )2 grid must resolve dt=dz/c - LD Pulse length z=0 Moving Window - # Grids NMW = (LD/L) NLF z=L Propagation distance

Boosted Frame J-L Vay, PRL 98, 130405 (2007) Lab Frame t=L/c dt' = gdt Boost - g t' =L/gc z' =-g. LD z=0 z=L Qualifications: -no backscatter. In boosted frame gives upshift. Requires smaller dz' -transverse grid size sets limit on dt' < dx' /c= dx /c z'=0 dz' = gdz Boosted frame - # grids NBF = g-2 NMW Optimum based on 1 D makes simulation “square”

Using conventional PIC techniques, 2 -3 orders of magnitude speedup reported in 2 D/3 D by various groups Osiris: trapped injection Vorpal: external injection w/ beam loading Warp: external injection wo/ beam loading J-L Vay Reported speedups limited by various factors: • laser transverse size at injection, • statistics (trapped injection), • short wavelength instability (most severe).

2. Full EM vs. Laser Envelope Driver Duration >> Laser Period TD >> TL • Required Approximation for Laser envelope: wlaser tpulse >> 1, rspot >> wp /wlaser <<1 • Advantages of envelope model: -Larger time steps Full EM stability: Dt < Dx/c Envelope accuracy: Dt < 2 p Dx 2/ c - No unphysical Cherenkov radiation - Further approximations • Advantages of full EM: Includes Stimulated Raman back-scattering Also direct acceleration

Laser Envelope Approximation • Laser + Wake field: E = E laser + E wake = A (x , t) exp ik x • Vector Potential: A laser 0 0 • Laser Frame Coordinate: + c. c. x = ct – z • Envelope equation: Necessary for: Raman Forward Self phase modulation vg< c Drop (eliminates Raman back-scatter)

Validity of Envelope Equation Extended Paraxial – True: Requires : 2 k c 2 + w 2 p v g(w ) = c / 1 + 2 w 2 k 2 c 2 + w 2 p v g(w ) = c 1 – 2 1/2 w 2 k c 2 , w 2 p < < w 2 Extended Paraxial approximation - Correct treatment of forward and near forward scattered radiation - Does not treat backscattered radiation - If grid is dense enough can treat Dw ~ w

3. Full Lorenz Force vs. Ponderomotive Description • Full Lorenz: d p i v ´B = q E+ i c dt • Separation of time scales E = E laser + E wake = x(t) + x(t) i dx = vi dt • Requires small excursion × ÑE x(t) laser < < Elaser • Ponderomotive Equations v ´ Bwake d p = q E wake + + Fp c dt 2 q A laser Fp = – mc Ñ mc 2 2 g g = q A laser p 2 1+ 2 2 + mc mc 2 2 2

Ponderomotive Guiding Center PIC Code: Turbo. WAVE D. F. Gordon, et al. , IEEE TRANSACTIONS ON PLASMA SCIENCE , V 28 , 1224 -1232 ( 2000 ) Fields are separated into high and low frequency components. The low frequency component is treated as in an ordinary PIC code. The high frequency component is treated using a reduced description which averages over optical cycles. Low Frequency Cycle High Frequency Cycle Deposit Sources j and r Deposit Source nq 2/<gm> Advance Fields (Maxwell’s Equations) Advance Laser (Envelope Equation) Lorentz Push + Ponderomotive Push

INF&RNO (INtegrated Fluid & pa. Rticle simulatio. N c. Ode) C. Benedetti at al. (LBNL) 2 D cylindrical + envelope for the laser (ponderomotive approximation) full PIC/fluid description for plasma particle (quasi-static approx. is also available) switching between PIC/fluid modalities (hybrid PIC-fluid sim. are possible) dynamical particle resampling to reduce on-axis noise nd 2 order upwind/centered FD schemes + RK 2/RK 4 (& Implicit) for time integration linear/quadratic shape functions force interpolation/charge deposition high order low-pass compact filter for current/field smoothing “BELLA”-like runs (10 Ge. V in ~ 1 m) become feasible in a few days on small machines FLUID PIC 35

4. Quasi - Static vs. Dynamic Wake P. Sprangle, E. Esarey and A. Ting, Phys. Rev. Lett. 64, 2011 (1990) Laser Pulse Plasma. Wake t pulse Electron transit time: t e = 1 – vz / c c Plasma electron c - vz Trapped electron Electron transit time << Pulse modification time d = ¶ + c –v ¶ +v ×Ñ z dt ¶t ¶x Advantages: fewer particles, less noise (particles marched in ct-z) Disadvantages: particles are not trapped

Quasi Static Simulation Code WAKE P. Mora and T. M. Antonsen Jr. - Phys Plasma 4, 217 (1997) x=ct-z Two Time Scales 1. “fast time” t ~ TD ~ wp-1 particle trajectories and wake fields determined 2. “slow time” laser pulse evolves diffraction self-focusing depletion r Particle trajectories

Quick. PIC: 3 D quasi-static particle code Description • Massivelly Parallel, 3 D Quasi-static particle-incell code • Ponderomotive guiding center for laser driver • 100 -1000+ savings with high fidelity • Field ionization and radiation reaction included • Simplified version used for e-cloud modeling • Developed by UCLA + UMaryland + IST • • • Examples of applications Simulations for PWFA experiments, E 157/162/164 X/167 (Including Feb. 2007 Nature) Study of electron cloud effect in LHC. Plasma afterburner design up to Te. V Efficient simulation of externally injected LWFA Beam loading studies using laser/beam drivers Chengkun Huang: huangck@ee. ucla. edu http: //exodus. physics. ucla. edu/ http: //cfp. ist. utl. pt/golp/epp New Features • Particle tracking • Pipelining • Parallel scaling to 1, 000+ processors • Beta version of enhanced pipelining algorithm: enables scaling to 10, 000+ processors and unprecedented simulation

QUASI-STATIC CODE STRUCTURE Laser “slow time - t” “fast time x=ct-z” Particles and Wake

PARTICLES CONTINUED constant of motion : • Hamiltonian: • Weak dependence on “t” in the laser frame Transverse Dynamics • Introduce potentials • Algebraic equation:

LOCAL FREQUENCY MODIFICATION M. Tsoufras, Ph. D Thesis, UCLA Intensity ne Frequency shift is proportional to propagation distance Assume complete cavitation Relative frequency shift is unity at the dephasing distance

Quasi Static Modeling of Pulse Depletion W. Zhu Wednesday WG 1

Comparison of Spectra; Modified Paraxial, Full-wave, and FDTD-PIC

Conclusions • Numerical simulation of Laser-Plasma interactions is a powerful tool • A variety of models and algorithms exist first principles reduced modles • Field is still advancing with new developments Boosted Frame Calculations speed-up ~ g 2 3 D Parallel Quasi-static speed-up ~ [w 0/wp]2 GPUs

WAKE - Cavitation P. Mora and T. Antonsen PHYSICAL REVIEW E, Volume: 53 R 2068 (1996) Intensity Density Complete cavitation Suppression of Raman instability Stable propagation for 30 Rayleigh lengths Trajectories

Wake simulation of pulse propagation in corrugated plasma channels See WG #1 B. Layer Wed. PM, J. Palastro Thu. AM =800 nm wch=15 m ao=. 35 t=1 ps no=7 x 1018 cm-3 =. 9 m=. 035 cm 3 cm ~ 90 corrugations ~ 40 TR 150 m density intensity t=0 TR Upper: corrugated Lower: smooth t=15 TR t=30 TR Suppression of Raman Instability

QUICKPIC Quasi-Static Plasma Particles Wake Fields Dynamic Beam Particles Laser Field

Second Order Accurate Split-Step Scheme time-t Laser Advance t Axial (dt) Transverse (dt) t+dt Axial (dt) Transverse (dt) Axial (dt/2) Plasma Particle Advance push particles time-s t t+dt

Verification : Full PIC vs. Quasi-static PIC Benchmark for different drivers: Quick. PIC vs. Full PIC Excellent agreement with full PIC. 100 to 10000 times savings in CPU needs e- driver with ionization e+ driver No noise issues and no unphysical Cerenkov radiation laser driver 100 to 10000 CPU savings with “no” loss in accuracy

Iteration of Electromagnetic Field Electromagnetic portion Parallel electric field Transverse electric field Iterate to find Ampere’s law Equation of motion

Quasi-Static Field Representation Lorenz Pro: Simple structure Compatible with 2 D PIC Con: A carries “electrostatic” field Transverse Coulomb Pro: A = 0 in electrostatic limit Con: non-standard field equations

Particle Promotion in WAKE S. Morshed Po. P, to be published Osiris Wake “Plasma particles” for which quasi -static violated are promoted to “beam particles”

New “Leapfrog” ADI-FDTD Formulation u Combine a half step forward with a half step back: u Subtracting equations gives a new scheme: Simple source terms Tri-diagonal operator “curl H” term Compare to FDTD: S. J. Cooke, M. Botton, T. M. Antonsen, Jr. , B. Levush, “A leapfrog formulation of the 3 -D ADI-FDTD algorithm” Int. J. Numer. Model. 2009; 22: 187– 200 54

NRL Plasma Physics Division Turbo. WAVE Framework PIC EM pusher, ponderomotive pusher, gather/scatter Nonlinear Optics “SPARC” Multi-species hydrodynamics Anharmonic Lorentz model, Quasistatic model Field Solver Modules Chemistry Explicit, envelope, direct fields, coulomb gauge, electrostatic. Lindman boundaries, simple conducting regions, PML media. Arbitrary species, reactions Numerical Framework Linear algebra, elliptical solvers, fluid advection… Turbo. WAVE Base Grid geometry, regions, domain decomposition, diagnostics, input/output…

Images from Turbo. WAVE Gridded Gun Modeling Contours are equipotential lines. Colors represent current density. Black circles represent grid wires. NRL Plasma Physics Division Blowout Wakefield Red = ion rich Blue = electron rich ve k 0 Utilizes 3 D electromagnetic PIC Utilizes 2 D electrostatic PIC (slab or cylindrical) with conducting regions

Envelope model in VORPAL enables full 3 D simulations of meter-scale LPA stages Speedup of 50, 000 shown for meter-scale parameters u Rigorously tested, showing second-order convergence and correct laser group velocity u Benchmarks against scaled Ben Cowan simulations show excellent [1] B. Cowan et al. , Proc. AAC 2008 agreement [2] P. Messmer and D. Bruhwiler, PRST-AB u 9, 031302 (2006) [3] D. Gordon, IEEE Trans. Plasma Sci. 35, 1486 (2007) [4] P. Mora and T. M. Antonsen, Phys. Plasmas 4, 217 (1997)

Charge Conserving Algorithm

Full Format: PIC Algorithm

VORPAL

OSIRIS 2. 0 osiris framework • • Massivelly Parallel, Fully Relativistic Particle-in-Cell (PIC) Code Visualization and Data Analysis Infrastructure Developed by the osiris. consortium: UCLA + IST Widely used: UCLA, SLAC, USC, Michigan, Rochester, IST, Imperial College, Max Planck Inst. • Examples of applications • • Mangles et al. , Nature 431 529 (2004). Tsung et al. , Phys. Rev. Lett. , 94 185002 (2004) Mangles et al. , Phys. Rev. Lett. , 96, 215001 (2006) Lu et al. , Phys. Rev. ST: AB, 10, 061301 (2007) I. A. Blumenfeld et al. , Nature 445, 241 (2007) • • • Frank Tsung: tsung@physics. ucla. edu Ricardo Fonseca: ricardo. fonseca@ist. utl. pt http: //exodus. physics. ucla. edu http: //cfp. ist. utl. pt/golp/epp New Features in v 2. 0 • • Bessel Beams Binary Collision Module Tunnel (ADK) and Impact Ionization Dynamic Load Balancing PML absorbing BC Optimized Higher Order Splines Parallel I/O (HDF 5) Boosted Frame in 1/2/3 D

+10 Ge. V self-injection in nonlinear regime Controlled self-guided a 0=5. 8 Smooth accelerating field Injected electrons Laser pulse Boosted frame 7000 x 256 cells ~109 particles 3 x 104 timesteps Υ=10 Samuel F. Martins et al. Nature Physics V 6, April 2010 ~300 x faster than lab 7 -12 Ge. V 1 -2 n. C