Classical trajectory Monte Carlo simulations of particle confinement using dual levitated coils

The particle confinement properties of plasma confinement systems that employ dual levitated magnetic coils are investigated using classical trajectory Monte Carlo simulations. Two model systems are examined. In one, two identical current-carrying loops are coaxial and separated axially. In the second, two concentric and coplanar loops have different radii and carry equal currents. In both systems, a magnetic null circle is present between the current loops. Simulations are carried out for seven current loop separations for each system and at numerous values of magnetic field strength. Particle confinement is investigated at three locations between the loops at different distances from the magnetic null circle. Each simulated particle that did not escape the system exhibited one of four modes of confinement. Reduced results are given for both systems as the lowest magnetic field strength that exhibits complete confinement of all simulated particles for a particular loop separation. C © 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http://dx.doi.org/10.1063/1.4890305]


I. INTRODUCTION
A number of early plasma confinement approaches have been studied that make use of one or more current-carrying coils immersed in the confined plasma. 1-3However, work on approaches using immersed current-carrying coils slowed as tokamaks became preferred in fusion energy research (see discussion in Ref. 4).Some recent experiments make use of a single coil immersed in the confined plasma.Examples include the Levitated Dipole Experiment (LDX) 5 and Ring Trap 1 (RT-1), 6,7 both of which demonstrate the use of a levitated superconducting coil to magnetically confine the surrounding plasma.
Two plasma confinement systems are considered here.Each system consists of dual, levitated, coaxial coils carrying identical current.Diagrams of the plasma confinement systems studied are shown in Fig. 1, where the orientations of Cartesian coordinate systems are indicated.In one model, two identical current loops are coaxial with the z axis, with one loop centered at a small distance above and the other the same distance below the coordinate origin.In the second, the two loops have different radii, are concentric with the origin, and are coplanar in the z = 0 plane.The model consisting of axially separated loops will be referred to as the stacked system, and the model with concentric and coplanar loops will be referred to as the coplanar system.The term 'system' will be used to designate either the coplanar or stacked model, throughout.
The particle confinement capabilities of the stacked and coplanar systems are investigated using classical trajectory Monte Carlo (CTMC) simulations.To conduct the simulations, the magnetic field is calculated for each system using exact analytical expressions for a current-carrying loop and the superposition principle.The trajectories of a set of simulated particles starting from one of three initial positions is calculated by numerically solving the equations of motion subject to randomized initial velocity orientations.The simulations are carried out using a CTMC simulation program, for which a preliminary description is provided in Ref. 8. a Author to whom correspondence should be addressed.Electronic mail: cao@unt.edu. 1.The two dual-levitated-loop plasma confinement systems studied.Both systems consist of two levitated loops, each carrying identical current, I.The stacked system, top, consists of two current loops with identical radii, S, coaxial with the z-axis and centered at (0, 0, ±ε).In the second system, the loops are coplanar and concentric with the origin, with inner loop radius S and outer loop radius S e .
Plasma production systems incorporating dual current loops have been described for plasma processing applications. 9, 10Levitation of a dual superconducting ring system presents many challenges. 11Levitation of a system consisting of two concentric superconducting rings has been demonstrated. 12Also, plasma confinement using two physically supported coplanar dipoles carrying oppositely directed currents has been studied both theoretically and experimentally. 13n ideal plasma confinement approach for producing antihydrogen would be capable of providing long confinement times for a cold, dense, antimatter plasma.The systems studied here represent alternatives that may ultimately provide such confinement.Cold sources of antihydrogen are needed for conducting antimatter gravity studies, such as those described in Refs.14 and 15.
The ALPHA 16,17 and ATRAP 18 collaborations use nested Penning traps with superimposed minimum-B magnetic fields to confine antihyhdrogen.Using multipole fields to produce the magnetic minimum has been studied. 19,20 he dual current loop systems studied here possess null circles between the loops.The addition of a toroidal field with lines parallel to the null circle could be used to convert the dual current loop systems to finite-minimum-B traps.Such a trap might be capable of confining antihydrogen.
The ALPHA, ATRAP, ASAUSA, 21 AEgIS, 22,23 and GBAR 24 collaborations and the former ATHENA 25 collaboration have all aimed to address fundamental questions in current theoretical physics by studying antimatter.These questions include CPT (charge conjugation, parity, time reversal) invariance, the gravitational interaction of antimatter with matter, and the ratio of the inertial mass to the gravitational mass of antimatter.Experiments with antihydrogen seeking insight into these questions may require larger amounts of antihydrogen to be confined for a significant amount of time.
Currently, the nested Penning traps used for production, via three-body recombination, and trapping of antihydrogen exhibit significant plasma drift, which represents a formidable problem. 26he strong interest in antihydrogen research may lead to experiments that require alternatives to nested Penning traps for antihydrogen synthesis.These alternatives will need to confine cold, dense antiproton and positron plasmas in significant quantities over long time scales.
The work presented here extends knowledge of the dual-levitated-coil confinement approach by evaluating its ability to confine particles throughout a region of parameter space.In Sec.II, expressions for the magnetic fields of the two systems are developed.The equations of motion for the systems are introduced and put in dimensionless form in Sec.III.Simulation procedure and consistency control are described in Sec.IV, and results are presented in Sec.V.The modes of confinement are described in Sec.VI, and a final summary is provided in Sec.VII.

