Drift velocity of charged particles in magnetic fields and its relation to the direction of the source current

Integrable motion of charged particles in magnetic fields produced by stationary current distributions is investigated. We treat motion in the magnetic field from an infinite flat current sheet, a Harris current sheath, an infinite rectilinear current, and a dipole in its equatorial plane. We find that positively charged particles as a rule will drift in the same direction as the current that is the source of the magnetic field in question. The conclusion is that charged particles moving under the influence of current distributions tend to enhance the current and that this indicates current self-amplification.


Introduction
Finding the motion of classical charged particles in a given magnetic field is a problem of great practical interest in many areas of physics and engineering [1][2][3]. One reason for investigating these problems is that they can throw some light on the much more difficult coupled problem of particles and fields in mutual interaction, a problem that can be approached in fluid mechanics using magnetohydrodynamics, in kinetic theory using the Vlasov equations, or using statistical mechanics based on the Darwin Hamiltonian. The latter approach, with its long range velocity dependent interactions, takes the energy lowering responsible for the attraction of parallel currents into account. Investigations by us [4][5][6][7][8][9] indicate that a plasma will generate currents and magnetic fields in its approach to equilibrium. Here we will study the drift of charged particles in the magnetic fields of given current distributions and conclude that this drift will, as a rule, enhance the current producing the field, thereby supporting the conclusions based on the Darwin Hamiltonian.
Drift velocities can be derived approximately for slowly varying fields by Alfvén's guiding center theory [1,10]. Analytical solutions for drift velocities published by Headland and Seymour [11] confirm the validity of the Alfvén formula The magnetic fields studied there are, however, difficult to realize in practice. Here we focus on analytical solution of drift velocities in magnetic fields that arise from idealized but reasonably realistic current distributions. Our a e-mail: hanno@mech.kth.se solutions start from the Lagrangian where m is mass, q electric charge, and A the vector potential.
In the examples we study the motion is periodic in one of the coordinates, say x. We denote this period with T . The drift velocity in another coordinate, say y, is then, i.e. the mean velocity over one period. This is the definition used below. We investigate the problem of the drift of a classical charged particle in a static inhomogeneous magnetic field produced by a given stationary current distribution. The following current distributions are treated: an infinite flat current sheet, an infinite current sheath found as a solution to the Vlasov equations by Harris [12], an infinite straight current carrying wire [13], current due to a rigidly rotating spherical charged shell i.e. the source of a dipole field.
The infinite flat current sheet case is mathematically trivial. It thus makes the physical reason for the drift especially clear, but we are not aware that this has been published before. Motion in the Harris sheath is integrable but the integrals require numerical solution (as far as we know). We are not aware that this problem has been studied before. Analytical results can be found for the limiting cases of slow and fast particles. The case of fast particles reduce to the flat sheet case.
In addition to Hertweck [13] the rectilinear current problem has been treated by Müller and Dietrich [14], Yafaev [15], Essén [16], and Aguierre et al. [17]. Here we focus on the drift velocity and how it is affected by the initial conditions. We point out that the maximum drift occurs for helical motion while zero angular momentum motion results in minimum drift. We believe that this is a new result.
The problem of particles in a magnetic dipole field is called the Störmer problem after Carl Störmer [18,19] who treated it in many publications from 1907 and onwards. There is an extensive literature (for a review see Dragt [20]) on this non-integrable problem, which is of interest as a model for the van Allen belts outside the Earth. When one restricts the motion to the equatorial plane perpendicular to the dipole the problem is integrable [21,22]. For certain initial conditions charged particles are bound between two radii and for these we calculate the angular drift velocity.

Charged particle near flat current sheet
In this section we calculate the motion of particles traversing a flat current sheet. We thus assume that there is a surface current density in the yz-plane k = σVŷ where σ is the surface charge density of current carriers and V their velocity in the y-direction. From the relation, wheren is normal to the surface of the current and B + is the magnetic field on the side thatn points to, while B − is the magnetic field on the opposite side, we find, using k = kŷ,n =x, B + = −Bẑ, and B − = Bẑ, that The vector potential is a limiting case of the potential in Figure 4 in Section 3, while the magnetic field can be written B = −B(x/|x|)ẑ.
In practice one can approximate such a sheet with many thin parallel current carrying wires. A small charged particle moving near the sheet can then be assumed to to pass essentially unhindered through such a sheet.
Since the magnetic field is constant on both sides of the sheet and simply changes sign in the sheet the trajectories of charged particles are easily calculated. Typical trajectories are illustrated in Figure 1. The circular trajectories of a charged particle have radii r determined by their speed v according to r = cm eB v. When the position of the circle is such that it is does not intersect the sheet the drift velocity is zero. For circular trajectories that intersect the sheet there will be a positive drift. For trajectories that are tangent to the sheet the drift velocity will be either v or −v depending on from which side the particle approaches the sheet.
Consider a particle initially moving with x > 0 in a circular trajectory with radius r and speed v centered at  a point with x-coordinate ξ. For 0 ≤ ξ < r we find that the drift velocity in the y-direction is from elementary calculations. If the center of the circle is on the other side of the x-axis −r ≤ ξ < 0 and one finds the drift velocity Since arccos(−x) = π − arccos x this expression is in fact identical to equation (7). The drift velocity in units of v as a function of ξ/r is shown in Figure 2. For ξ = 0 the center of the circle is on the sheet and the drift velocity is (2/π)v.

