Comparison of Numerical Schemes to Solve the Newton-Lorentz Equation Jabus van den Berg1 1 Centre for Space Research, North-West University, Potchefstroom, South Africa 24182869@nwu.ac.za Abstract To better understand the transport of solar energetic particles, a conceptual understanding of the micro-physics of charged particle propagation in electric and magnetic ﬁelds is needed. The movement of charged particles are governed by the Newton-Lorentz equation, which becomes increasingly diﬃcult to solve analytically for complicated electric and magnetic ﬁelds. In this work, the methods of Boris (1970) and Vay (2008), as well as a forth-order Runge-Kutta scheme, are investigated for numerically solving the Newton-Lorentz equation. This work focuses on accuracy rather than computational speed, since these numerical schemes will be used to analyse the motion of particles in turbulent electric and magnetic ﬁelds. The Larmor radius, deviation between the ﬁnal numerical and analytical position, as well as the execution time, are calculated and recorded for diﬀerent time steps. 1 The Newton-Lorentz Equation The motion of a particle with mass m and charge q, moving with velocity ⃗v in an electric ﬁeld ⃗E and magnetic ﬁeld ⃗B, is governed by the Newton-Lorentz equation md(γ⃗v) dt = q  ⃗E + ⃗v × ⃗B  , (1) where γ = 1/ p 1 −(v/c)2 is the Lorentz factor, with c the speed of light in vacuum (Griﬃths, 1999). This can be split into two, ﬁrst-order diﬀerential equations d⃗r(t) dt = ⃗v(t) and md⃗u(t) dt = q h ⃗E (⃗r(t); t) + ⃗v(t) × ⃗B (⃗r(t); t) i , (2) where ⃗u(t) = γ(t)⃗v(t). This is solved in Cartesian coordinates since the vector product is easily calculated and it is diﬃcult to ﬁnd a symmetry axis for a general magnetic ﬁeld. The numerical methods used to solve this equations should be accurate to simulate charged particles in turbulent electric and magnetic ﬁelds, as shown in Fig. 1. The Larmor radius of the particle’s gyration is deﬁned as rL = v⊥/ωc, (3) where v⊥= v sin[arccos(⃗v · ⃗B/vB)] is the particle’s speed perpendicular to the magnetic ﬁeld and ωc = |q|B/γm (4) is the cyclotron frequency. s 1 3 2 1 0 1 s 2 3 2 1 0 1 2 3 s 0 2 4 6 8 x [m] 6 5 4 3 2 1 0 1 y [m] 1 0 1 2 3 4 5 z [m] 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Fig. 1: Left panel: Background (black) and turbulent (green) magnetic ﬁeld lines of a composite slab-2D turbulence model at time t = 0. The magnetic ﬁeld lines can be though of as the path of the particle’s guiding centre if the particle is stuck to the ﬁeld line. Right panel: Trajectory (red) of an electron in the composite turbulence model, simulated using the Runge-Kutta method. A single background magnetic ﬁeld line (black dashed) and the particle’s guiding centre (blue dotted) is also shown. 2.1 The Method of Boris (1970) The position can be calculated by a time-centred approximation of Eq. 2a ⃗rn+1/2 −⃗rn−1/2 ∆t ≈⃗vn =⇒ ⃗rn+1/2 = ⃗rn−1/2 + ⃗vn∆t, (5) where a subscript n denotes that quantity at time tn = t0 + n∆t, with t0 the initial time and n = 0; 1; 2; . . . , and it should be remembered that ⃗vn = ⃗un/γn. Boris (1970) suggested a decoupling of the electric and magnetic ﬁeld to solve Eq. 2b by ﬁrst applying half of the electric ﬁeld, then applying the rotation of the velocity vector by the magnetic ﬁeld and lastly applying the other half of the electric ﬁeld: ⃗u−= ⃗un−1 + q∆t 2m ⃗En−1/2 ⃗un+1/2 = ⃗u−+ f1⃗u−× ⃗Bn−1/2 (6a) ⃗u+ = ⃗u−+ f2⃗un+1/2 × ⃗Bn−1/2 ⃗un+1 = ⃗u+ + q∆t 2m ⃗En−1/2, (6b) where f1 = tan qBn−1/2∆t/2mγ−  /Bn−1/2, (7a) f2 = 2f1/  1 + (f1Bn−1/2)2 and γn = p 1 + (un/c)2, (7b) with ⃗En−1/2 = ⃗E(⃗rn−1/2; tn−1/2), ⃗Bn−1/2 = ⃗B(⃗rn−1/2; tn−1/2) and Bn−1/2 = |⃗Bn−1/2|. 2.1.1 An Approximation Since tan θ ≈θ to ﬁrst order in θ if |θ| ≪1, as should hold for a small time step, Eqs. 6a to 7b can be simpliﬁed (Birdsall & Langdon, 1991): ⃗u−= ⃗un−1 + q∆t 2m ⃗En−1/2 ⃗un+1/2 = ⃗u−+ ⃗u−× ⃗f1 (8a) ⃗u+ = ⃗u−+ ⃗un+1/2 × ⃗f2 ⃗un+1 = ⃗u+ + q∆t 2m ⃗En−1/2, (8b) where ⃗f1 ≈q ⃗Bn−1/2∆t/2mγ− and ⃗f2 = 2⃗f1/(1 + f 2