A. Field of a single loop
The magnetic field produced by a single current loop of radius S, in cylindrical coordinates (r, φ, z), lying in the z = 0 plane with the origin at the center of the loop, is (see, for example, Ref. 27) where the magnitude of the loop's magnetic field at the origin is and the coefficients are Here, m ≡ 4Sr (r +S) 2 +z 2 is a dimensionless parameter, K(m) is a complete elliptical integral of the first kind, E(m) is a complete elliptical integral of the second kind, I is the current in the loop, and μ 0 is the vacuum permeability.The complete elliptical integrals are The coefficients are used in the simulations and include an artificial parameter δ, which is given a sufficiently small value to remove the singularity at the origin without otherwise affecting the computations.

B. Fields of dual loop systems
Both systems consist of two current loops separated by a small distance.The superposition principle gives the total field of each system as the sum of the individual loop fields.In the stacked system, the loops carry identical current, have the same radius, and are centered at (x, y, z) = (0, 0, ±ε).Therefore, the expressions in Eq. (1) give the total field equations Here, the subscript s indicates that these equations apply only to the stacked system.In the coplanar system the loop radii are different and are labeled S for the smaller loop and S e for the larger loop.The field components for the inner loop, B ri (r, z) and B zi (r, z), are calculated by using S in Eq. (1).Similarly, using S e yields the expressions for the outer loop, B re (r, z) and B ze (r, z).Thus, the total magnetic field expressions for the coplanar system are B rc (r, z, S, S e ) = B ri (r, z) + B re (r, z), (6a) Here, the total field dependence on the loop radii is stated explicitly, and the subscript c indicates that the equations apply only to the coplanar system.

A. Normalizing parameters
The simulations are conducted within a normalized, dimensionless parameter space.A subscript n is attached to a parameter's symbol to indicate that it has been normalized.The quantities used for normalization are the radius of the current loops in the stacked system or the radius of the inner current loop for the coplanar system S, particle mass M, charge q, and kinetic energy K.Each of these parameters is normalized to unity, such that, Other parameters are normalized by writing them in terms of the un-normalized parameters.The separation between the loops is characterized by either ε n = ε/S or S en = S e /S, depending on the system.A particle's position, velocity, and acceleration are normalized as ) a, respectively.The normalized time and magnetic field are written as