Charged particle in Harris sheath
The Harris sheath [12] is a plasma current density that is a solution of the collisionless Vlasov equations. The current is in the y-direction and is located near the yz-plane. The number density of charged particles n(r) is given by  (9)), of the Harris current sheath, with a peak at x = 0 (red curve). This charge density moves in the y-direction near the yz-plane. Also shown is the resulting magnetic field B(x), in the z-direction (Eq. (13)). The field is zero at x = 0 (green curve).
where V is the speed of the particles in the y-direction, n 0 the number density at x = 0, c the speed of light, and L D is the Debye length, with θ = k B T , Boltzmann's constant times the absolute temperature,while q is the charge of the particles, see Figure 3. The vector potential becomes where The magnetic field B corresponding to this vector potential is in the z-direction, B = B(x)ê z , and the z-component is We now study the dynamics of a charged particle in this magnetic field. The Lagrangian for a particle of charge q and mass m moving in the magnetic field is given by equation (2). We first divide the Lagrangian with the particle mass. We then choose 2θ mV = 1 (14) as our unit of speed and as unit of length. The unit of time is then mcL D /(2θ) = 1. The resulting Lagrangian is now One sees that the equation of motion in the z-direction (z = 0,ż = const.) decouples from the xy-motion and is trivial. We will thus discuss only the xy-motion in what follows and use the Lagrangian There are two constants of the motion, the generalized momentum, and the energy, so the problem is integrable, but the integrals are are not expressible in terms of well known functions, it seems.
To start with one notices that there is no magnetic field in the plane x = 0 (yz-plane). A particle in this plane withẋ = 0 will therefore stay in this plane and move as a free particle. The sign ofẏ depends only on the initial condition.
Note also that the positive sign in front of ln(cosh(x)) corresponds to a positive charge. A graph of this function is shown in Figure 4. It is clear that for positive κv the particle will always have a positive velocity in the y-direction.
For negative κv the y-velocity is negative for values of |x| < arccosh(e −κv ) but further away from the yz-plane it will become positive. We now investigate the x-motion. We see from (19) that plays the role of an effective potential energy for the x-motion. Graphs of this potential for different values of κv are shown in Figure 5. We therefore have from (23) From this one obtains The period T of the x-motion is then given by For the y motion is we have from (20) that dy/dt = κv + ln(cosh(x)) so that or, using (24), that The top curve (red) is for large speed, v 1, and is given by equation (40). The lowest curve (yellow) is for small speed, v 1, and is given by equations (33) and (34). For intermediate speeds the drift velocities are between these two. The middle curve (blue) is for v = 1 and is calculated numerically by doing the integrals (25) and (28). As a decreasing κ-value reaches −1 the trajectories bifurcate into two separate trajectories outside the sheath.
The change in y-coordinate during one period of the xmotion is thus (28) The y-drift velocity can now be calculated as function of v and κ as v y = Δy/T but the integrals must be done numerically (to the best extent of our knowledge) after finding the relevant turning points x − and x + . This drift velocity as a function of κ for v = 1 is shown in Figure 6.

Approximations for small and large v
By introducing the scaled variablex = x/ √ v one can go to the limit of small speeds in the integrals (25) and (28) and do them analytically. This then gives explicit expression for the drift velocity.
When v 1 it turns out to be advantageous to use the scaled variablex = x/v.

