AA Magnetic Velocity Verlet Method
A. Chambliss ∗ and J. Franklin † Department of Physics, Reed College, Portland, Oregon 97202, USA
Abstract
We discuss an extension of the velocity Verlet method that accurately approximates the kinetic-energy-conserving charged particle motion that comes from magnetic forcing. For a uniform mag-netic field, the method is shown to conserve both particle kinetic energy and magnetic dipolemoment better than midpoint Runge-Kutta. We then use the magnetic velocity Verlet method togenerate trapped particle trajectories, both in a cylindrical magnetic mirror machine setup, and fordipolar fields like the earth’s magnetic field. Finally, the method is used to compute an exampleof (single) mirror motion in the presence of a magnetic monopole field, where the trajectory canbe described in closed form. a r X i v : . [ phy s i c s . c l a ss - ph ] A ug . INTRODUCTION The motion of particles under the influence of magnetic forces is complicated. We knowthat a particle’s speed does not change but that does not help us predict or visualize thetwists and turns induced by non-uniform fields. Our one concrete analytically tractableexample is motion in a uniform magnetic field, hardly an advertisement for the exotic,spirograph-like trajectories that can appear in general field configurations. Some particletrajectories can be described using special functions, but these are not always accessible tostudents who are unfamiliar with those functions.Absent a closed-form solution, there are various types of predictions one can make aboutthe behavior of particles moving in non-trivial fields. These are almost always either approx-imate or incomplete. As an example of the latter, one of our favorites is Griffiths’s Problem5.43, inviting students to show that a particle that starts at the center of a circular flux-freefield region exits perpendicular to its boundary. That demonstration is essentially an exer-cise in angular momentum conservation. But what if the particle doesn’t exit; what doesit do instead? Another example, the “magnetic mirror,” in which particles are deflectedfrom a region of converging magnetic field lines, is explored later in this article. Beyondthese almost qualitative descriptions, one must use either canned numerical methods, ora homemade implementation of a method that can handle forces that depend on velocity inorder to generate arbitrarily accurate approximations to particle trajectories.One such approach, the Runge-Kutta method, while almost universally applicable tosystems of ODEs, is not informed by the underlying physical problem that it solves. In thispaper we extend a different numerical integrator, the Verlet method, to include magneticforces. This update to the method, while known, is not widely used or taught at theundergraduate level. Its derivation is already of interest to students thinking about theLorentz force, and its ease of implementation, accuracy, and speed make it a desirable toolin a physicists’ numerical toolbox.After deriving the method, we demonstrate its energy conservation superiority over mid-point Runge-Kutta using a uniform magnetic field as a test case. Then we proceed to bothmotivate and numerically generate trajectories for a magnetic mirror configuration. In par-ticular, we use the magnetic moment, the “first adiabatic invariant,” to show that therecan exist oscillatory motion along a field line. Then we generate numerical trajectories that2ealize this oscillatory motion using the magnetic Verlet extension, and test the constancyof the kinetic energy and adiabatic invariant. Having established the numerical existenceof oscillatory trajectories in a field with cylindrical symmetry, we apply the method to adipolar field to see if it can exhibit the oscillatory behavior observed in the cylindrical case.For the magnetic field of a dipole, the equation of particle motion defines “St¨ormer’s prob-lem,” and we can use the magnetic Verlet extension to solve it numerically. The dipolefield can produce particle motion that spirals tightly along a field line while moving up andback along it and a nearby one. As a final example, we demonstrate mirror motion for aparticle moving in a magnetic monopole field, where we have an exact solution with whichto compare.
II. THE VERLET METHOD
The Verlet method generates approximate solutions to Newton’s second law for position-dependent forces. Its simplest variant, position Verlet, can be obtained quickly from Taylorexpansion. The expansion of a particle’s position vector at time t ± ∆ t for small ∆ t is x ( t ± ∆ t ) ≈ x ( t ) ± ˙ x ( t )∆ t + 12 ¨ x ( t )∆ t ±
16 ... x ( t )∆ t + 124 .... x ( t )∆ t ± . . . (1)Adding together x ( t + ∆ t ) and x ( t − ∆ t ) in order to cancel the ˙ x ( t ) term, we get x ( t + ∆ t ) + x ( t − ∆ t ) = 2 x ( t ) + ¨ x ( t )∆ t + O (∆ t ) , (2)where O (∆ t ) means that there are errors of size ∆ t . For a particle of mass m movingunder the influence of a force F ( x , t ), Newton’s second law is m ¨ x ( t ) = F ( x ( t ) , t ), and wecan write ¨ x ( t ) in terms of x ( t ). Then the sum in Eq. (2) can be solved for x ( t + ∆ t ), x ( t + ∆ t ) = 2 x ( t ) − x ( t − ∆ t ) + 1 m F ( x ( t ) , t )∆ t + O (∆ t ) . (3)Viewed as an update method, this equation gives an approximation to the position of aparticle at time t + ∆ t based on its position at the current time t and previous t − ∆ t , bothof which are known. The method is manifestly form-invariant under time-reversal. Moving x ( t + ∆ t ) to the right hand side, and x ( t − ∆ t ) over to the left, we can run the dynamics3ackwards, from t to t − ∆ t : x ( t − ∆ t ) = 2 x ( t ) − x ( t + ∆ t ) + 1 m F ( x ( t ) , t )∆ t + O (∆ t ) , (4)and this mimics the time-reversibility of Newton’s second law itself. The method is anexample of a “symplectic integrator,” a class of numerical methods that preserve certainstructural properties of the underlying Hamiltonian system. Such methods do a good jobconserving total energy, and it is this feature that recommends Verlet to our attention here.If we wanted to test the energy conservation numerically, we would need access to thevelocity at time t in order to construct the kinetic energy. There are a variety of ways toextract that information, but one of the more popular approaches is to re-order the positionupdate in Eq. (3) into separate position and velocity updates. There is nothing new here,just a rearrangement and labelling, leading to the velocity Verlet method defined by theupdates: x ( t + ∆ t ) = x ( t ) + v ( t )∆ t + 12 m F ( x ( t ) , t )∆ t , and (5) v ( t + ∆ t ) = v ( t ) + 12 m [ F ( x ( t ) , t ) + F ( x ( t + ∆ t ) , t + ∆ t )] ∆ t. (6)The velocity approximation is less accurate than the position one, making the method asa whole similar in accuracy to midpoint Runge-Kutta, which makes the same error inboth position and velocity components. In terms of timing, velocity Verlet requires twoevaluations of the force at each update step, just as midpoint Runge-Kutta does. For thesereasons, we will compare the velocity Verlet method to the midpoint Runge-Kutta methodwhen evaluating the former’s numerical properties.The velocity Verlet method is easily applied to conservative forces that depend onlyon position, but unlike Runge-Kutta methods, it is more difficult to introduce velocity-dependent forces like damping. This difficulty is clear from the form of the updates in Eq.(5) and Eq. (6): while the updated positions are known going into the velocity update, ifthere were a force that depended on velocity, F ( x , v , t ), Eq. (6) would become: v ( t + ∆ t ) = v ( t ) + 12 m [ F ( x ( t ) , v ( t ) , t ) + F ( x ( t + ∆ t ) , v ( t + ∆ t ) , t + ∆ t )] ∆ t, (7)and it is not clear how to solve for v ( t + ∆ t ) since it appears on both the left and right-4and sides of this equation. For a pure magnetic force, it is possible to isolate v ( t + ∆ t )algebraically, as will be shown in the next section. III. MAGNETIC VELOCITY VERLET METHOD
Consider a particle with charge q moving in a given magnetic field B ( x ), so that theLorentz force is F ( x , v ) = q v × B ( x ). The velocity update in Eq. (7) can be written v ( t + ∆ t ) = v ( t ) + q ∆ t m [ v ( t ) × B ( x ( t )) + v ( t + ∆ t ) × B ( x ( t + ∆ t ))] . (8)The dependence of the magnetic field on the updated positions, B ( x ( t + ∆ t )), is not aproblem because we have access to those from the position update in Eq. (5). Let’s focuson isolating v ( t + ∆ t ). For visual clarity, define v ≡ v ( t ), B ≡ B ( x ( t )), w ≡ v ( t + ∆ t ), C ≡ B ( x ( t + ∆ t )), and α ≡ q ∆ t/ (2 m ), then the vector update from Eq. (8) reads w = ( v + α v × B ) (cid:124) (cid:123)(cid:122) (cid:125) ≡ d + α w × C , (9)and the term in parenthesis on the right is a constant independent of the target w ; callit d as shown. Dotting C into both sides of Eq. (9) gives w · C = d · C since w × C isperpendicular to C . Crossing C into both sides of Eq. (9), we get w × C = d × C − α (cid:0) w C − C ( d · C ) (cid:1) . (10)The expression on the right depends linearly on w . Using this form for w × C back in theoriginal Eq. (9) yields an equation we can use to isolate w on one side, with known vectorquantities appearing on the other, w = 11 + α C (cid:2) d + α d × C + α C ( d · C ) (cid:3) . (11)5he full velocity Verlet update in this specialized setting is d ≡ (cid:18) v ( t ) + q ∆ t m v ( t ) × B ( x ( t )) (cid:19) , (12) x ( t + ∆ t ) = x ( t ) + d ∆ t, (13) C ≡ B ( x ( t + ∆ t )) , (14) v ( t + ∆ t ) = 11 + (cid:0) q ∆ t m (cid:1) C · C (cid:34) d + q ∆ t m d × C + (cid:18) q ∆ t m (cid:19) C ( d · C ) (cid:35) . (15)If we are given initial values, x (0) = x and v (0) = v , we can use Eq. (13) and Eq. (15) togenerate approximations to x (∆ t ) and v (∆ t ), then use those values to obtain x (2∆ t ) and v (2∆ t ), and so on, up to any desired final time.The pseudocode for the magnetic velocity Verlet method is shown in Algorithm 1 below.You provide the function “MVVerlet” with: 1. the initial position vector, x , 2. the initialvelocity vector, v , 3. the particle mass m , 4. the particle charge q , 5. the time step size ∆ t ,6. the number of steps to take, n , and 7. a function B ( x ) that returns the magnetic field atthe point x . The method returns a list of the particle’s position ( X ) and velocity ( V ) wherethe list index j is associated with time ( j − t for j = 1 → n . Algorithm 1
MVVerlet( x , v , m, q, ∆ t, n, B ) X ← length n list of zeroes X ← x V ← length n list of zeroes V ← v α ← q ∆ t/ (2 m ) x ← x v ← v for j = 2 → n dod ← v + α v × B ( x ) x ← x + d ∆ t C ← B ( x ) v ← ( d + α d × C + α C ( d · C )) / (1 + α C ) X j ← x V j ← vend forreturn { X, V } V. UNIFORM CIRCULAR MOTION
To test the method and display its energy conservation, take a uniform magnetic field B = B ˆ z . A particle of mass m and charge q > x = R ˆ x with initial velocity v = − qB R/m ˆ y and moves in a circle with period T = 2 πm/ ( qB ).We’ll approximate the trajectory for ten cycles, with a time step of ∆ t = T /
50 using boththe magnetic Verlet approach from above and a midpoint Runge-Kutta method with thesame initial conditions and time step. A plot of the kinetic energy for each method as afunction of time is shown in Fig. 1. The difference between the maximum kinetic energyand minimum kinetic energy, divided by the initial kinetic energy, is ∼ − for the Verletmethod, compared with 10 − for midpoint Runge-Kutta. Kinetic Energy (in units of ( qB R ) / (2 m ))time (in units of T = 2 πm/ ( qB )) .
03 10 midpoint RK Verlet
FIG. 1. The kinetic energy as a function of time for a charged particle moving in a uniform magneticfield as calculated by the magnetic Verlet method and the midpoint Runge-Kutta method.
The uniform circular motion here defines a constant magnetic dipole moment. In general,the dipole moment for a current density J ( x ) is defined by the volume integral over all space, m = 12 (cid:90) x × J ( x ) dτ. (16)For a particle at location r ( t ) at time t , the current density is J ( x ) = qδ ( x − r ( t )) ˙ r ( t ),and the dipole moment becomes m = q r ( t ) × v ( t ) /
2. In cylindrical coordinates, { s, φ, z } ,with motion in the xy plane and r ⊥ v as it is here, the moment is m = − ( qv φ s/ z wherethe minus sign comes from the direction of circulation, v = − v φ ˆ φ (clockwise motion). The7niform circular trajectory has s = R and v φ = qB R/m so the particle’s dipole moment is m = − ( qR ) B m ˆ z . (17)The magnitude of the magnetic dipole moment is constant, but how well do the numer-ical methods preserve its value? Working in units of ( qR ) B / (2 m ), the magnitude of themagnetic moment as a function of time is calculated for both the magnetic Verlet methodand the midpoint Runge-Kutta method. The results are shown in Fig. 2, where it is clearthat the Verlet method has dipole magnitude with bounded error as time goes on, whilethe midpoint Runge-Kutta method has error that grows linearly with time, similar to itstreatment of the kinetic energy. Overall, the magnetic velocity Verlet method is superior inpreserving these constants of the motion for this simplest test case. time (in units of T = 2 πm/ ( qB )) .
03 10 midpoint RK Verlet | m | (in units of ( qR ) B / (2 m )) . FIG. 2. The magnitude of the magnetic dipole moment of the particle undergoing uniform circularmotion as a function of time. The midpoint Runge-Kutta method shows linear growth in | m | ,while the Verlet method has much smaller, bounded oscillatory error with period that is the sameas the period of the circular motion. The dipole moment magnitude can be written in terms of the particle’s kinetic energy inthe direction perpendicular to the magnetic field. For the uniform circular motion in thisexample, the perpendicular kinetic energy is K ⊥ = 12 m (cid:18) qB Rm (cid:19) . (18)Comparing this expression with the magnitude of the dipole moment in Eq. (17), it is thecombination K ⊥ /B that captures | m | .In settings where the longitudinal magnetic field does not change magnitude much over8ome region, particle trajectories follow roughly circular motion in the perpendicular plane.Then µ ≡ K ⊥ B , (19)which is constant for uniform B , is approximately constant, and is called the “first adiabaticinvariant.” An equivalent expression for this constant, motivated by uniform circularmotion where it follows from Eq. (19), is µ = 12 qv ⊥ s (20)for s , the radius of the circular motion with speed v ⊥ . We can use the approximatelyconstant value of µ to develop the mirror motion of a particle moving through a non-uniformlongitudinal field. V. CYLINDRICAL MAGNETIC MIRROR MACHINE
Take a longitudinal magnetic field component pointing in the z direction with magnitudethat changes as a function of z , B z ( z ). As a charged particle moves in the z direction, itis presented with a series of field magnitudes that are uniform in the xy plane, and willundergo roughly circular motion in that plane. To get an approximately constant value of µ , the longitudinal motion must be slow enough that multiple cycles of the approximatelycircular motion occur before the longitudinal magnetic field changes significantly. Thisis the “adiabatic assumption,” that the particle’s longitudinal motion is slower than itsperpendicular motion. We will check that this is satisfied, and also that µ is constant, inthe numerical solutions below.A magnetic field cannot consist of B z ( z )ˆ z alone; that violates ∇ · B = 0. In orderto preserve azimuthal symmetry, introduce a magnetic field component pointing in the s direction that depends on s and z , B = B z ( z )ˆ z + B s ( s, z )ˆ s . To make the field divergenceless, B s ( s, z ) must be related to B z ( z ) by B s ( s, z ) = − s dB z ( z ) dz . (21)The longitudinal force on the particle is no longer zero. For the clockwise circulation asso-9iated with a positive charge, v ⊥ = − v φ ˆ φ , and the equation of motion in the z directionis m ¨ z ( t ) = qv φ B s = − qv φ s dB z dz . (22)The term sitting out front in the second equality is precisely µ from Eq. (20) with v ⊥ = v φ ,so that m ¨ z ( t ) = − µ dB z dz . (23)Multiplying both sides of this equation by ˙ z ( t ), m ˙ z ( t )¨ z ( t ) = − µ dB z dz dz ( t ) dt , (24)we can write both sides as total time derivatives provided µ is constant, ddt (cid:18) m ˙ z ( t ) (cid:19) = − ddt ( µB z ( z ( t ))) . (25)The integration is easy to carry out, giving12 m ˙ z ( t ) + µB z ( z ( t )) = C (26)where C is a constant of integration to be set by the initial conditions.If a particle starts off at z (0) = z , with ˙ z (0) = 0, then C = µB z ( z ), and Eq. (26)becomes 12 m ˙ z ( t ) = − µ ( B z ( z ( t )) − B z ( z )) . (27)In order for the longitudinal speed to be real, B z ( z ( t )) must be less than the initial value of B z ( z ). The particle will move towards regions of smaller field, increasing its longitudinalspeed while circulating in the xy plane. If the longitudinal magnetic field magnitude increasestowards the value of B z ( z ) at some location, the particle’s longitudinal speed will decrease.For a location z with B z ( z ) = B z ( z ), the particle must have zero longitudinal speed. Thelongitudinal kinetic energy in Eq. (27) plays a role similar to a one dimensional energyconservation equation in classical mechanics. We can use it to predict some features ofmotion in the longitudinal direction, like turning points and points of maximum speed.Using Eq. (27), it is easy to see how to construct B z ( z ) so as to get oscillatory longi-10udinal motion: Make a field that is symmetric about z = 0 with maxima at ± z . Then aparticle that starts with no longitudinal speed at z will also have zero longitudinal velocitycomponent at − z . So the motion of the particle will be confined between these two points.This “mirror machine” configuration was originally introduced as a way to confine plasmaswithout having mechanical pieces in contact with the plasma. A simple, physically inspired way to achieve a magnetic mirror machine field is to put apair of current loops of radius a carrying steady current I at locations ± d along the z axis.The field produced by this configuration, at z along the axis is B = µ Ia (cid:124) (cid:123)(cid:122) (cid:125) ≡ b ˆ z (cid:34) d − z ) + a ) / + 1(( d + z ) + a ) / (cid:35) , (28)with field magnitude shown in Fig. 3. zB z ( z ) − d d b ! a + 1( a + 4 d ) / " b ( a + d ) / FIG. 3. The magnitude of the magnetic field from Eq. (28). The current rings are at locations ± d , and the maximum and minimum field values depend on both d and the radius of the rings, a ,in addition to the size of the steady current, enapsulated in b . For motion occurring near the z axis, this field is approximately valid, and we’ll take itto define the longitudinal component of the target mirror field. The requirement in Eq. (21)then provides the radial component. Thus, the idealized magnetic field is B = b (cid:34) d − z ) + a ) / + 1(( d + z ) + a ) / (cid:35) ˆ z − bs (cid:34) d − z ( a + ( d − z ) ) / − d + z ( a + ( d + z ) ) / (cid:35) ˆ s . (29)11tarting a particle off at z (0) = z with ˙ z (0) = 0, we expect it to move back and forthbetween ± z . We can take a look at one such trajectory to see how well the kinetic energyand the value of µ are conserved numerically. For the initial position and velocity, take r (0) = R ˆ x + d z (30) v (0) = − qB z ( d/ Rm ˆ y , (31)where R = a/
20 is chosen to be small compared to the radius of the current loops. Withthese initial conditions, if the longitudinal field were uniform with constant value B z ( d/ R with period T = 2 π mqB z ( d/ . (32)Running the magnetic Verlet method for a total time of 100 T in steps of ∆ t = T / z is periodic, going backand forth between ± d/
4. That period of oscillation is much larger than the period of theperpendicular circular motion, with roughly five full cycles along the z axis occurring overthe 100 T time frame. The longitudinal motion is slower than the circular motion, so theadiabatic assumption is satisfied here. Another piece of our assumption was that the radiusof the circular motion does not change much over the course of the trajectory, and that istrue here, as shown in Fig. 4(b).In Fig. 5, we’ve plotted the kinetic energy and value of µ from Eq. (19) using B z ≈ B for the denominator since B s is small. In both of these plots, it is the ratio with the initialvalue that gives us a dimensionless measure of the numerical constancy. The kinetic energy,which is strictly conserved by the equations of motion, sets the numerical standard for a“constant of the motion” with ∼ . µ is approximately conserved exhibiting larger ∼ .
01% change over the 100 T timeframe.The initial velocity in Eq. (31) used a very specific tuning for its perpendicular compo-nent. That velocity was taken to be the value that would ensure uniform circular motion ofradius R in a uniform magnetic field. But changing the initial value of v ⊥ just changes theradius of the circular motion. Starting from an initial position r (0) = ( d/ z with velocity12 (in units of T ) − z ( t ) / ( d/ (a) The particle’s z position as a fractionof the initial z (0) = d/ t (in units of T ) s ( t ) /R . (b) The particle’s radial distance to the z axis as a fraction of the initial radius R .(c) The particletrajectory. FIG. 4. Properties of the numerical solution using the initial conditions in Eq. (30) and Eq. (31). v (0) = − v ˆ x , we can define the time-scale R ≡ mv qB z ( d/ , (33) T ≡ πRv . (34)Then using a time step of ∆ t = T / v = 2 πR/T . The first, with v = α ≡ πR /T for R ≈ .
06 m, has properties shownin Fig. 6, and the second, with a value of v = 3 α , has properties shown in Fig. 7. Thelarger speed leads to a larger radius for the circular motion of the particle shown in Fig. 7.In both cases, the kinetic energy is well conserved, while the first adiabatic invariant showsbetter conservation in the first case where the motion is closer to the z axis.13 (in units of T ) K ( t ) /K (0) t (in units of T ) . . µ ( t ) /µ (0) FIG. 5. The kinetic energy (top) and adiabatic constant µ (bottom), both displayed as fractionsof their initial values. These show that the kinetic energy is quite well conserved, and the valueof µ approximately conserved as a particle oscillates between − d/ d/ z axis whileexecuting (roughly) circular motion in the xy plane. VI. DIPOLE MAGNETIC MIRROR
The spiraling, oscillatory motion we saw in the cylindrical setting persists in more com-plicated field geometries. Imagine bending a magnetic field line from the previous sectionso that it curves. We can still get particle motion that follows the curving field line, circlingaround it, while encountering regions of increasing magnetic field that act as mirrors. Con-sider, for example, a dipolar magnetic field like the one outside the earth. For a magneticdipole pointing in the z direction, the field is B = k π ( x + y + z ) / (cid:2) xz ˆ x + 3 yz ˆ y − ( x + y − z )ˆ z (cid:3) , (35)14 a) The particle trajectory t (in units of T ) . K/K (0) t (in units of T ) . . µ/µ (0) (b) The kinetic energy (top) andadiabatic constant (bottom) as fractionsof their initial values. FIG. 6. Numerical trajectory information for a particle with v = α . where the constant k sets the magnitude. We ran the magnetic Verlet method for a particlewith initial position in the yz plane using a polar angle of θ = 50 ◦ , with initial velocity inthe x direction. The result is shown in Fig. 8, where we can see one “cycle” of the mirroredtrajectory, and then multiple cycles, going all the way around.It is interesting that the particle trajectory follows the field lines, spiraling tightly aroundthem as it goes up and back. This trajectory allows us to map out the dipole field visuallyusing the motion of the charged particle. For electric forcing, field lines point in the directionof acceleration of the particle, and if we start a particle from rest, charges travel along theselines. We build intuition about electric fields using this property. Magnetic field lines areperpendicular to the velocity vector, which points in the direction of motion, so we arenot used to mapping magnetic field lines using the trajectories they generate. The type ofmotion shown in Fig. 8, induced by the magnetic field of the earth, is what causes the VanAllen belts of particles trapped by the earth’s field, traveling up and back along its fieldlines. a) The particle trajectory. t (in units of T ) . K/K (0) t (in units of T ) µ/µ (0) . . (b) The kinetic energy (top) andadiabatic constant (bottom) as fractionsof their initial values. FIG. 7. Numerical trajectory information for a particle with v = 3 α . VII. MONOPOLE FIELD
As our final example, we study the motion of a charged particle in the presence of amagnetic monopole field (or, if you prefer, very close to the north pole of a dipole field).We will again see the behavior that has been the focus of this paper: a charge follows afield line while corkscrewing around it. As the particle encounters a region of increasingfield magnitude, its motion along the field line slows, eventually stopping and reversingdirection along the line. This time, since the field lines converge radially, the geometry ofthe circulation perpendicular to a field line is conical. The magnetic monopole case has theadvantage that there is a closed form expression for the motion of the charge (see Griffiths’Problem 5.45 and references there ), giving us a rare exact result with which to comparethe numerical trajectory.For a magnetic monopole, the field is B = k ˆ r / (4 πr ) where k ≡ µ q m is a constant thatis set by the charge of the monopole. As usual, the kinetic energy of a charged particle16 IG. 8. The mirror effect for a dipolar field. On the left we see a cycle of the particle moving fromthe northern to southern hemisphere and back. On the right, many cycles of the same motion.The sphere is shown just for scale; it has nothing to do with the dipolar field source. moving under the influence of this field is constant. In addition, there is a conserved vector Q ≡ m ( r × v ) − kq π ˆ r , (36)and this vector can be aligned so that it points along the z axis: Q = Q ˆ z . Using theconstancy of Q and the kinetic energy, one can show that θ ( t ) ≡ θ is a constant of themotion, and develop expressions for the time derivatives of r ( t ) and φ ( t ):˙ r ( t ) = ± (cid:115) v − (cid:18) Q sin θmr ( t ) (cid:19) , (37)˙ φ ( t ) = Q mr ( t ) , (38)where v is the constant speed of the particle. Dividing ˙ r ( t ) by ˙ φ ( t ), the ODE governing thespherical r coordinate parametrized by φ is dr ( φ ) dφ = ± (cid:115) v − (cid:18) Q sin θmr ( φ ) (cid:19) (cid:18) mr ( φ ) Q (cid:19) . (39)This equation can be solved, taking the minus sign, r ( φ ) = − Q mv sin θ cos (( φ − α ) sin θ ) , (40)17here α is the constant of integration. We have lost the temporal evolution, but thisequation can be used to draw a picture of the trajectory, and we can compare that withthe solution to the equations of motion that we get numerically from the magnetic Verletmethod. The numerical solution requires a complete set of initial conditions. We can getthose by first choosing the constants of motion, θ (we took θ = 4 ◦ ), α = − Q (negative one, related to the sizeof the magnetic monopole). For the initial value of φ , it is convenient to start at φ = 0, andthen the initial value of r ( φ = 0) is given by the solution in Eq. (40). Finally, the initialvalues for ˙ r ( t = 0) and ˙ φ ( t = 0) can be obtained from Eq. (37) and Eq. (38) respectively,using the rest of the initial coordinate values and constants.The numerical solution is shown plotted as points on top of the positions obtainedfrom Eq. (40) in Fig. 9. There we can see the mirror effect as the particle moves downat first, then reverses direction and goes back up. There is good agreement between the nu-merical solution and the exact one, with the two overlapping. The constant Q is preserved bythe numerical method with norm that varies by only (max( Q ) − min( Q )) / (min( Q )) ≈ − over the portion of the trajectory shown.Unlike the cases we have considered so far, the component of the motion that circlesaround a field line is not cylindrical here, even in approximation. Instead, the charge movesaround a cone with constant polar angle θ , with the tip of the cone at the monopole (shown onthe right in Fig. 9), but this change from cylindrical to conical doesn’t change the qualitativepicture much. We still have a particle that follows a field line while moving around it, slowingand reversing its direction along the line as it encounters increasing field strength. VIII. CONCLUSION
The extension of velocity Verlet to include magnetic fields provides a simple and usefultool for calculating the trajectories of charged particles moving in those fields. While velocityVerlet is less accurate than some other ODE solving techniques (notably higher-order Runge-Kutta methods), it is easy to implement, and its derivation highlights some of the vectorgeometry associated with the Lorentz force in the context of Newton’s second law. Onecan easily extend the method to include additional position-dependent forces. Going backto Eq. (8), an additional force ¯ F ( x ) would introduce a term like ¯ F ( x ( t )) + ¯ F ( x ( t + ∆ t )),18 tartmagneticmonopolestart FIG. 9. A particle starts (as shown) moving towards a magnetic monopole along a radially directedfield line. The mirror effect, produced by converging field lines, causes the particle to changedirection and head away from the monopole. On the left, the points are the numerical solution,with the solid gray line the solution from Eq. (40). On the right, we see the cone superimposedwith the numerical solution. The tip of the cone is at the monopole, and it has the opening angleused in the numerical solution, θ = 4 ◦ . where both force evaluations rely on quantities that are known by the time the update isperformed, so that these terms would get lumped into the vector d in Eq. (9), at whichpoint the derivation proceeds as in the text. The method is fast, and can handle multipleparticles (making it a favorite of molecular dynamics solvers ). We hope this paper servesto increase its use in the undergraduate E&M curriculum. ACKNOWLEDGMENTS
The authors thank David Griffiths for useful commentary and for suggesting the monopolemagnetic field example. ∗ [email protected] † [email protected] Jane M. Repko, Wayne W. Repko, Allan Saaf, “Charged particle trajectories in simple non-uniform magnetic fields,”
Am. J. Phys. , , 652–655 (1991). David J. Griffiths,
Introduction to Electrodynamics , 4th ed. (Cambridge University Press, 2017). M. Kaan ¨Ozt¨urk, “Trajectories of charged particles trapped in Earth’s magnetic field,”
Am. J.Phys. , , 420–428 (2012). George C. McGuire, “Using computer algebra to investigate the motion of an electric charge inmagnetic and electric dipole fields,”
Am. J. Phys. , , 809–812 (2003). Elisha R. Huggins, Jeffrey J. Lelek, “Motion of electrons in electric and magnetic fields; intro-ductory laboratory and computer studies,”
Am. J. Phys. , , 992–999 (1979). J. Franklin, K. C. Newton, “Classical and quantum mechanical motion in magnetic fields,”
Am.J. Phys. , (4), 263–269 (2016). Q. Spreiter, M. Walter, “Classical molecular dynamics simulation with the Velocity Verletalgorithm at strong external magnetic fields,”
J. Comp. Phys. , 1, 102–119 (1999). John David Jackson,
Classical Electrodynamics , 3rd ed. (Wiley, 1998). Herbert Goldstein, Charles P. Poole Jr. and John L. Safko,
Classical Mechanics , 3rd ed. (Pear-son, 2001). Loup Verlet, “Computer ‘Experiments’ on Classical Fluids. I. Thermodynamical Properties ofLennard-Jones Molecules,”
Phys. Rev. (1), 98–103 (1967). Denis Donnelly, “Symplectic Integrators: An introduction,”
Am. J. Phys. , (10), 938– 945(2005). M. P. Allen, D. J. Tildesley,
Computer Simulations of Liquids , (Oxford University Press, NewYork, 1987). Joel Franklin,