B. Normalized magnetic field
The symbol r m is defined as the cyclotron radius of a particle moving in a circle in a uniform magnetic field.If the particle is moving in a field with the same magnitude as the field at the center of a single current loop given by Eq. (2), then the cyclotron radius and field magnitude are which, when normalized, is B mn = sign(q) √ 2/r mn .Here, sign(q) = q/|q|, and r mn = r m /S.Substituting Eq. ( 7) into the expressions for the magnetic field of a single loop Eq. (1), and writing the fields in terms of dimensionless parameters, gives the normalized, dimensionless magnetic field expressions for a single loop, The field expressions, Eq. ( 8), are used in Eqs. ( 5) and ( 6) with normalized parameters to obtain the normalized, total magnetic field expressions for the stacked system, and for the coplanar system.As a matter of preference, the simulations are conducted employing a Cartesian coordinate system.The total normalized magnetic field in Cartesian coordinates is Here, î, ĵ , and k are the unit vectors for the x, y, and z directions.The final form of the total magnetic field expressions used in the CTMC simulations is obtained by substituting Eqs. ( 9) or (10) into Eq.(11).The fields are depicted for both systems in Fig. 2.

C. Equations of motion
The CTMC simulations are conducted by solving the equations of motion for each particle subject to random initial conditions.The equations of motion are obtained by normalizing Newton's second law M a(t) = qv(t) × B, which yields Expanding the cross product gives the equations of motion in Cartesian coordinates, Here, the variables x n , y n , z n , and their derivatives are all functions of normalized time, t n .The simulated particles are mono-energetic with velocity Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions.Download to IP: 129.120.246.14On: Tue, 30 Aug 2016 21:51:37 2. The magnetic fields of the systems studied.A contour plot of the magnitude of the magnetic field is shown for the stacked system (a) and the coplanar system (c).Also shown are plots of the field lines produced by the stacked system (b) and the coplanar system (d).The locations of the saddle (P 1 ), interior (P 2 ), and null (P 3 ) points are labeled on each plot and are defined in Sec.IV.

Reuse of AIP
The initial position is specified by a selection of (x n0 , y n0 , z n0 ) for each particle.The initial velocity is calculated from the random numbers θ 0 and φ 0 using The magnitude of the normalized velocity is |v n0 | = √ 2 such that the normalized kinetic energy is Here, φ 0 is sampled uniformly on the interval 0 ≤ φ 0 < 2π and θ 0 is sampled on the interval 0 ≤ θ 0 ≤ π via the sampling expression θ 0 = 2 arcsin √ R θ .The sampling expression is obtained by inverting R θ = ( θ 0 0 sin θ dθ )/( π 0 sin θ dθ ), where R θ is a random number sampled uniformly between zero and one.Sampled in this manner, the velocity vectors generated are equally likely to take any orientation.

A. Simulation procedure
The trajectories are simulated by solving the equations of motion, Eq. ( 12), starting from a specified initial position, Eq. ( 13), with a randomized initial heading given by Eq. ( 14).Each trajectory is simulated until it reaches a boundary or a maximum time has elapsed.The two conditions are considered to be an escape or a confinement, respectively.The simulation repeats with the same parameter set and initial position, but a new randomized initial velocity for a specified number of particles.
The process repeats for seven values each of the separation between the two loops in the stacked system 2ε n and the normalized radius of the larger loop in the coplanar system S en .The magnetic field strength, characterized by the normalized cyclotron radius r mn , is varied incrementally for each value of ε n or S en to find the strength necessary for confinement of all simulated particles.Hereafter, a system with a specific value of ε n or S en is referred to as a 'configuration.' Confinement is tested at three initial particle positions.Graphic representations of their relative locations, labeled P 1 , P 2 , and P 3 , are shown in Fig. 2. The initial positions are defined such that they have analogous and unambiguous locations in all configurations.The initial positions are selected in the y = 0 plane and are specified by a choice of x and z coordinates.One initial position is located at a saddle point in the magnitude of the magnetic field and is labeled P 1 .
Another, labeled P 3 , is located along a line connecting the loops giving x n = 1 for the stacked system and z n = 0 for the coplanar.The second coordinate, z for the stacked system and x for the coplanar system, is the same as the corresponding coordinate of the saddle point in that configuration.In the stacked system, the initial position labeled P 3 is located halfway between the loops due to symmetry.In the coplanar system the saddle point in the magnetic field is shifted from being directly between the loops.Its coordinates were solved for numerically.The initial position labeled P 3 is less than 0.03 normalized distance units from the magnetic null circle in all configurations and hereafter is referred to as being (approximately) at the null point.
Also, the point halfway between the saddle and null points is used as an initial position and is labeled P 2 .Table I lists the coordinates of the three initial positions for each configuration.