The nature of the sheath drift
In summary we have found that a positively charged particle will drift in the positive y-direction for all positive values of κ and also for most negative values of κ > −1.
Note that for κ < 0 the function V eff (x) develops a bump at x = 0 (see Fig. 5) and this means that particles spend more time inside the current sheath. For v = 1 drift in the negative y-direction starts at κ ≈ −0.80747. For κ < −1 the drift is always in the negative y-direction and the motion is bound in a potential minimum on one side of the sheath. The particles that drift in the positive y-direction are those that, not too slowly, pass through the sheath. One notes that the guiding center formula (1) correctly predicts the drift in the negative y-direction for particles staying well outside the sheath i.e. in a direction opposite to the current in the sheath. In our following examples the corresponding drift will be found to be in the same direction as the current. The main reason for this difference is that for a current localized near a point or near a line the magnetic field strength will fall off away from the source, whereas for a current localized near a surface it will increase away from the source.

Particle outside line current
First some notation and kinematics. We will use cylindrical coordinates, ρ, ϕ, z defined through, x = ρ cos ϕ, y = ρ sin ϕ, z = z, in terms of cartesian coordinates.
The general equations of motion can be found from the Lagrangian (2). For an (infinite) line current I along the z-axis one can use, with corresponding magnetic field B(r) = 2I cρê ϕ . The Lagrangian, divided by particle mass, is then of dimension velocity, represents a characteristic velocity of the system. Note that if q and I are of opposite sign this is actually a negative quantity. We use it as unit for velocity and put v 0 = 1 below. Since the coordinates ϕ and z are cyclic one finds that the corresponding generalized momenta are constants of the motion, while the ρ-equation of motion is, In terms of the constants of the motion (44)-(45) the conserved (kinetic) energy (per mass) can be written (47) We now assume initial conditions The energy is then Using (47) we findρ which using (49) is easily rearranged tȯ Inserting ρ = 1 we see thatρ = 0, as assumed. This means that ρ = 1 is the inner turning point in the ρ-motion. If we denote the outer turning point by ρ 1 we find that the period of the ρ-motion is given by Use of equation (44) and (46) giveṡ Integrating this over one period of the ρ-motion gives The average (or drift) velocity in the z-direction is thus, and it is positive if v 0 of equation (43) is positive (q and I of the same sign), otherwise negative. The same result has been derived by Yafaev [15] and a similar one can be found in reference [14]. Below we find explicit expressions for this drift velocity in the two special cases of helical and plane motion.

Helical trajectories
The equation for a helix of radius ρ = R can be written (φ = ω = const.), where we use time t as parameter. Introducing the pitch angle θ through, for the tangent (velocity) vector of this helix. Helical trajectories must be solutions of the equations of motion (44)-(46) withρ =ρ = 0, and ρ = R = constant. These equations give and comparison with equation (59) then gives, for the pitch angle θ of the helical trajectory. This can be in units of v 0 of equation (43), so that a negative v 0 means that v z is in the opposite direction. One notes that for large v most of the velocity is in the z-direction.

Plane trajectories
We now consider plane trajectories with ϕ = constant anḋ ϕ = L z = 0. Equation (51) at the outer turning point ρ 1 then gives, for the case of positive v 0 , which we assume in what follows. Using this (51) gives Changing the integration variable to ρ we thus find We now change variables in the integrals and put so that ln ρ(ln ρ 1 − ln ρ) = k y(1 − y) while ρ = 1 corresponds to y = 0, and ρ = ρ 1 to y = 1. The integrals then become These integrals can be expressed in terms of modified Bessel functions I ν (z), see Olver et al. [23], as follows Now, returning to dimensional quantities, one finds that where v = 1 2 ln(ρ 1 /ρ 0 ) is the dimensionless speed, ρ 0 the inner turning point, and v 0 is given by equation (43). A graph of this function is shown in Figure 8 where it is compared to the helix result (64).

Results for all solutions
Numerical results for arbitrary solutions show that the ratio of drift velocity along the current to speed, v z /v, is maximal for the helical trajectories and minimal for the plane ones. This means that all solutions will lie between the two curves in Figure 8. We illustrate this in Figure 9 where results for solutions between the two extremes are shown. On the horizontal axis is the (constant) and on the vertical axis the (pitch) angle θ of equation (60) at the outer turning point ρ 1 . This angle is zero for the plane trajectories. It obeys v = cos θ/ sin 2 θ for the helical ones according to (62). This corresponds to the dashed curve at the top of the diagram. For the other trajectories the θ-value represents an initial condition. The numbers on the contours are the values of v z /v. The conclusion is that for any solution of the equations of motion v z /v → 1 for v → ∞, for positive charge, otherwise to −1.

