title: A Magnetic Velocity Verlet Method
authors: A. Chambliss J. Franklin
text: 
  
  
  
  
  chambgam@alumni.reed.edu
  
  
  jfrankli@reed.edu
  Department of Physics, Reed College, Portland, Oregon 97202,
  USA
  
  
  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 magnetic field, the method is shown to conserve both particle kinetic energy and magnetic dipole moment better than midpoint Runge-Kutta.  We then use the magnetic velocity Verlet method to generate trapped particle trajectories, both in a cylindrical magnetic mirror machine setup, and for dipolar fields like the earth's magnetic field.  Finally, the method is used to compute an example of (single) mirror motion in the presence of a magnetic monopole field, where the trajectory can be described in closed form.
  
  
  A Magnetic Velocity Verlet Method
      J. Franklin
      December 30, 2023
  =================================
  
  
  
  
  
  
  § INTRODUCTION
  
  
  The motion of particles under the influence of magnetic forces is complicated.  We know that a particle's speed does not change but that does not help us predict or visualize the twists and turns induced by non-uniform fields.  Our one concrete analytically tractable example is motion in a uniform magnetic field, hardly an advertisement for the exotic, spirograph-like trajectories that can appear in general field configurations.  Some particle trajectories can be described using special functions,<cit.> but these are not always accessible to students who are unfamiliar with those functions.
  
  
  Absent a closed-form solution, there are various types of predictions one can make about the behavior of particles moving in non-trivial fields.  These are almost always either approximate or incomplete.  As an example of the latter, one of our favorites is Griffiths's Problem 5.43,<cit.> inviting students to show that a particle that starts at the center of a circular flux-free field region exits perpendicular to its boundary.  That demonstration is essentially an exercise in angular momentum conservation.  But what if the particle doesn't exit; what does it do instead? <cit.>  Another example, the “magnetic mirror," in which particles are deflected from a region of converging magnetic field lines, is explored later in this article.  Beyond these almost qualitative descriptions, one must use either canned numerical methods, <cit.> or a homemade implementation of a method that can handle forces that depend on velocity in order to generate arbitrarily accurate approximations to particle trajectories.
  
  
  One such approach, the Runge-Kutta method, while almost universally applicable to systems of ODEs, is not informed by the underlying physical problem that it solves.  In this paper we extend a different numerical integrator, the Verlet method, to include magnetic forces.  This update to the method, while known, <cit.> is not widely used or taught at the undergraduate level.  Its derivation is already of interest to students thinking about the Lorentz force, and its ease of implementation, accuracy, and speed make it a desirable tool in a physicists' numerical toolbox.
  
  
  After deriving the method, we demonstrate its energy conservation superiority over midpoint Runge-Kutta using a uniform magnetic field as a test case.  Then we proceed to both motivate and numerically generate trajectories for a magnetic mirror configuration. <cit.>  In particular, we use the magnetic moment, the “first adiabatic invariant," <cit.> to show that there can exist oscillatory motion along a field line.  Then we generate numerical trajectories that realize this oscillatory motion using the magnetic Verlet extension, and test the constancy of the kinetic energy and adiabatic invariant.  Having established the numerical existence of oscillatory trajectories in a field with cylindrical symmetry, we apply the method to a dipolar 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örmer's problem," and we can use the magnetic Verlet extension to solve it numerically.  The dipole field can produce particle motion that spirals tightly along a field line while moving up and back along it and a nearby one.  As a final example, we demonstrate mirror motion for a particle moving in a magnetic monopole field, where we have an exact solution with which to compare.
  
  
  
  § THE VERLET METHOD
  
  
  
  The Verlet method generates approximate solutions to Newton's second law for position-dependent forces.  Its simplest variant, position Verlet, <cit.> can be obtained quickly from Taylor expansion.  The expansion of a particle's position vector at time t ±Δ t for small Δ t is
  
      x(t ±Δ t) ≈ x(t) ±ẋ(t) Δ t + 1/2ẍ(t) Δ t^2 ±1/6⃛ x(t) Δ t^3 + 1/24⃜ x(t) Δ t^4 ±…
  
  Adding together x(t+Δ t) and x(t - Δ t) in order to cancel the ẋ(t) term, we get
  
      x(t + Δ t) +  x(t - Δ t) = 2  x(t) + ẍ(t) Δ t^2 + O(Δ t^4),
  
  where O(Δ t^4) means that there are errors of size Δ t^4.  For a particle of mass m moving under the influence of a force F( x, t), Newton's second law is m ẍ(t) =  F( x(t),t), and we can write ẍ(t) in terms of x(t).  Then the sum in PVpre can be solved for x(t+Δ t),
  
      x(t+Δ t) = 2  x(t) -  x(t - Δ t) + 1/m F( x(t),t) Δ t^2 + O(Δ t^4).
  
  Viewed as an update method, this equation gives an approximation to the position of a particle at time t + Δ t based on its position at the current time t and previous t - Δ t, both of 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 dynamics backwards, from t to t - Δ t:
  
      x(t-Δ t) = 2  x(t) -  x(t + Δ t) + 1/m F( x(t),t) Δ t^2 + O(Δ t^4),
  
  and this mimics the time-reversibility of Newton's second law itself.  The method is an example of a “symplectic integrator," a class of numerical methods that preserve certain structural properties of the underlying Hamiltonian system. <cit.>  Such methods do a good job conserving 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 the velocity at time t in order to construct the kinetic energy.  There are a variety of ways to extract that information, but one of the more popular approaches is to re-order the position update in posVerlxup into separate position and velocity updates.  There is nothing new here, just a rearrangement and labelling, leading to the velocity Verlet method defined by the updates: <cit.> 
  
      x(t + Δ t)    =    x(t) +  v(t)Δ t  + 1/2 m F( x(t),t) Δ t^2,  and 
       v(t + Δ t)    =     v(t) + 1/2 m[  F( x(t),t) +  F( x(t + Δ t), t+Δ t)] Δ t.
  
  The velocity approximation is less accurate than the position one, making the method as a whole similar in accuracy to midpoint Runge-Kutta, <cit.> which makes the same error in both position and velocity components.  In terms of timing, velocity Verlet requires two evaluations of the force at each update step, just as midpoint Runge-Kutta does.  For these reasons, we will compare the velocity Verlet method to the midpoint Runge-Kutta method when evaluating the former's numerical properties.
  
  
  The velocity Verlet method is easily applied to conservative forces that depend only on 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 vvdefx and vvdefv:  while the updated positions are known going into the velocity update, if there were a force that depended on velocity, F( x,  v, t), vvdefv would become:
  
      v(t + Δ t) =  v(t) + 1/2 m[  F( x(t), v(t), t) +  F( x(t + Δ t),  v(t + Δ t),  t+Δ t)] Δ t,
  
  and it is not clear how to solve for v(t + Δ t) since it appears on both the left and right-hand 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.
  
  
  
  
  § MAGNETIC VELOCITY VERLET METHOD
  
  Consider a particle with charge q moving in a given magnetic field B( x), so that the Lorentz force is F( x,  v) = q  v× B( x).  The velocity update in vimp can be written 
  
      v(t + Δ t) =  v(t) + q Δ t/2 m[  v(t) × B( x(t)) +   v(t + Δ t) × B( x(t + Δ t)) ].
  
  The dependence of the magnetic field on the updated positions, B( x(t + Δ t)), is not a problem because we have access to those from the position update in vvdefx.  Let's focus on 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 LorentzUp reads
  
      w =(  v + α v× B)_≡ d + α w× C,
  
  and the term in parenthesis on the right is a constant independent of the target w; call it d as shown.  Dotting C into both sides of wstart gives w· C =  d· C since w× C is perpendicular to C.  Crossing C into both sides of wstart, we get
  
      w× C  =  d× C - α(  w C^2 -  C ( d· C)).
  
  The expression on the right depends linearly on w.  Using this form for w× C back in the original wstart yields an equation we can use to isolate w on one side, with known vector quantities appearing on the other,
  
      w = 1/1 + α^2 C^2[  d + α d× C + α^2  C ( d· C) ].
  
  The full velocity Verlet update in this specialized setting is
  
      d   ≡    v(t) + q Δ t/2 m v (t) × B( x(t)),  
       x(t + Δ t)    =    x(t) +  dΔ t   , 
       C   ≡    B( x(t + Δ t)), 
       v(t + Δ t)    =   1/1 + q Δ t/2 m^2  C· C[ d + q Δ t/2 m d× C +   q Δ t/2 m^2  C d· C] .
  
  If we are given initial values, x(0) =  x_0 and v(0) =  v_0, we can use mvxup and mvvup to generate 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_0, 2. the initial velocity vector, v_0, 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 at the point x.  The method returns a list of the particle's position (X) and velocity (V) where the list index j is associated with time (j-1) Δ t for j = 1 → n.
  
  
  
  
  § UNIFORM CIRCULAR MOTION
  
  
  
  To test the method and display its energy conservation, take a uniform magnetic field B = B_0 ẑ.  A particle of mass m and charge q > 0 starts at initial position x_0 = R x̂ with initial velocity v_0 = -q B_0 R/m ŷ and moves in a circle with period T = 2 π m/(q B_0).  We'll approximate the trajectory for ten cycles, with a time step of Δ t = T/50 using both the magnetic Verlet approach from above and a midpoint Runge-Kutta method with the same initial conditions and time step.  A plot of the kinetic energy for each method as a function of time is shown in fig:KEcomp.  The difference between the maximum kinetic energy and minimum kinetic energy, divided by the initial kinetic energy, is ∼ 10^-14 for the Verlet method, compared with 10^-2 for midpoint Runge-Kutta.
  
  
  
  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 = 1/2∫ x× J( x) dτ.
  
  For a particle at location r(t) at time t, the current density is J( x) = q δ^3( x -  r(t)) ṙ(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 = -(q v_ϕ s/2) ẑ where the minus sign comes from the direction of circulation, v = -v_ϕϕ̂ (clockwise motion).  The uniform circular trajectory has s = R and v_ϕ = q B_0 R/m so the particle's dipole moment is 
  
      m =-(q R)^2 B_0/2 mẑ.
  
  
  
  The magnitude of the magnetic dipole moment is constant, but how well do the numerical methods preserve its value?  Working in units of (q R)^2 B_0/(2 m), the magnitude of the magnetic moment as a function of time is calculated for both the magnetic Verlet method and the midpoint Runge-Kutta method.  The results are shown in fig:moment, where it is clear that the Verlet method has dipole magnitude with bounded error as time goes on, while the midpoint Runge-Kutta method has error that grows linearly with time, similar to its treatment of the kinetic energy.  Overall, the magnetic velocity Verlet method is superior in preserving these constants of the motion for this simplest test case. 
  
  
  
  
  The dipole moment magnitude can be written in terms of the particle's kinetic energy in the direction perpendicular to the magnetic field.  For the uniform circular motion in this example, the perpendicular kinetic energy is
  
      K_⊥ = 1/2 m q B_0 R/m^2.
  
  Comparing this expression with the magnitude of the dipole moment in ucmdip, it is the combination K_⊥/B_0 that captures | m|.   
    
  
  In settings where the longitudinal magnetic field does not change magnitude much over some region, particle trajectories follow roughly circular motion in the perpendicular plane.  Then
  
      μ≡K_⊥/B,
  
  which is constant for uniform B_0, is approximately constant, and is called the “first adiabatic invariant." <cit.> 
  An equivalent expression for this constant, motivated by uniform circular motion where it follows from mudef, is
  
      μ = 1/2 q v_⊥ s
  
  for s, the radius of the circular motion with speed v_⊥.
   We can use the approximately constant value of μ to develop the mirror motion of a particle moving through a non-uniform longitudinal field.
  
  
  
  
  
  § CYLINDRICAL MAGNETIC MIRROR MACHINE
  
  
  
  Take a longitudinal magnetic field component pointing in the z direction with magnitude that changes as a function of z, B_z(z).  As a charged particle moves in the z direction, it is presented with a series of field magnitudes that are uniform in the xy plane, and will undergo 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 approximately circular motion occur before the longitudinal magnetic field changes significantly.  This is the “adiabatic assumption," that the particle's longitudinal motion is slower than its perpendicular motion.  We will check that this is satisfied, and also that μ is constant, in the numerical solutions below.
  
  
  A magnetic field cannot consist of B_z(z) ẑ alone; that violates ∇· B = 0.  In order to preserve azimuthal symmetry, introduce a magnetic field component pointing in the s direction that depends on s and z, B = B_z(z) ẑ + B_s(s,z) ŝ.  To make the field divergenceless, B_s(s,z) must be related to B_z(z) by
  
      B_s(s,z) = -s/2dB_z(z)/d z.
  
  The longitudinal force on the particle is no longer zero.  For the clockwise circulation associated with a positive charge, v_⊥ = -v_ϕϕ̂, and the equation of motion in the z direction is
  
      m z̈(t) = q v_ϕ B_s = -q v_ϕ s/2d B_z/d z.
  
  The term sitting out front in the second equality is precisely μ from mualt with v_⊥ = v_ϕ, so that 
  
      m z̈(t) = -μd B_z/d z.
  
  Multiplying both sides of this equation by ż(t), 
  
      m ż(t) z̈(t) =-μd B_z/d zd z(t)/d t,
  
  we can write both sides as total time derivatives provided μ is constant, 
  
      d/dt1/2 m ż(t)^2  = -d/dtμ B_z(z(t)).
  
  The integration is easy to carry out, giving
  
      1/2 m ż(t)^2 + μ B_z(z(t) ) = C
  
  where C is a constant of integration to be set by the initial conditions.  
  
  
  If a particle starts off at z(0) = z_0, with ż(0) = 0,  then C = μ B_z (z_0), and Kcons becomes
  
      1/2 m ż(t)^2 = -μ B_z(z(t)) - B_z(z_0).
  
  In order for the longitudinal speed to be real, B_z(z(t)) must be less than the initial value of B_z(z_0).  The particle will move towards regions of smaller field, increasing its longitudinal speed while circulating in the xy plane.  If the longitudinal magnetic field magnitude increases towards the value of B_z(z_0) at some location, the particle's longitudinal speed will decrease.  For a location z_1 with B_z(z_1) = B_z(z_0), the particle must have zero longitudinal speed.  The longitudinal kinetic energy in longKE plays a role similar to a one dimensional energy conservation equation in classical mechanics.  We can use it to predict some features of motion in the longitudinal direction, like turning points and points of maximum speed.  
  
  
  Using longKE, it is easy to see how to construct B_z(z) so as to get oscillatory longitudinal motion:  Make a field that is symmetric about z = 0 with maxima at ± z_0.  Then a particle that starts with no longitudinal speed at z_0 will also have zero longitudinal velocity component at -z_0.  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 plasmas without having mechanical pieces in contact with the plasma.<cit.>
  
  
  
  A simple, physically inspired way to achieve a magnetic mirror machine field is to put a pair 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 = μ_0 I a^2/2_≡ bẑ[ 1/ (d - z)^2 + a^2^3/2 + 1/ (d + z)^2 + a^2^3/2],
  
  with field magnitude shown in fig:Bzplot.
  
  
  For motion occurring near the z axis, this field is approximately valid, and we'll take it to define the longitudinal component of the target mirror field. The requirement in Bsdiv then provides the radial component.  Thus, the idealized magnetic field is <cit.>
  
      B   =  b [ 1/ (d - z)^2 + a^2^3/2 + 1/ (d + z)^2 + a^2^3/2]  ẑ
         - 3 b s/2[ d-z/a^2 + (d-z)^2^5/2 - d+z/a^2 + (d+z)^2^5/2] ŝ.
  
  
  
  Starting a particle off at z(0) = z_0 with ż(0) = 0, we expect it to move back and forth between ± z_0.
  We can take a look at one such trajectory to see how well the kinetic energy and the value of μ are conserved numerically.  For the  initial position and velocity, take
  
      r(0)    =    R x̂ + d/4ẑ
       v(0)    =    -q B_z(d/4) R/mŷ,
  
  where R = a/20 is chosen to be small compared to the radius of the current loops.  With these initial conditions, if the longitudinal field were uniform with constant value B_z(d/4), the particle would undergo circular motion of radius R with period
  
      T = 2 πm/q B_z(d/4).
  
   Running the magnetic Verlet method for a total time of 100 T in steps of Δ t = T/200, it is clear from the position plot in fig:d4data that the motion in z is periodic, going back and forth between ± d/4.  That period of oscillation is much larger than the period of the perpendicular circular motion,  with roughly five full cycles along the z axis occurring over the 100T time frame.  The longitudinal motion is slower than the circular motion, so the adiabatic assumption is satisfied here.  Another piece of our assumption was that the radius of the circular motion does not change much over the course of the trajectory, and that is true here, as shown in Fig. 4(b).
   
   
  
  
  
  
  
   In fig:d4em3, we've plotted the kinetic energy and value of μ from mudef using B_z ≈ B for the denominator since B_s is small.  In both of these plots, it is the ratio with the initial value 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 ∼ .005% change over the timescale shown.  The adiabatic constant μ is approximately conserved exhibiting larger ∼ .01% change over the 100 T time frame. 
    
   
    
  
  
  The initial velocity in velzero used a very specific tuning for its perpendicular component.  That velocity was taken to be the value that would ensure uniform circular motion of radius R in a uniform magnetic field.  But changing the initial value of v_⊥ just changes the radius of the circular motion. <cit.>  Starting from an initial position r(0) = (d/2) ẑ with velocity v(0) = -v_0 x̂, we can define the time-scale
  
      R    ≡   m v_0/q B_z(d/2), 
      
      T    ≡   2 π R/v_0.
  
  Then using a time step of Δ t = T/200, we made two trajectories with different values of v_0 = 2 π R/T.  The first, with v_0 =α≡ 2 π R_0/T for R_0 ≈ .06 m, has properties shown in fig:trajes2, and the second, with a value of v_0 = 3 α, has properties shown in fig:trajes3.  The larger speed leads to a larger radius for the circular motion of the particle shown in fig:trajes3.  In both cases, the kinetic energy is well conserved, while the first adiabatic invariant shows better conservation in the first case where the motion is closer to the z axis.
  
  
  
  
  
  
  
  
  
  
  § DIPOLE MAGNETIC MIRROR
  
  The spiraling, oscillatory motion we saw in the cylindrical setting persists in more complicated field geometries.  Imagine bending a magnetic field line from the previous section so that it curves.  We can still get particle motion that follows the curving field line, circling around it, while encountering regions of increasing magnetic field that act as mirrors.  Consider, for example, a dipolar magnetic field like the one outside the earth.  For a magnetic dipole pointing in the z direction, the field is 
  
      B = k/4 πx^2 + y^2 + z^2^5/2[ 3 x z x̂ + 3 y z ŷ - (x^2 + y^2 - 2 z^2) ẑ],
  
  where the constant k sets the magnitude.  We ran the magnetic Verlet method for a particle with initial position in the yz plane using a polar angle of θ_0 = 50^∘, with initial velocity in the x direction.  The result is shown in fig:dipole, where we can see one “cycle" of the mirrored trajectory, and then multiple cycles, going all the way around.
  
  
  
  
  
  It is interesting that the particle trajectory follows the field lines, spiraling tightly around them as it goes up and back.  This trajectory allows us to map out the dipole field visually using the motion of the charged particle.  For electric forcing, field lines point in the direction of acceleration of the particle, and if we start a particle from rest, charges travel along these lines.  We build intuition about electric fields using this property.  Magnetic field lines are perpendicular to the velocity vector, which points in the direction of motion, so we are not used to mapping magnetic field lines using the trajectories they generate.  The type of motion shown in fig:dipole, induced by the magnetic field of the earth, is what causes the Van Allen belts of particles trapped by the earth's field, traveling up and back along its field lines. <cit.>
   
  
  
  § MONOPOLE FIELD
  
  
  
  As our final example, we study the motion of a charged particle in the presence of a magnetic 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 a field line while corkscrewing around it.  As the particle encounters a region of increasing field magnitude, its motion along the field line slows, eventually stopping and reversing direction along the line.  This time, since the field lines converge radially, the geometry of the circulation perpendicular to a field line is conical.  The magnetic monopole case has the advantage that there is a closed form expression for the motion of the charge (see Griffiths' Problem 5.45 and references there <cit.>), giving us a rare exact result with which to compare the numerical trajectory.
  
  
  For a magnetic monopole, the field is B = k r̂/(4 π r^2) where k ≡μ_0 q_m is a constant that is set by the charge of the monopole.  As usual, the kinetic energy of a charged particle moving under the influence of this field is constant.  In addition, there is a conserved vector
  
      Q≡ m (  r× v) - k q/4πr̂,
  
  and this vector can be aligned so that it points along the z axis: Q = Q_0 ẑ.  Using the constancy of Q and the kinetic energy, one can show that θ(t) ≡θ is a constant of the motion, and develop expressions for the time derivatives of r(t) and ϕ(t):
  
      ṙ(t)    =    ±√( v^2 - (Q_0 sinθ/m r(t))^2 ),  
      ϕ̇(t)    =    Q_0/m r(t)^2,
  
  where v is the constant speed of the particle.  Dividing ṙ(t) by ϕ̇(t), the ODE governing the spherical r coordinate parametrized by ϕ is
  
      d r(ϕ)/d ϕ =±√(v^2 - (Q_0 sinθ/m r(ϕ))^2)(m r(ϕ)^2/Q_0).
  
  This equation can be solved, taking the minus sign, 
  
      r(ϕ) =- Q_0/m vsinθ/cos( (ϕ - α) sinθ),
  
  where α is the constant of integration.  We have lost the temporal evolution, but this equation can be used to draw a picture of the trajectory, and we can compare that with the solution to the equations of motion that we get numerically from the magnetic Verlet method.  The numerical solution requires a complete set of initial conditions.  We can get those by first choosing the constants of motion, θ (we took θ = 4^∘), α = -8 (chosen to give a trajectory that moved radially inward initially), and Q_0 (negative one, related to the size of the magnetic monopole).  For the initial value of ϕ, it is convenient to start at ϕ = 0, and then the initial value of r(ϕ = 0) is given by the solution in rofphi.  Finally, the initial values for ṙ(t = 0) and ϕ̇(t=0) can be obtained from dotr and dotphi respectively, using the rest of the initial coordinate values and constants.
  
  
  The numerical solution is shown plotted as points on top of the positions obtained from rofphi in fig:monopole.  There we can see the mirror effect as the particle moves down at first, then reverses direction and goes back up.   There is good agreement between the numerical solution and the exact one, with the two overlapping.  The constant Q is preserved by the numerical method with norm that varies by only (max(Q) -min(Q))/(min(Q)) ≈ 10^-7 over the portion of the trajectory shown. 
  
  
  Unlike the cases we have considered so far, the component of the motion that circles around a field line is not cylindrical here, even in approximation.  Instead, the charge moves around a cone with constant polar angle θ, with the tip of the cone at the monopole (shown on the right in fig:monopole), but this change from cylindrical to conical doesn't change the qualitative picture much.  We still have a particle that follows a field line while moving around it, slowing and reversing its direction along the line as it encounters increasing field strength.
  
  
  
  
  § CONCLUSION
  
  
  
  The extension of velocity Verlet to include magnetic fields provides a simple and useful tool for calculating the trajectories of charged particles moving in those fields.  While velocity Verlet 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 vector geometry associated with the Lorentz force in the context of Newton's second law.  One can easily extend the method to include additional position-dependent forces.  Going back to LorentzUp, an additional force F̅(x) would introduce a term like F̅( x(t)) + F̅( x(t + Δ t)), where both force evaluations rely on quantities that are known by the time the update is performed, so that these terms would get lumped into the vector d in wstart, at which point the derivation proceeds as in the text.  The method is fast, and can handle multiple particles (making it a favorite of molecular dynamics solvers <cit.>).
  We hope this paper serves to increase its use in the undergraduate E&M curriculum.
  
  
  
  
  
  
  
  The authors thank David Griffiths for useful commentary and for suggesting the monopole magnetic field example.
  
  
  
  18
  APMOT Jane M. Repko, Wayne W. Repko, Allan Saaf, “Charged particle trajectories in simple non-uniform magnetic fields," Am. J. Phys., 59, 652–655 (1991).
  GRIFFITHS David J. Griffiths, Introduction to Electrodynamics, 4th ed. (Cambridge University Press,  2017).  
  VANALLEN M. Kaan Öztürk, “Trajectories of charged particles trapped in Earth's magnetic field," Am. J. Phys., 80, 420–428 (2012).
  MAPLES George C. McGuire, “Using computer algebra to investigate the motion of an electric charge in magnetic and electric dipole fields," Am. J. Phys., 71, 809–812 (2003).
  COMPUTERSTUDIES Elisha R. Huggins, Jeffrey J. Lelek, “Motion of electrons in electric and magnetic fields; introductory laboratory and computer studies," Am. J. Phys., 47, 992–999 (1979).
  FRANKLINAJP J. Franklin, K. C. Newton, “Classical and quantum mechanical motion in magnetic fields," Am. J. Phys., 84(4), 263–269 (2016).
  MDVERLET Q. Spreiter, M. Walter, “Classical molecular dynamics simulation with the Velocity Verlet algorithm at strong external magnetic fields," J. Comp. Phys. 152, 1, 102–119 (1999).
  JACKSON John David Jackson,  Classical Electrodynamics, 3rd ed. (Wiley, 1998).
  GOLDSTEIN Herbert Goldstein, Charles P. Poole Jr. and John L. Safko, Classical Mechanics, 3rd ed.  (Pearson, 2001).
  VERLET Loup Verlet, “Computer `Experiments' on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules," Phys. Rev. 159(1), 98–103 (1967).
  SYMP Denis Donnelly, “Symplectic Integrators: An introduction," Am. J. Phys., 73 (10), 938– 945 (2005).
  AT M. P. Allen, D. J. Tildesley, Computer Simulations of Liquids, (Oxford University Press, New York, 1987).
  FRANKLINCMP Joel Franklin, Computational Methods for Physics, (Cambridge University Press, New York, 2013).
  POST  R. F. Post, “Summary of UCRL Pyrotron (Mirror Machine) Program," Proceedings of the Second International Conference on Peaceful Uses of Atomic Energy, Geneva, Paper 31, Vol 32, 245–265 (1958).  Available at: 
  
  
  
  
  
  
  INTNOTE  We encourage curious readers to find the source current density that produces this field and compare that source with the one that produced the field in Bzonly.  In addition, one can compare the field in mirrorfield with the full, off-axis, field produced by a pair of current loops placed at ± d along the z axis.
  PROVIDED There is also the possibility that the particle is not trapped at all.  One can develop, from the expression for μ, constraints on the ratio of the perpendicular to longitudinal velocity components that prevent trapping.  We leave this interesting opportunity, and its numerical verification, for the reader <cit.>.
  MD E. della Valle, P. Marracino, S. Setti, R. Cadossi, M. Liberti and F. Apollonio, “Magnetic molecular dynamics simulations with Velocity Verlet algorithm," 2017 XXXIInd General Assembly and Scientific Symposium of the International Union of Radio Science (URSI GASS), Montreal, QC, 2017, pp. 1-4, doi: 10.23919/URSIGASS.2017.8105168.
  
  