B. Simulation duration and repetition
The particle trajectories are simulated starting at normalized time t n = 0 and continuing until each particle escapes or the maximum tracking time, denoted t n,max , elapses.To escape, a particle must reach the specified boundary with t n < t n,max .The boundary is defined as a finite-length cylinder coaxial with the z axis and centered at the coordinate origin.The cylinder has a radius of seven normalized units and extends twelve normalized units in the z dimension.The boundary is defined such that the ratio of the magnitude of the magnetic field at the boundary to the magnitude at a saddle point does not exceed 0.05 for any configuration.At the boundary, the calculated value of the cyclotron radius of a particle exceeds the particle's distance from the origin.
Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions.Download to IP: 129.120.246.14On: Tue, 30 Aug 2016 21:51:37 TABLE I.The coordinates of the three initial positions used, at a saddle point (P 1 ), near the null point (P 3 ), and halfway between the two (P 2 ), are listed.The coordinates are given as (x n0 , z n0 ), with y n0 = 0, for seven values of ε n for the stacked system and for seven values of S en for the coplanar system.An initial sample set of 500 trajectories is used with t n,max = 10 for all runs.Here and throughout, the term 'run' is used to designate a series of trajectory simulations for the same configuration using the same initial position, value of t n,max , and value of r mn but a different randomly oriented velocity vector for each simulation.If fewer than five trajectories reach the boundary in a run, the run is repeated with 1000 trajectories.Thus, each run exhibiting complete confinement consists of 1000 trajectories.
Each run repeats with the value of t n,max increased by a factor of 1.5 for each repetition, until the number of escape trajectories varies by no more than 10% of the sample size between two consecutive repetitions.To keep the computation time manageable, the maximum number of repetitions is limited to five.If the consistency condition is not met by the fifth repetition, the data from the final run is recorded and flagged.

V. RESULTS
Similar simulations were conducted for both systems.The fraction of trajectories which reach a system boundary with t n < t n,max was determined at each of the three initial positions.This fraction is referred to as the escape fraction and is denoted f esc .Magnetic field strength parameters in the range 0.375 ≤ r mn ≤ 5 were used for seven configurations of each system.In regions of parameter space where f esc 0.25, r mn was changed in increments of 0.0125.The size of the increment was increased for larger f esc values, and the largest increment used was 0.5.The escape fraction is shown versus r mn in Fig. 3 for trajectories starting at the saddle point P 1 with 2ε n = 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6 and with S en = 1.05, 1.1, 1.2, 1.3, 1.4, 1.5, and 1.6 for the stacked and coplanar systems, respectively.
Reduced results are given as the lowest magnetic field strength (highest r mn ) for which there were no escape trajectories.This value is referred to as the cutoff r mn and denoted R n .The R n values are listed in Table II.None of the data points at or very near the cutoff values were flagged for failure to converge.
The cutoff r mn is fit using least squares regression with a function of the form where α is either 2ε n or (S en − 1), depending on the system.These fit functions are plotted for both systems in Fig. 4. The fit coefficients are listed in  All distances listed in the tables are in normalized units.They are converted to unnormalized units through multiplication by the radius of the loops in the stacked system or the radius of the inner loop in the coplanar system S. Thus, ε = Sε n , S e = SS en , r m = Sr mn , and (x 0 , y 0 , z 0 ) = S(x n0 , y n0 , z n0 ).The r mn and R n values are related to the unnormalized magnitude of the magnetic field at the center of a single loop with radius S by B m = √ 2M K /(|q|Sr mn ).This expression is obtained from Eq. (7).With Eq. (2), the loop current necessary to produce the confinement condition, r mn ≤ R n , is If the particles follow a Maxwellian velocity distribution, an additional condition for single particle confinement is kT K, where T is the plasma temperature and k is the Boltzmann constant.With Eq. ( 7), the confinement condition, r mn ≤ R n , gives