Charged particle in magnetic dipole field
Here we study the non-relativistic motion of trapped particles in the equatorial plane of a magnetic dipole. We calculate the angular drift velocity exactly and reach the same conclusion as above: the drift of the trapped charged particles is in the same direction as the current producing the field. The vector potential of a dipole m = m zêz in the z-direction is, when restricted to the equatorial plane, using cylindrical coordinates. The magnetic field is B = −(m z /ρ 3 )ê z . The Lagrangian (2) is The angular momentum, is conserved. The energy is, as always in magnetostatic problems, just the kinetic energy and can be written For > 0 the function has a minimum at ρ = where f = 0 and a maximum at ρ = 2 where f = 1/(16 2 ) (see Fig. 10). If we put the turning points in the ρ-motion are the solutions of For λ > 4 there are three positive roots of this equation given by: and one negative (un-physical) root. From these one can show that, Since on a circle with ρ = constant, the magnetic field is constant, there must be circular solutions to the equations of motion, as long as the speed is adapted the strength of the field. These are easily found and turn out to have radius ρ c = 2 and λ c = 4 . This means that they are located at the maximum of the curve f (ρ) in Figure 10. They are thus highly unstable. Orbits with λ > 4 and with ρ-values between the two smaller roots will remain there and will have an angular drift velocity around the dipole, which we now calculate. Using (79) and dt = dρ/ρ, one finds that, is the period of the ρ-motion. Using E = m 2 (ρ 2 + ρ 2φ2 ) we also find . (87) Fig. 11. The angular drift velocity (88) for particles trapped by a dipole field in its equatorial plane. In the units used the angular velocity increases with k where k = 0 corresponds to zero velocity and k = 1 corresponds to the maximum trapped velocity.
By doing these integrals one can calculate, the angular drift velocity. After some hard work the integrals turn out to give and where K(k), E(k) and Π(−k, k) denote complete elliptic integrals, see Olver et al. [23], and where Algebra shows that λ = 2 (1/k + k). Also ρ 1 and ρ 2 can then be expressed as functions of k. Using |L z |/E as unit of time we can now plot the angular drift velocity Δϕ/T as a function of 0 < k < 1 (see Fig. 11). The function is undefined at the endpoints of the interval but at k = 0 the limit is 3 and at k = 1 it is 4. For the case λ > 4 the charged particles can either be trapped between the two ρ-values ρ 1 and ρ 2 , or restricted to the region ρ > ρ 3 (see Fig. 10). Should λ < 4 or < 0 there will not be any region with trapped particles. Note that > 0 means that L z and qm z have the same sign. This means that positive (negative) trapped particles move round the dipole in the same (opposite) direction to the circulating current producing the dipole.

Conclusions
We have calculated the drift velocity of charged particles in the magnetic field of three different current distributions exactly: current located near a two-dimensional plane, current located near a one-dimensional line, and current in the neighborhood of a (zero-dimensional) point. In the two latter cases we find that charges will drift in a direction that enhances the current producing the field. In the dipole case this is valid for particles that remain near the dipole. In the case of current in the neighborhood of a plane the situation is complicated by the fact that the field increases away from the current but also in this case faster particles will drift to enhance the source current.
In the Darwin statistical mechanical approach [4] currents are predicted to be correlated in thermal equilibrium. This however does not explain this behavior within the present model: charged test particles moving in timeindependent external magnetic fields. The fact that there are exceptions to the rule for slow particles in the Harris sheath shows that a simple general explanation probably does not exist. For slowly varying fields one can find an explanation in the derivation of the the Alfvén formula (1) as given by e.g. Jackson [24] by noting how the field arises from the current distribution. For the straight wire case one has that the current is parallel to the z-axis so j ẑ. The magnetic field is then B ẑ ×ρ. Finally the gradient of |B| is ∇|B| −ρ. Putting this into (1) gives v D (ẑ ×ρ) × (−ρ) =φ × (−ρ) =ẑ j. For the dipole equatorial plane case we similarly have j φ, B −ẑ, and ∇|B| −ρ, so v D (−ẑ) × (−ρ) =φ j. For a current sheet on the other hand the magnetic field is not slowly varying. The explanation for this case and for fast particles in a Harris sheath is instead found in Figure 1. For slow particles outside the Harris sheath we get for the two regions outside the sheath, j ŷ, B ∓ẑ, but now ∇|B| ±x, so the Alfvén drift formula gives v D (∓ẑ) × (±x) = −ŷ −j. This explains the negative drift for slow particles outside the sheath.
The study was motivated by the much more difficult problem of finding the behavior of systems where charged particles and fields move according to the coupled equations of Lorentz and Maxwell. This difficult problem is at the heart of plasma physics. It seems that understanding how current is induced in a plasma in the neighborhood of a given current distribution, due to the inhomogeneous magnetic field from the source, should improve our qualitative understanding of this problem.