Abstract Field-reversed configuration (FRC) plasmas are of interest because of their potential application in magnetic confinement fusion reactors. Here we investigate the possibility of forming an FRC in a new way simply by allowing a hot region of plasma to expand against a background magnetic field in a conducting cavity. This differs from existing methods requiring input of large currents into the plasma by external sources. The expanding plasma case is studied via particle-in-cell computer simulations in both two and three dimensions. From these simulations it is found that at best the magnetic field in the center of the plasma can be reduced to fluctuations around zero as a result of diamagnetic currents at the plasma’s edge. A coherent FRC was never formed. Simulations with a finite volume ideal magnetohydrodynamics code are also run and compared to the particle-in-cell simulations. These results lead to the conclusion that this means of FRC formation is unlikely feasible. Particle-in-cell investigation of diamagnetic field reversal in expanding plasmas Jeffrey Kollasch INTRODUCTION Field-reversed configuration The field-reversed configuration (FRC) is a classification of magnetic field structure that occurs in plasmas and falls under the general category of compact tori (CT). The FRC magnetic field line topology is closely related to that of the streamlines in an incompressible flow structure called the Hill vortex which has a natural tendency to occur in situations such as skillfully blown smoke rings, unwanted recirculation in helicopter rotors, and the flow of blood into the human heart. The analogy to the incompressible Hill vortex is so strong that it’s magnetohydrodynamic (MHD) counterpart has been used by Zakharov and Shafronov to study axially symmetric current-carrying plasmas in equilibrium like the FRC [1]. The simple Hill vortex model of an FRC is plotted in Figure 1. This basic model of an FRC is cylindrically symmetric, has closed magnetic field lines inside an ellipsoidal region enclosed by a separatrix surface. Inside this region the field lines form closed loops pointing in a reversed direction in the center with respect to the external field. It is for this reason the configuration is called field-reversed. Two types of magnetic nulls exist in the FRC. The first is at the x-points which occur at both ends of the major axis of the ellipse. One knows the field is zero here because this is a condition for permitting field line crossing. The other magnetic null is at the o-point marked in Figure 1. The o-point extends to a ring in three dimensions and marks roughly the center of the area where a toroidal (azimuthal) current flows creating poloidal (transverse) field lines in accordance with Maxwell’s equations. The FRC shape is determined by zS, the distance from the center to an x-point, and rS, the radial distance from the center to the separatrix. The ratio of these is called the separatrix elongation; E = zS / rS. [Figure 1: Graph of FRC field-lines in r-z plane plane. Arrows indicate locations of magnetic nulls.] While FRCs can occur naturally in astrophysical contexts [2] the primary research motivation is their potential to supplement other magnetic confinement concepts for fusion energy such as the tokamak. There are currently over a dozen active experimental FRC research programs worldwide [3]. Potential engineering advantages to an FRC fusion reactor include lack of a center stack, simple cylindrical device design, possibility of a smaller less expensive machine, and potential to burn advanced fuels (i.e. those which produce fewer neutrons). Current FRC fusion concepts include the smashing of the FRC plasmoid in an imploding metal liner, the colliding of two plasmoids to produce shock heating, and the steady-state sustainment of a single FRC using odd-parity rotating magnetic fields (RMF) [3]. FRC proponents believe these advantages could provide a faster path to clean, reliable fusion energy. However, before this is possible understanding more thoroughly the processes by which FRCs are formed and sustained is of fundamental importance. Forming an FRC While the stabilization of the FRC plasma against instabilities such as the rotational, tearing, and tilt modes [3] is perhaps more intensely studied, in this work we focus on methods of formation. A number of formation methods already exist and are used to create laboratory FRCs. The standard method is called a θ-pinch. In a θ-pinch a linear device starts with field lines all pointing in one direction. Then an enormous discharge of azimuthal current is made in a cylindrical column around the device creating field in the opposite direction at the edges of the plasma. This provides the conditions necessary for the onset of magnetic reconnection forming the FRC separatrix with closed field lines inside. The plasma then contracts toward equilibrium. The θ-pinch is often used in magnetized target fusion (MTF) experiments. Other means of FRC formation include merging spheromaks, coaxial guns, and RMF current drive [3]. In this paper an alternative idea for creating an FRC is explored. Consider a vacated conducting cavity with a much smaller region of hot plasma released at the center. For simplicity let both the cavity and the plasma have cylindrical symmetry as in Figure 2. We now ask ourselves if the subsequent expansion of this out-of-equilibrium system might, given the right parameters, produce the currents necessary to reverse the background magnetic field in the centermost part of the plasma and cause an FRC to form. Although the statement of the initial conditions is simple, the ensuing dynamics are not. We therefore employ a particle-in-cell (PIC) numerical computer code to study this system. [Figure 2: Case where a plasma inside a conducting cavity expands with a background magnetic field producing diamagnetic currents.] METHODS AND MATERIALS The Lsp PIC code For this study our primary tool is Large-Scale Plasma, or Lsp, which is a leading commercial particle-in-cell (PIC) code designed to be used for simulating plasmas with time and length scales large enough to be relevant to fusion devices. Lsp is commercial software developed by Mission Research Corporation, ATK, and Voss Scientific. The particle-in-cell description of a plasma is particularly valuable to FRC research because the magnetic nulls make kinetic effects important which are difficult to capture with MHD or gyrokinetics models [4]. Both the basic principles of PIC simulation, and the particular features of merit possessed by Lsp will be briefly discussed. Lsp uses the particle-in-cell philosophy [5]. In a PIC simulation one models the plasma as a finite collection of particles moving against a background grid. In the case of Lsp this grid is a simple structured mesh based on Cartesian, cylindrical, or spherical coordinate systems. Even a small laboratory plasma may contain more than 10^15 particles so it is necessary to use a smaller number of particles in a PIC simulation and have each superparticle weighted in such a way as to represent a large number of ions or electrons. Another key approximation is storing the fields at grid points and then pushing particles based on an interpolation of these fields. This is the feature that makes PIC drastically faster than summing a full N^2 binary interactions between each particle and every other particle. These approximations are required to make the Lagrangian, particle-based method more tractable. It is also very common to use non-physical ion to electron mass ratios, however this approximation is not used in the simulations presented in this paper. With these approximations PIC can be a valuable tool for simulating plasma physics inadequately captured by fluid models. During a given time step in a PIC simulation the following occurs. 1) Given an initial distribution of electric and magnetic fields one calculates the Lorentz force for each particle by interpolating these fields off of the grid. 2) With the Lorentz force we now can numerically solve the ODE to advance particles’ positions and velocities. This step is often called the particle pusher. The algorithm can either be explicit (e.g. the very common leap-frog Boris mover) or implicit which allows for a relaxed temporal stability condition but requires costly matrix inversions. 3) The new electric current and charge densities are calculated from the updated positions and velocities of the superparticles. 4) Update the fields at the grid points according to Maxwell’s equations While numerous techniques exist to accomplish these tasks, the essence of PIC simulation is this simple. A full-featured code like Lsp might have over 200,000 lines of code; however a basic pedagogical PIC code like ES1 typically contains fewer than 2,000 lines [5]. One of the appealing features of Lsp is it has been written using extensible C programming in such a way that adding new physics and numerical packages is done very easily. Lsp contains models for things such as scattering (using the standard MCC method [5]), fusion, electrical circuits, ionization, and recombination. It also employs novel time advancing algorithms (Step 2 above) such as Energy Conserving Explicit (ECE), Direct Implicit (DI), and an presently under-development algorithm called Magnetic Implicit (MI). The parallel performance of these algorithms is vital for making simulations of FRC devices possible. Lsp can be run in serial or in parallel using the standard Message Passing Interface (MPI) for inter-processor communication. It has been benchmarked by the distributor using the QUICKSILVER PBENCH1 and PBENCH2 cases [6] with respectable 1 to 32 processor speedups of 26 and 20 respectively. Computational resources Lsp has been compiled and run on computing clusters at PPPL including the larger Dawson and Kestrel clusters, each having several hundred processors. Dawson will soon have over 800 processors. Lsp has also recently been set up to run on the Cray XT4 and XE6 architectures of the Franklin and Hopper national supercomputers. As of June 2011 these are the 8th and 27th fastest civilian supercomputers in the world respectively. A scaling study of Lsp was run on Franklin in which one of the 3D cases described below was run on 16, 32, 64, and 128 processors. The resulting speedups are less than linear as shown in Figure 3. This test case does not necessarily reflect the parallel behavior of Lsp on all problems, however it does indicate that good (i.e. near linear) scaling up to a more massive number of processors might pose a serious challenge. While some profiling was done, an attempt at speeding up the parallel implementation of the algorithms was not made. Many of the challenges in parallelization rest in the complexity of the physics models and relatively high number of communications required each time step. Also the long time step allowed by the implicit particle push means particles cross cell boundaries after a fewer number of time steps and the net result is a code which is communication-bound as opposed to compute-bound (as is the case with many explicit PIC codes without complex physics models). [Figure 3: Scaling analysis for Lsp on NERSC's Franklin XT4 machine.] RESULTS Numerous cases were run with Lsp in order to determine if the configuration described in Section 1 might produce an FRC. The first simulations were run in three dimensions in cylindrical coordinates in order to study physics more relevant to real plasma devices. While these simulations were enlightening they were also somewhat demanding in terms of computational resources, so a switch was made to two-dimensional simulations which revealed the essential physics results faster. The key result is that the magnetic field at the center of the expanding plasma cloud can be reduced to fluctuations around zero in some cases but never fully reversed to form an FRC. An ideal MHD simulation was also run to simulate a similar case with results that confirm those of the PIC simulation. Results from 3D PIC, 2D PIC, and 2D MHD will now be presented. Three-dimensional simulations Numerous three-dimensional simulations were performed with double precision arithmetic. Electrons were treated as an implicit species for the Direct Implicit (DI) particle advance algorithm. Ions were hydrogen nuclei with a full mass of 1836 times that of the electron. A typical case involved a conducting cylindrical cavity of dimension R=8 cm and Z=16 cm. The plasma is initialized hot with dimensions Rp=3 cm and Zp=8 cm. Figure 4 shows roughly to scale what the conducting domain and the initialized plasma particles look like at time t=0. [Figure 4: Initial state of typical 3D simulation at t=0.] The conducting wall is grounded. When particles hit the conducting wall they are simply removed from the simulation domain and therefore the system as a whole might have a net charge once particles begin to leave. We expect the particles to leave on a time scale inversely proportional to the ion acoustic speed. The ion acoustic speed is given by Cs = (γZkTe/mi)^(1/2) and has a value of 1.4E+7 cm/s for this particular hydrogen plasma. For the typical simulation we are presenting here the plasmas density is n=1.0E+12 cm^-3, the ion temperature is 1 eV, the electron temperature is 100 eV, and the initial BZ field is a constant that will be varied parametrically for the cases shown in this section. Simulation results are now presented. Figure 5 and Figure 6 illustrate the diamagnetic currents and the magnetic field due to the diamagnetic currents (not including background field) respectively. Figures 7 and 8 show contours of the electron density at times 160 ns and 320 ns respectively. Figures 9, 10, 11, and 12 show contours of the diamagnetic magnetic field (i.e. only that part of the field caused by the plasma currents and not including the background field) for cases with initial background fields of 0 G, 200 G, 300 G, and 500 G respectively. These particular cases have a 24x16x24 grid and 324,000 superparticles, with a time step of 1.8 ps for the DI algorithm. This time step is roughly the highest allowed for stability of the alternating direction implicit (ADI) field solver used by Lsp in conjunction with the DI algorithm. [Figure 5: Streamlines of the diamagnetic current produced in the plasma for the case where background B field was 500 G. They resemble a current loop.] [Figure 6: The component of the magnetic field due to the magnetic field corresponding to the current in Figure 5. This does not include the background magnetic field which is much stronger.] [Figure 7: Electron density profile at t=160 ns for the case where the initial background magnetic field was 0 G.] [Figure 8: Electron density profile at t=320 ns for the case where the initial background magnetic field was 0 G. The plasma is just now beginning to hit the wall and leave the system.] [Figure 9: The z component of the magnetic field due to currents in the plasma for the case where there was no initial background magnetic field.] [Figure 10: The z component of the magnetic field due to diamagnetic currents (not including the background field) for the case with an initial background field of 200 G. Only 13 percent of the field reduction needed for reversal was achieved.] [Figure 11: The z component of the magnetic field due to diamagnetic currents (not including the background field) for the case with an initial background field of 300 G. Only 6 percent of the field reduction needed for reversal was achieved.] [Figure 12: The z component of the magnetic field due to diamagnetic currents (not including the background field) for the case with an initial background field of 500 G. Only 16 percent of the field reduction needed for reversal was achieved.] Two-dimensional Cartesian simulations Since the simulations in three-dimensions were very time consuming it was decided it would be more direct to perform a higher number of less costly 2D simulations. A wide range of parameters were studied. Densities, ion temperatures, electron temperatures, magnetic field strengths, shape of the initial plasma, and the relative size of the initial plasma were all varied over at least one order of magnitude. Parameters affecting numerical accuracy such as grid size, number of particles, and time step were also varied. Many cases were possible as a typical simulation took between 1 and 30 minutes for the domain size described below. For the case which we will present shortly the plasma density was n=2.0E+12. The ion temperature and electron temperature were 1 eV and 100 eV respectively as before. We only show one case with a background B-field of 100 G in the positive x direction. The grid spacing is 200x100, a total of 400,000 particles are present at the start, and the time step is 10 ps. This time step is nominally in violation of the Courant stability condition based on the speed of light, however in practice it was found to be stable with the Energy Conserving Explicit (ECE) algorithm used for these cases. The wall is again conducting. The plasma parameters associated with the initial state of this case are given in Table 1. Table 1: Plasma parameters in standard CGS units for presented 2D Cartesian case at simulation start. - Ion acoustic speed (CS): 1.4000E+07 - electron cyclotron freq (fce): 2.8000E+08 - electron plasma freq (fpe): 7.9762E+10 - debye length (λd): 5.2538E-03 - thermal over magnetic pressures (β): 8.06E-01 - electron gyroradius (ρe): 0.238 - ion gyroradius (ρi): 43.705542 - electron inertial length (c/ω_pe): 3.75E-01 - ion inertial length (c/ω_pi): 1.61E+01 - Alfven speed (vA): 1.54E+07 The locations of a small fraction of the particles at initialization is shown in Figure 13 to show the system to scale. Figure 14 is a plot of the x component of the magnetic field measured at the very center of the domain over the course of the simulation. This is the key result and tells us that the initial 100 G magnetic field is reduced all the way to fluctuations around zero but does not continue past zero. Figures 15 shows the x component of the total magnetic field (now including the initial 100 G superposed) and Figure 16 shows the out of plane currents. [Figure 13: A scale representation of the system at t=0 showing a small fraction of the 400,000 total particles.] [Figure 14: A plot of the x component of the magnetic field measured at the center of the simulation domain.] [Figure 15: A contour plot of the x component of the magnetic field at time t=50 ns. The center are fluctuates around zero.] [Figure 16: The out of plane diamagnetic currents for the 2D case at t=50 ns.] [Figure 17: The magnetic field lines for the 2D case at time t=50 ns. The field is very near zero in the center and the field line integration algorithm is not accurate in this region.] Comparison to ideal MHD Related cases were run using an ideal MHD code called MHD2D (a code written by Peter Norgaard: http://www.princeton.edu/~norgaard/) which uses a Godunov finite volume method. Many cases were tried and it was similarly discovered that the magnetic field could never be reversed, but could be reduced to zero for cases where β is high. For an example consider a case with periodic boundary conditions and a domain that is 100 x 100 (these and all subsequent units are non-dimensionalized) in dimension with a 200x200 grid. We set the initial outward pointing magnetic field to 0.01. The pressure is 1.0 in a centrally located square of dimension 20x20 and 0.1 in the region outside (this bears similarity to the standard MHD blast wave benchmark). In the central plasma β=2P/B^2 = 2.0E+4. This is intentionally high. The time step was set to continuously adapt to a maximum CFL number of 0.75. We show a case in which the central plasma was given an initial counterclockwise rotation of ω=0.2. This is just one example of many cases, all of which never produce field reversal in the center even though some reduce the field to near zero. A graph of out of plane field is shown for time t=100 in Figure 18. [Figure 18: Ideal MHD case for which the out of plane magnetic field is reduced from 0.01 to near zero in the center at time t=100.] DISCUSSION AND CONCLUSION The fundamental result, as stated earlier, is that no simulations have been able to produce a coherent FRC strictly due to the diamagnetic currents caused when a plasma expands against a background magnetic field. However, under certain conditions, such as high β for example, we can completely reduce the field to near zero (MHD) fluctuations around zero (PIC) for a whole variety of parameters. It however does appear as though there is some force preventing the field to go below zero at the center and form an FRC. We are not surprised by this result in the case of ideal MHD since the frozen in condition on the magnetic field lines allows for exclusion of lines from the center but not for reversal in the center at least without considering numerical resistivity. The time scale for Spitzer resistivity is on the order of a microsecond for the systems we have presented while the events of the plasma evolution occur on a time scale roughly an order of magnitude smaller. We might expect other forms of resistivity to manifest themselves in PIC simulation however so this is not a complete explanation of why the field has not become reversed. It has been suggested that the near zero field in the center is a minimum energy state and going below zero would increase the energy of the system and therefore it does not naturally evolve in this way. Further consideration will be needed to fully understand this result. A logical next step is to incorporate initial currents in the plasma which might aid the diamagnetic currents in reversing the field. ACKNOWLEDGEMENTS Jeff Kollasch gratefully acknowledges the advising of Dr. Sam Cohen and Dr. Stephane Ethier at PPPL. Dr. Dale Welch at Voss Scientific also provided valuable assistance with Lsp. Dr. Elena Belova, Dr. Roscoe White, and Dr. Steve Jardin are also thanked for taking time to think about a possible explanation of the results. Joshua Blumenkopf, a fellow SULI student, is gratefully acknowledged for setting up a complementary 2D simulation with OOPIC to support the 2D simulations of Lsp. This work was supported by the US Department of Energy – SULI Program. REFERENCES [1] L. D. Zakharov and L. D. Shafranov, “Equilibrium of Current-Carrying Plasmas in Toroidal Configurations,” Reviews of Plasma Physics, vol. 11, p.153, 1986. [2] P. M. Bellan, “Spheromaks: a practical application of magnetohydrodynamic dynamos and plasma self-organization”, Icp, 2000. [3] L. C. Steinhauer, “Review of field-reversed configurations”, Physics of Plasmas, vol. 18, p. 070501, 2011. [4] D.R. Welch, S.A. Cohen, T.C.Genoni, and A.H. Glasser, “Formation of field-reversed-configuration plasma with punctuated-betatron-orbit electrons,” Phys. Rev. Lett., vol. 105, p. 015002, 2010. [5] C. K. Birdsall, A. B. Langdon, “Plasma physics via computer simulation,” Taylor and Francis, 2004. [6] M. F. Pasik et al., QUICKSILVER 2.4p Programmer’s Guide, Internal Document (available from M. Kiefer, Electromagnetics and Plasma Physics Analysis Dept.), 1998, Sandia National Laboratories.