A. Statistical uncertainty in simulated values
If each run is considered a series of Bernoulli trials, the probability of a particle escaping should be near the simulated fraction of escaping trajectories.The simulated values of f esc were found to increase by at least 0.25 between the runs with r mn = R n and r mn = R n + 0.5 for all configurations.(Recall that f esc = 0 when r mn = R n .)Thus, the particle's probability of escaping when r mn = R n + 0.0125 is greater than 0.006, assuming the probability increases linearly over this interval.(Recall that r mn was changed in increments of 0.0125 in the regions of parameter space where f esc 0.25.)The cumulative distribution function of the binomial distribution gives that the probability of at least one success with p ≥ 0.006 and 1000 trials is greater than or equal to 0.9976.The uncertainty is conservatively estimated as ±0.0125 for the R n values.
The normalized kinetic energy of a simulated particle K n was found to vary by less than 4 × 10 −6 during a simulation.Thus, numerical error is not expected to be significant.FIG. 4. The highest r mn value for which the fraction of escaping particles is zero versus loop separation.These values, which are given in Table II, are shown together with a least squares fit using a third-order polynomial.The symbol size is larger than the stated uncertainty in R n of ±0.0125.See Sec.VI A for details.
Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions.Download to IP: 129.120.246.14On: Tue, 30 Aug 2016 21:51:37 Here, α is 2ε n for the stacked system and (S en − 1) for the coplanar system.The fit coefficients are given for each of the three initial positions (P 1 , P 2 , and P 3 ) used for both systems.

B. Modes of confinement
Four modes of confinement are observed in trajectories for which r mn < R n .Trajectories can be confined near the null circle, can follow a field line encircling one or both loops, or can remain between the loops but away from the null circle.Simulated trajectories in the large radius limit of the two loop systems that demonstrate the modes of confinement are shown in Fig. 5.The large radius limit of both the coplanar and stacked systems consists of two infinite parallel wires carrying the same current, with a null line instead of a null circle.
Particles at the null line encounter an increasing magnetic field magnitude in every direction except along the null line.Most particles can only travel relatively short distances perpendicular to the null line before being reflected.Those that do leave the vicinity of the null line follow a closed magnetic field line and eventually return to the vicinity of the null line.The particle may become confined near the null line again, or it may follow a closed field line around the same or the other loop.These modes of confinement are demonstrated with two single-particle trajectories in Figs.5(a) and 5(b).
While traveling along a field line, away from the null line, a particle can be reflected by the increasing field magnitude near the loops and confined in the region between the loops, as demonstrated in Fig. 5(c).Some particles can travel along a magnetic field line that encircles both loops.This mode of confinement is shown in Fig. 5(d).
The confinement modes represented in Figs.5(a) and 5(c) indicate a magnetic mirror effect.Particles can switch between the modes represented in (a) and (b).However, with a sufficient mirror ratio the probability of this transition can become negligibly small.The mirror ratio is unbound approaching the null circle but is limited by the cyclotron radius of the trapped particle.Particles confined near the null circle can only be lost due to transport antiparallel to the gradient of the magnetic field.

C. Antihydrogen experiment
The stacked system is envisioned for conducting a gravity experiment, which is described in Ref. 14.The initial goal of the experiment would be to determine the direction of gravitational acceleration for antimatter.To date, no direct measurements have been reported that indicate the direction of free-fall gravitational acceleration of antimatter, such as antihydrogen.
The stacked system could be used for simultaneous confinement of antiprotons and positrons, under conditions whereby the two charged species may combine to form antihydrogen. Specific parameter values are given here by way of example: S = 5 cm, B m = 1 T, and 2ε n = 0.3.The value R n = 0.6 is obtained by taking the smallest of the three values reported in Table II.The right hand side of Eq. ( 17) has a value of 43 keV for antiprotons and a larger value for positrons.The confined antiprotons and positrons are both assumed to have temperatures of 4 K, which would satisfy the Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions.Download to IP: 129.120.246.14On: Tue, 30 Aug 2016 21:51:37 FIG. 5. Confinement modes of a single particle in the large-radius limit of the dual loop model.The limit is the same for both systems, and the locations of the wires are marked with +.The locations of the saddle and null points are marked with P 1 and P 3 , respectively, and the point halfway between them is marked with P 2 .The confinement modes are described in Sec.VI B. confinement condition indicated by Eq. (17).Antihydrogen atoms would be produced with a thermal distribution that is nearly identical to that of the antiprotons.
In Ref. 14, a setup was studied for carrying out an antihydrogen gravity experiment.Antihydrogen was considered to be produced from a 4 K source region for which radial streaming of the antihydrogen would be unimpeded.The stacked system, with the axis of symmetry oriented vertically, may provide such a source.The configuration would consist of the stacked system positioned at the geometrical center and two circular, parallel plates that have an axis of symmetry directed away from the center of the earth.The plates are separated by a small vertical distance, and include one or more pairs of circular barriers that protrude from the upper and lower plates, thereby forming an aperture between the plates.Antihydrogen annihilations that occur just beyond each barrier, within a "shadow" region, are detected.An asymmetry in the annihilation rate on the upper shadow region relative to the lower shadow region would occur when the effect of gravity is significant.
Antihydrogen atoms would not need to be trapped for conducting a gravity experiment similar to that described in Ref. 14. Simulations with an imposed toroidal field haven't been conducted in the work reported here.Nevertheless, it is speculated that adding a toroidal field and confining a sufficiently cold antihydrogen plasma may make it possible to use the stacked or coplanar system to simultaneously contain the two charged antimatter species together with a neutral antimatter species.

VII. SUMMARY
Two plasma confinement systems employing dual levitated coils were studied using classical trajectory Monte Carlo (CTMC) simulations.The systems were modeled as two levitated current loops.In one system, the loops are arranged coaxially and separated axially.For the second system, the loops have different radii and are coplanar and concentric with the coordinate origin.In both systems a magnetic null circle exists between the current loops, and closed magnetic field lines encircle both loops.
The particle confinement properties of the systems were investigated by simulating the trajectories of a large sample set of particles starting from three specified points in the system subject to a randomly-oriented, constant-magnitude initial velocity vector.The fraction of simulated particles that escaped the system by reaching a boundary within a given time interval was determined.The three initial positions used were at or near saddle and null points in the magnetic field and at a point halfway between the two.
The numerical results were summarized as the lowest magnetic field strength for which no simulated particles escape confinement.The results were given for seven different separations between the loops for each system and were fit with polynomials.
When confined, the simulated particles exhibited one of four modes.Trajectories starting near the null circle encountered an increasing magnetic field magnitude in every direction except along the null circle and tended to be confined near it.A fraction of simulated particles traveled more than a relatively short distance away from the null circle.As the magnetic field lines are all closed, these particles eventually returned to the vicinity of the null circle.Some trajectories starting near the saddle point in the magnetic field or halfway between the saddle point and the null circle followed a magnetic field line encircling both loops.Others were reflected by the increasing magnetic field magnitude near the loops.
FIG.1.The two dual-levitated-loop plasma confinement systems studied.Both systems consist of two levitated loops, each carrying identical current, I.The stacked system, top, consists of two current loops with identical radii, S, coaxial with the z-axis and centered at (0, 0, ±ε).In the second system, the loops are coplanar and concentric with the origin, with inner loop radius S and outer loop radius S e .

TABLE II .
The highest r mn value (weakest magnetic field strength) for which no particles escape (R n ) for trajectories starting at each of three initial positions for seven values of ε n for the stacked system and seven values of S en for the coplanar system.(Further details are in Sec.IV B.) The uncertainty in the R n values are estimated to be ±0.0125.(See Sec.VI A for details.)

TABLE III .
The cutoff r mn values from Table II are fit with the function R n