Extended space and time correlations in strongly magnetized plasmas
EExtended space and time correlations in strongly magnetized plasmas
Keith R. Vidal and Scott D. Baalrud a) Department of Physics and Astronomy, University of Iowa, Iowa City, Iowa 52242,USA Department of Nuclear Engineering and Radiological Sciences, University of Michigan, Ann Arbor, MI 48109,USA (Dated: 9 February 2021)
Molecular dynamics simulations are used to show that strong magnetization significantly increases the spaceand time scales associated with interparticle correlations. The physical mechanism responsible is a channelingeffect whereby particles are confined to move along narrow cylinders with a width characterized by thegyroradius and a length characterized by the collision mean free path. The predominant interaction is 180 ◦ collisions at the ends of the collision cylinders, resulting in a long-range correlation parallel to the magneticfield. Its influence is demonstrated via the dependence of the velocity autocorrelation functions and self-diffusion coefficients on the domain size and run time in simulations of the one-component plasma. A verylarge number of particles, and therefore domain size, must be used to resolve the long-range correlations,suggesting that the number of charged particles in the collection must increase in order to constitute aplasma. Correspondingly, this effect significantly delays the time it takes to reach a diffusive regime, inwhich the mean square displacement of particles increases linearly in time. This result presents challengesfor connecting measurements in non-neutral and ultracold neutral plasma experiments, as well as moleculardynamics simulations, with fluid transport properties due to their finite size. I. INTRODUCTION
Plasmas are collections of charged particles largeenough to exhibit collective behavior. Correspondingly,the assumption of a scale separation between the short-range interactions responsible for collective motion at themicroscopic scale, and the longer-range macroscopic evo-lution is inherent in the development of plasma kinetictheory. This scale is commonly characterized by the De-bye length because shielding limits the range over whichparticles interact. Thus, the collection of charges mustbe much larger than the Debye length in order to satisfythe definition of a plasma. However, this expectation isbased on the assumption that the plasma is weakly mag-netized in the sense that gyromotion is negligible at themicroscopic correlation scale.In this work, we show that the microscopic correlationscale characterizing interactions between particles can in-crease substantially if the plasma is strongly magnetized.In this context, strong magnetization refers to plasmas inwhich the typical gyrofrequency of particles significantlyexceeds the plasma frequency: β = ω c /ω p (cid:29)
1, where ω c = qB/m and ω p = (cid:112) e n/(cid:15) o m . Equivalently, thiscorresponds to particle gyroradii that are smaller thanthe Debye length, indicating that magnetization influ-ences inter-particle collisions. In particular, our molec-ular dynamics (MD) simulations of the one-componentplasma (OCP) show that strong magnetization causesparticles to become confined to cylinders with a lengththat is characterized by a mean free path for 180 ◦ col-lisions, and a width that is characterized by the gyro-radius. The length of these cylinders increases as β in- a) Electronic mail: [email protected] creases. This long-range correlation dramatically extendsthe scale separating microscopic interactions from macro-scopic behavior when the magnetization strength is large.It has implications for connecting experimental measure-ments , as well as MD simulations , to fluid trans-port processes.Strongly magnetized plasmas arise in a number ofcontexts. Non-neutral plasmas, such as those used totrap antimatter, can be very strongly magnetized β (cid:29) Electrons in a number of other laboratory experi-ments, including ultracold neutral plasmas, dusty plas-mas, high-field inertial confinement fusion experi-ments, and even to some extent in conventional mag-netic confinement fusion experiments, reach regimesin which β >
1. In nature, strongly magnetized plas-mas are encountered in several astrophysical contexts,including the magnetospheres of some planets and exo-planets, and the atmosphere of neutron stars. Manyof these examples are influenced by moderate or strongCoulomb coupling simultaneously with strong magneti-zation. Coulomb coupling can be quantified using theCoulomb coupling parameter, Γ = ( q /a ) / (4 π(cid:15) o k B T ),which is the ratio of the Coulomb potential energy atthe average interparticle spacing, a = (3 / πn ) / , to theaverage kinetic energy. Understanding how the combinedeffects of strong magnetization and strong coupling influ-ence transport, and the system scale required to exhibita plasma regime, is important to developing models todescribe these systems.Here, we show how coupling and magnetization in-fluence the space and time requirements for reaching aplasma regime using MD simulations of the OCP. Thespace requirement was determined through convergencetests of the parallel, perpendicular, and transverse ve-locity autocorrelation functions ( Z (cid:107) , Z ⊥ and Z ∧ ), andself-diffusion coefficients ( D (cid:107) , D ⊥ and D ∧ ), with respect a r X i v : . [ phy s i c s . p l a s m - ph ] F e b to the size of the simulation domain at a fixed tempera-ture and density. In the regime of weak coupling (Γ (cid:28) β < λ D = (cid:112) (cid:15) o k B T / ( q n ). However,if the plasma is weak to moderately coupled (Γ (cid:46)
1) butstrongly magnetized ( β (cid:29) (cid:38) β < ∼ parti-cles due to computational expense). It also implies thatsome experiments, such as non-neutral and ultracoldneutral plasmas, may not be large enough to rep-resent a three-dimensional plasma regime. These experi-ments access a wide range of particle numbers 10 − ,as well as coupling and magnetization strengths, but en-counter the same situation as the simulations in that verylarge particle numbers are required to represent a plasmaat high β . For example, particle numbers on the orderof 10 , Γ values near (cid:39) . β values ( >
10) arecommon in non-neutral plasma experiments.
The extended scale of spatial correlations implies a cor-responding extension of their timescale. Since the chan-neling effect generated by strong magnetization causesparticles to travel much further in the parallel dimensionbefore colliding with another particle, the time betweenthese collisions increases as β increases. This is demon-strated through the relaxation time of the velocity auto-correlation functions and the mean square displacementof particles from their initial positions in each the paral-lel, perpendicular and transverse directions. By the Ein-stein relation, the fluid process of diffusion correspondsto the regime in which the mean square displacement in-creases linearly in time. The time it takes to reach thisregime is related to the time it takes particles to interact,which is several plasma periods ( ω − p ) in a weakly mag-netized plasma ( β (cid:46) whereby the meansquare displacement in the perpendicular direction in-creases at a rate faster than proportional to time; i.e.,as t α where α >
1. This was later observed experimen-tally, but thought to be associated with not reaching thetime asymptotic limit associated with the hydrodynamicregime . We also observe this behavior in our 3D sim-ulations when they are not run sufficiently long. How-ever, eventually, the mean square displacement reachesthe normal diffusive regime in which it increases lin-early in time. This suggests that true diffusion persistsin strongly magnetized plasmas in three-dimensions, butthe time at which the diffusive regime is reached is de-layed significantly. Experimentally, this has implicationsfor the length of time over which the plasma must evolve,or measurements must be taken, in order to probe fluidproperties. For example, ultracold neutral plasmas ex-pand at a rate determined by ion dynamics, which couldbecome comparable to the time it takes electrons to reacha fluid regime if they are sufficiently strongly magnetized.The remainder of this paper is organized as follows.Section II describes the MD simulations and the calcula-tion of diffusion coefficients. Section III summarizes re-sults for the correlations scales in unmagnetized plasmasas coupling strength varies. Section IV shows how thesespace and time scales dramatically extend in response toan applied magnetic field. Conclusions are discussed inSec. V. II. SIMULATION SETUP
Simulations were prepared by randomly assigning po-sitions and velocities to N charged particles in a cubicbox with periodic boundary conditions. Velocities andpositions were updated by using the Velocity-Verlet al-gorithm. In simulations with an applied magnetic fielda modified Velocity-Verlet algorithm provided by Spre-iter and Walter that allows for an arbitrarily stronguniform external magnetic field was used. The simula-tion was performed in two stages. In the equilibriumstage the system was brought into constant energy con-figuration using a velocity rescaling thermostat. Afterequilibrium was achieved, particles evolved in time in themicrocanonical ensemble during the main run stage. Insimulations with an external magnetic field, the field wasturned on only during the main run stage and was notpresent during the equilibrium stage in order to reducecomputational time. The presence of an external mag-netic field in a classical system does not change the equi-librium state via the Bohr-van-Leeuwen theorem, thusturning on the magnetic field during the main run stagedoes not change the energy or configuration of the equi-librium state.The Coulomb potential was implemented using anEwald Summation technique through the Particle-Particle Particle-Mesh (P M) algorithm. This splits theCoulomb potential into short-range (Particle-Particle)and long-range (Particle-Mesh) components. This rep-resentation of the potential allows for high resolution ofclose interactions and a computationally efficient long-range force which is computed on a mesh using three-dimensional fast Fourier transforms.Positions and velocities of a random subset of particleswere recorded during the main run stage of the simula-tion, and used to compute the velocity autocorrelationfunction in each of three vector directions Z (cid:107) ( t ) = (cid:104) v z ( t ) v z (0) (cid:105) (1a) Z ⊥ ( t ) = 12 (cid:20) (cid:104) v x ( t ) v x (0) (cid:105) + (cid:104) v y ( t ) v y (0) (cid:105) (cid:21) (1b) Z ∧ ( t ) = 12 (cid:20) (cid:104) v x ( t ) v y (0) (cid:105) − (cid:104) v y ( t ) v x (0) (cid:105) (cid:21) (1c)where (cid:104) . . . (cid:105) denotes an average over the ensemble of savedtrajectories and the coordinate orientation is r (cid:107) = (ˆ z · r )ˆ z (2a) r ⊥ = r − r (cid:107) (2b) r ∧ = ˆ z × r . (2c)Here, B = B ˆ z , and the symbols (cid:107) , ⊥ , and ∧ refer tothe parallel, perpendicular, and transverse directions. Inthe absence of a magnetic field, Z ( t ) = (cid:104) v ( t ) · v (0) (cid:105) = Z (cid:107) ( t ) + Z ⊥ ( t ) and Z ∧ ( t ) = 0. Diffusion coefficientswere computed from the velocity autocorrelation func-tions using the Green-Kubo relations D (cid:107) = k B Tm (cid:90) ∞ Z (cid:107) ( t ) dt (3a) D ⊥ = k B Tm (cid:90) ∞ Z ⊥ ( t ) dt (3b) D ∧ = k B Tm (cid:90) ∞ Z ∧ ( t ) dt. (3c)In the absence of a magnetic field D = k B Tm (cid:82) ∞ Z ( t ) dt = D (cid:107) + D ⊥ , and D ∧ = 0. Since the simulations can onlybe run for a finite time, the upper limit in Eq. (3) is thesimulation run time t run , which must be sufficiently longfor the integral to converge. Determining the length oftime necessary for this to occur is a focus of this study.In establishing the run time necessary to reach a diffu-sive regime, it is also informative to consider Einstein’sdescription relating the diffusion coefficients to the meansquare displacement of a particle from its initial position D (cid:107) = lim t →∞ (cid:104)| z ( t ) − z (0) | (cid:105) t (4a) D ⊥ = lim t →∞ (cid:104)| x ( t ) − x (0) | (cid:105) + (cid:104)| y ( t ) − y (0) | (cid:105) t (4b) D ∧ = lim t →∞ (cid:104)| x ( t ) − y (0) | (cid:105) − (cid:104)| y ( t ) − x (0) | (cid:105) t . (4c)In the unmagnetized limit D = lim t →∞ (cid:104)| r ( t ) − r (0) | (cid:105) / (6 t ). In practice, the mean square displace-ment was computed not from the particle positions, but TABLE I. Unmagnetized OCP Simulation ParametersΓ δt ( ω − p ) t run ( ω − p ) ≥ . .
01 1638 . . − .
45 0 . . . − .
225 10 − . . − .
05 10 − − − . − − . β t run ( ω − p )0 . . −
10 3276 . − .
601 0 . − . − . . − . . −
10 5242 . − . from their velocities using r ( t ) − r (0) = (cid:82) t v ( t (cid:48) ) dt (cid:48) , so (cid:104)| r ( t ) − r (0) | (cid:105) = 2 (cid:82) t dt (cid:48) (cid:82) t (cid:48) dt (cid:48)(cid:48) (cid:104) v ( t (cid:48) ) · v ( t (cid:48)(cid:48) ) (cid:105) , which is di-rectly related to the velocity autocorrelation function (cid:104)| r ( t ) − r (0) | (cid:105) = 6 t (cid:90) t (cid:18) − st (cid:19) Z ( s ) ds. (5)Each of the vector components in the magnetized casewere obtained from using the appropriate component ofthe velocity autocorrelation function in this expression.The time-step, δt , and simulation duration, t run , usedfor simulations of the unmagnetized OCP are shown inTable I. The time step was chosen in order to conserveenergy to (cid:46) − . Nominally, the simulation durationmust be much greater than the decay time of the ve-locity autocorrelation functions, and converge at a ratefaster than t ( − / . This condition was satisfied forsimulations with Γ ≥ .
05, but for those with Γ < . β ≤ β > δt = 0 . ω − c =0 . ω − p /β in order to resolve the gyromotion. The rangeof run times indicated in the table are a result of the timeit takes to reach convergence in the cumulative integralof the velocity autocorrelation functions, and were foundto depend significantly on β . III. UNMAGNETIZED PLASMA
First, we review the correlation length and timescalesfor unmagnetized plasmas in order to set the expectations
FIG. 1. Percent change of the computed diffusion coefficientversus the number of particles used in the simulation. Per-cent change was calculated in reference to the result from thesimulation that used the largest number of particles. For (a),this was from the relaxation time model (“R.T. model,” redcircles), and for (b)-(f) from the cumulative integrals (“C.I.”,black circles). for comparison to strongly magnetized plasmas. Figure 1shows a convergence test of the computed diffusion coef-ficients with respect to the simulation size, quantified interms of the number of particles. Recall that since thedensity is fixed, the domain volume is proportional to theparticle number, and the domain is cubic: L = N/n .The expectation is that the simulation domain mustsignificantly exceed the range over which particles inter-act. At strongly correlated conditions (Γ (cid:38) (cid:46) N ≈ nλ D ∝ Γ − / . This is consistentwith the data shown in Fig. 1 where approximately 10 FIG. 2. Diffusion coefficients calculated from MD simulationacross a range of coupling strengths, Γ. Dots indicate MDsimulations that used 100 particles, and red circles 50,000 par-ticles. The dash-dotted line represents predictions Landau-Spitzer theory, and the dashed line Effective PotentialTheory (EPT). particles were required for convergence at Γ = 0 .
1, butmore than 10 appear to be required to achieve betterthan ∼
15% accuracy at Γ = 0 . δt (cid:46) r L /v T ∝ Γ / /ω p , where r L = e /k B T is the thermal distanceof closest approach, or Landau length. Simultaneously,the Coulomb collision relaxation time increases so thatthe simulation time requirement also scales steeply withdecreasing coupling strength: τ ∝ Γ − / ω p . Puttingthese together, the number of timesteps required scalesas τ /δt ∼ Γ − . The combination of the large particlenumber requirement and the long simulation durationrequirement lead to the conclusion that MD simulationsquickly become infeasible for Γ (cid:28)
1. Our simulationswere not able to reach convergence based on the cumu-lative integral method for Γ (cid:46) . ∼
20% ofthe expected result from standard theory. Figure 2 showsa surprising result that this trend extends to very weakcoupling, as convergence to the tens of percent level wereobtained with only 100 particles down to Γ = 10 − . Here,the data is compared with the accepted theoretical pre-dictions resulting from the lowest order Chapman-Enskog FIG. 3. The normalized velocity autocorrelation function, Z ( t ) /Z (0), and the cumulative integral, D ( t ), obtained from MDsimulations with varied number of particles, N , for Γ = 0 .
01 ((a) and (e)), Γ = 0 . Z ( t ) data and is only present in theΓ = 0 .
01, Γ = 0 .
1, and Γ = 1 plots. solution to the Boltzmann equation Dω p a = (cid:112) π/ / Ξ , (6)where Ξ → ln Λ = ln ( λ D /r L ) corresponds to the tradi-tional Landau-Spitzer formula applicable to the weaklycoupled limit. Here, the Effective Potential Theory(EPT) has been applied in order to extend the traditionaltheory to higher coupling (Γ (cid:46)
20) by calculating Ξ us-ing the potential of mean-force, but the predictions areequivalent for Γ (cid:46) . m d v ( t ) dt = − mξ v ( t ) + R ( t ) (7)with a constant friction coefficient, ξ , and an assumptionthat successive collisions are uncorrelated, (cid:104) R ( t ) · v (0) (cid:105) =0. With these assumptions, Eq. (1) reduces to a simpleexponential form Z ( t ) = Z (0) exp ( − tξ ) , (8)and the resulting self-diffusion coefficient is Dω p a = ω p ξ . (9)Figure 3 shows that the simulated velocity autocorrela-tion functions can be well fit by Eq. (8) with ξ as the fitparameter when Γ (cid:46)
1, and Fig. 2 shows that the dif-fusion coefficients resulting from using the fit ξ value inEq. (9) agree well with the theoretical predictions. Fig-ure 1 shows this relaxation time approximation model TABLE III. Summary of transport regimes. Here, λ col = denotes the collision mean-free-path, a = the average interparticlespacing, and r L = the thermal distance of closest approach (Landau length).Region Degree of Magnetization Gryoradius Criteria β − Γ Criteria1 Unmagnetized r c (cid:38) λ col β (cid:46) . / Ξ(Γ)2 Weakly Magnetized max { λ D , a/ √ } (cid:46) r c (cid:46) λ col . / Ξ(Γ) (cid:46) β (cid:46) min { , (cid:112) / (3Γ) } r L (cid:46) r c (cid:46) max { λ D , a/ √ } (cid:46) β (cid:46) / √ r c (cid:46) min { r L , a/ √ , λ col } max { / √ , (cid:112) / (3Γ) , . / Ξ } (cid:46) β agrees well with diffusion coefficients obtained from thecumulative integral when Γ (cid:46) IV. MAGNETIZED PLASMAS
Previous MD simulations of the magnetized OCP haveidentified four regimes in which different underlying phys-ical processes lead to qualitative changes in the diffusioncoefficients. These are summarized in Table III, andcan be understood based on a comparison of the gyrora-dius to other length scales relevant to collision physics:(1) Unmagnetized regime: When the gyroradius is largerthan the collision mean free path, particles do not gy-rate before scattering, so the magnetic field does notsignificantly influence collisional transport. (2) Weaklymagnetized regime: Here, the magnetic field influencestransport because the gyroradius is smaller than the col-lision mean free path, but it does not influence collisionsat the microscopic scale because the gyroradius is largerthan the Debye length. This is the traditional regimeof magnetized weakly coupled plasmas, as described byBraginskii transport theory. (3) Strongly magnetizedregime: Here, interactions between particles at the mi-croscopic scale are influenced by the magnetic field, sincethe gyroradius is smaller than the Debye length, but thegyroradius is much larger than the distance of closestapproach in a binary collision. Boltzmann kinetic the-ory does not apply in this regime, and generalized colli-sion models must be used. (4) Extremely magnetizedregime: When the gyromotion is smaller than any otherscale relevant to collision physics, particles move almost in one-dimension, but undergo nearly 180 ◦ degree colli-sions at the ends of collision cylinders.Previous work has explored the diffusion coefficient ineach of these regimes. Here, we focus primarily on thecorrelation scales associated with particle interactions inregions (3) and (4).
A. Correlation length
Strong magnetization causes a channeling effect inwhich particles move a long distance parallel to the mag-netic field within channels with a width characterized bythe gyroradius. This causes a long-range spatial correla-tion that can be observed by studying the convergence ofvelocity autocorrelation functions and self-diffusion coef-ficients with respect to the simulation domain size.Figure 4 shows computations of the velocity autocor-relation functions for Γ = 10 and β = 4 (region 4) as thenumber of particles in the simulation increases from 100to 500,000. Here, Z ⊥ and Z ∧ have been gyroaveraged inorder to visualize the behavior of the decay, since each ofthese components has a rapid oscillation at the gyrofre-quency that otherwise obscures the general trends; seeRef. 12. The figure shows that a very large particle num-ber is required for convergence in both the parallel andperpendicular directions, exceeding 10 particles. Thisis much larger than the approximately 500 particles re-quired to reach convergence in the unmagnetized regimeat the same Γ value; see Fig. 1. This slow convergence isalso clear in the cumulative integrals shown in the middlecolumn. The self-diffusion coefficients are computed fromthe plateau of these curves, demonstrating that the er-rors in the velocity autocorrelation functions add so thatthe error associated with insufficient particle number isstark. This error is exacerbated by the additional obser-vation that the true plateau associated with diffusion isdelayed until very late times, as shown in the right col-umn of the figure. Section IV B discusses how this delayin the approach to a diffusive timescale is related to theextended spatial correlation.Figure 5 shows the self-diffusion coefficients computedfrom the late-time plateau in the cumulative integral ofthe velocity autocorrelation function over a broad rangeof Coulomb coupling strengths (Γ = 0 . , ,
10 and 100)and magnetization strengths ( β = 10 − − FIG. 4. The velocity correlation functions ( Z (cid:107) ( t ) /Z (cid:107) (0), (cid:104) Z ⊥ ( t ) /Z ⊥ (0) (cid:105) g , and (cid:104) Z ∧ ( t ) (cid:105) g ) and cumulative integral ( D (cid:107) ( t ), (cid:104) D ⊥ ( t ) (cid:105) g , and (cid:104) D ∧ ( t ) (cid:105) g ) at Γ = 10 and β = 4 for the (a) parallel, (b) perpendicular, and (c) transverse directions. Theright panel plots show the cumulative integrals at an extended time scale. The (cid:104) ... (cid:105) g indicates the quantity was gyroaveraged.The dash-dotted line is a horizontal reference line. The arrow indicates the increasing N trend of the plots. identified in table III, which are indicated as the regionsbetween the vertical dotted lines in the figure. These re-sults confirm the expectation from Sec. III that only a fewhundred particles are required to obtain convergence inregions 1 and 2. However, as the magnetic field strengthincreases into regions 3 or 4, the D (cid:107) and D ⊥ coefficientsdepend significantly on the number of particles used andmore particles are required as β increases. This trend ismost pronounced in D (cid:107) at weak coupling, but it is alsosignificant in the D ⊥ component. For example, over 10 particles are required to reach convergence at Γ = 0 . β = 10. The spread of computed coefficients overthe range of N values narrows as Γ increases, becom-ing insensitive to the number of particles used over thisrange of β when Γ = 100. The D ∧ component is primar-ily associated with the Hall effect, which is insensitive tointeractions, and this is likely why this component is also insensitive to the number of particles used.The figure also shows a comparison with theory (redline). The unmagnetized and weakly magnetized regimes(1 and 2) are expected to be well-described by the tra-ditional Braginskii transport theory when Γ (cid:28)
1. Thelimitation to small coupling strengths is because theplasma kinetic theory upon which this solution is derivedapplies only to this regime. However, the recent effectivepotential theory (EPT) has provided a method to ex-tend plasma kinetic theory to the regime of Γ (cid:46) To lowest order in the Chapman-Enskog expansion, the
FIG. 5. Parallel ( D (cid:107) ), perpendicular ( D ⊥ ), and transverse ( D ∧ ) diffusion coefficients for (a) Γ = 0 .
1, (b) Γ = 1, (c) Γ = 10,and (d) Γ = 100. The red solid line shows the prediction of the EPT theory, and the open red squares, labeled BD, are resultsfrom Ref. 12. diffusion coefficients predicted by this model are D (cid:107) ω p a = (cid:112) π/ / Ξ (10a) D ⊥ = D (cid:107) (cid:18) π β Γ Ξ (cid:19) − (10b) D ∧ = √ πβ Γ / Ξ D ⊥ (10c)where the generalized Coulomb logarithm, Ξ, is calcu-lated from the potential of mean force as described inRef. 41. The results of this model are shown as the red line in Fig. 5.As expected, good agreement is observed in theregimes in which the theory is expected to apply: re-gions 1 and 2, and for Γ (cid:46)
20. For stronger magnetiza-tion strengths, in regions 3 and 4, the MD results divergefrom these predictions, with the D (cid:107) coefficient becom-ing dependent on β and the D ⊥ coefficient scaling moregradually with β than the β − dependence predicted bythe theory. These observations regarding D (cid:107) and D ⊥ have previously been made in Refs. 9 and 12. Figure 5shows the first MD computations of D ∧ over this rangeof parameters. The good agreement with theory even in FIG. 6. Percent change in diffusion coefficients in referenceto the values calculated from the simulation with the largestnumber of particles, N , for (a) Γ = 0 .
1, (b) Γ = 1, (c) Γ = 10,and (d) Γ = 100 at β = 10. regions 3 and 4 is likely because the D ∧ coefficient is de-termined entirely by the Hall effect in this region, andbecomes independent of the collision model.Next, we return to the question of why the MD simu-lations require a very large domain to reach convergencewhen the plasma is strongly magnetized, and weak tomoderately coupled. Figure 6 shows this trend from a dif-ferent perspective, plotting the self-diffusion coefficientsreferenced to the value computed from the simulationwith the largest particle number. This shows even moreclearly than Fig. 5 that an extremely large number ofparticles is required at strong magnetization when Γ ismoderate to small (Γ (cid:46) particles; extrapolatingthe observed trends suggests that it may take more than10 .The physical mechanism responsible for this long-rangecorrelation is revealed by looking at the particle trajec-tories. Figure 7 shows example trajectories from a sim-ulation at strongly magnetized conditions: Γ = 1 and β = 10. The most apparent feature is that particles movea long distance in the direction along the magnetic field,before scattering with a neighboring particle on a nearbyfield line, at which point it undergoes a 180 ◦ reflection.In this way, particles are largely confined to a “collisioncylinder” with a length of approximately the mean freepath for a 180 ◦ scattering event and a width of approxi- FIG. 7. Trajectories of a select few particles at β = 10, t run =10 ω − p , and N = 100 particles with B = B ˆ z for (a)-(b) Γ = 1,(c)-(d) Γ = 10, and (e)-(f) Γ = 100. The trajectories in the z - x plane ((b),(d), and (f)) change color when the particleunderwent a 180 ◦ collision. The dimension of the simulationbox was L = 7 . a . (This figure was created using OVITO ). mately a gyroradius. However, particles are also observedto continuously and gradually drift across the magneticfield. This aspect of the motion is associated with thecontinuous weak interaction between charged particlesthrough the long-range Coulomb force. The result is thatthe “collision cylinders” broaden, and even migrate, intime. Thus, both large-angle reflective scattering associ-ated with close interactions, and the small-angle continu-ous drifting associated with long-range interactions con-tribute to the overall transport. As the coupling strengthincreases, the plasma becomes more collisional and thedistance between 180 ◦ scattering events shrinks consider-ably. In the strongly coupled regime, particles are contin-uously undergoing strong inter-particle interactions andthese prevent particles from migrating far along the mag-netic field, even in a strong magnetic field regime. Thisreduces the spatial correlation length in the parallel di-mension, and a correspondingly smaller simulation do-0main is required to reach convergence.In the strongly magnetized regime, the most stringentrequirement for convergence appears to be that the par-allel dimension of the simulation domain be much largerthan the parallel correlation scale separating 180 ◦ scat-tering events. Since the boundaries are periodic, part ofthis requirement is that the same two particles do notinteract as a consequence of traveling through the do-main boundary, as this would cause an artificial corre-lation. However, determining the domain size require-ment is not as simple as estimating the intersection of acollision cylinder of gyroradius width intersecting near-est neighbor particles because the cross-field drift due tolong-range interactions causes collision cylinders to mi-grate.A test of this hypothesis is provided in Fig. 8. Thisshows the mean free path for 180 ◦ scattering events alongthe magnetic field for β = 10 and coupling strengthsΓ = 0 . , ,
10 and 100 as a function of the simulationdomain size. Data labeled “total” included all particlesin the mean free path calculation, whereas those labeled0 . . /
10 of the domain (dashed line). The notionthat the vast majority of collision mean free paths ( (cid:46) .
1) have to be much smaller than the simulation domain( L/
10) leads to an expectation that is consistent withthat observed in the convergence tests of Fig. 6.These results indicate that in a strongly magnetizedplasma, the correlation strength depends not only on theCoulomb coupling parameter (Γ), but also on the mag-netization parameter ( β ). Indeed, it has been pointedout previously that strong magnetization can cause bothpositive and negative correlations that are indicative ofa Γ > < Although the cor-relation scale of these interactions is long, they are alsoassociated with large-angle (180 ◦ ) scattering events, un-like the weak long-range interactions observed for Γ (cid:28) m d v ( t ) dt = − mξ · v ( t ) + q ( v × B ) + R ( t ) (11)where the constant friction coefficients are now repre- FIG. 8. The parallel mean-free-path for (a) Γ = 0 .
1, (b)Γ = 1, (c) Γ = 10, (d) Γ = 100. All simulations were at a β = 10 with t run = 1000 ω − p . The red dotted line shows thedimension of the the simulation box, L , and the red dashedline shows L/ sented as a tensor ξ = ξ ⊥ ξ ∧ − ξ ∧ ξ ⊥
00 0 ξ (cid:107) . (12)The associated parallel, perpendicular, and transversevelocity correlation functions obtained from Eq. (11) are Z (cid:107) ( t ) = Z (cid:107) (0) exp ( − tξ (cid:107) ) (13a) Z ⊥ ( t ) = Z ⊥ (0) exp ( − tξ ⊥ ) cos (cid:2) ( ω c − ξ ∧ ) t (cid:3) (13b) Z ∧ ( t ) = Z ⊥ (0) exp ( − tξ ⊥ ) sin (cid:2) ( ω c − ξ ∧ ) t (cid:3) (13c)and the resulting diffusion coefficients are D (cid:107) ω p a = ω p ξ (cid:107) (14a) D ⊥ ω p a = ω p (cid:20) ξ ⊥ ξ ⊥ + ( ω c − ξ ∧ ) (cid:21) (14b) D ∧ ω p a = ω p (cid:20) ω c − ξ ∧ ξ ⊥ + ( ω c − ξ ∧ ) (cid:21) . (14c)1 FIG. 9. The parallel [(a) and (d)], perpendicular [(b) and (e)], and transverse [(c) and (f)] velocity correlation functions, Z ( t ),and cumulative integrals, D ( t ), at Γ = 1, β = 3, and N = 5 × . The red line is the result of the fit to the MD data usingEqs. (13). The friction coefficients were obtained by fitting the sim-ulated velocity autocorrelation functions with Eq. (13)and the diffusion coefficients obtained from using the re-sults in Eq. (14). To reduce the accumulation of errors,the fit was applied to early times ( ≈ [0 , ω − p ).The relaxation time approximation assumes correla-tions between particles are weak. Figure 9 shows thatalthough this is an excellent approximation at small Γwhen β (cid:28)
1, it cannot capture the correlation causedby strong magnetization at the same Γ value. This canbe seen in both the parallel and perpendicular cumula-tive integrals in Fig. 9. It is particularly stark in theperpendicular direction, where the Z ⊥ predicted by therelaxation time approximation looks nearly indistinguish-able from the MD simulations. However, there is a slightphase-shift which, when integrated, leads to very largedifferences in the cumulative integral, and therefore, thepredicted diffusion coefficient. It is expected that theprimary cause of the failure of the relaxation time modelis the increased correlation scale introduced by strongmagnetization. However, it should also be noted that ithas recently been shown that the perpendicular frictionforce ( ξ ⊥ ) in magnetized plasma depends on the angle between the particle’s velocity and the magnetic field,and therefore cannot be characterized by a constant co-efficient even in the low-speed limit. This may couple Z (cid:107) and Z ⊥ , and is beyond the limitations of the modelpresented here. B. Correlation time
The simultaneous influence of large momentum trans-fer collisions at the ends of the collision cylinder, and thecontinuous drifting aspect of the motion, leads to a decayof the velocity autocorrelation function that is character-ized by two timescales. In a strongly magnetized plasma,the timescale associated with collisions from the ends ofthe collision cylinder can become very long; resulting ina longer time to reach the fluid (diffusive) regime. This isillustrated in Fig. 10, which shows the mean-square dis-placement in the parallel and perpendicular directions asa function of time over a wide range of magnetic fieldstrengths ( β = 0 . −
10) for Γ = 1 ,
10 and 100. In eachcase, there is an initial stage associated with the time forinteractions to be established within the collision volume.2
FIG. 10. Mean-squared-displacement (MSD) in the (cid:107) and ⊥ directions at several magnetic field strengths, β . The dashedand dash-dotted lines are linear fits to the β = 0 .
01 and 10data.
In a weakly coupled plasma, this timescale is character-ized by the plasma period because the collision diameteris the Debye length: τ ≈ λ D /v T = ω − p . In a stronglycoupled plasma, collective interactions determine the in-teraction range, which also occur near the plasma periodin the OCP. In each case shown in Fig. 10, there is aninitial slope that is steeper than linear in time associ-ated with this regime that lasts for the first few plasmaperiods.Qualitative differences between the weakly andstrongly magnetized regimes are observed at subsequenttimes, as the time it takes to reach a linear regime ismuch longer when β >
1, and the intermediate timescaleis steeper than linear in the perpendicular direction andshallower than linear in the parallel direction. The extentof this intermediate timescale extends with increasing β .Considering panel (b) of Fig. 10, the time to reach a lin-ear regime when Γ = 1 and β = 10 extends to morethan 1000 ω − p . Further evidence that this arrested de-cay is associated with the long-range parallel correlationis shown in panels (d) and (f), which illustrates that alinear regime is established sooner at higher Γ. As de-scribed in the previous section, and Fig. 7, the parallelcorrelation scale shrinks at high Γ due to the increasedcollisionality.These two distinct timescales can be observed directlyin the velocity autocorrelation functions and their cumu-lative integrals; see Fig. 4. In an unmagnetized plasma, the velocity autocorrelation function decays steadily overthe collision timescale so that a single plateau is observedin the cummulative integral at a time directly followingthe time at which the velocity autocorrelation functionhas visibly decayed to near zero; see Fig. 3. In contrast,two plateaus are observed at the strong magnetizationconditions shown in Fig. 4. The first is associated withthe time subsequent to the initial decay of the velocityautocorrelation function and may appear as though thisvalue should correspond to the self-diffusion coefficient.However, if the simulation is run much longer in time, asecond flatter plateau is observed much later, as shown inthe right column of Fig. 4. Comparing the cumulative in-tegral from Fig. 4 to the mean square displacement fromFig. 10(d) clearly shows that it is the late-time plateauthat is associated with the diffusive regime. The appar-ent plateau in the intermediate time interval is not quiteflat, though it can appear so to within the simulationresolution.These time profiles indicate that strong magnetizationarrests the approach to a fluid regime. If a simulation isnot run long enough, or a measurement is taken too earlyin time, this can leave the impression that the perpendic-ular diffusion is superdiffusive. Superdiffusion refers toevolution in which the mean square displacement growsfaster than linearly in time, such as in the interval from100-1000 ω − p in Fig. 10(b). However, we observe thatafter enough time has evolved to incorporate motion as-sociated with the long-range parallel correlation, that adiffusive (linear in time) regime is reached.A very wide range of timescales must be resolved inthe strongly magnetized regime, which presents a chal-lenge for MD simulations. Figure 11 shows the Fouriertransform of the velocity correlation functions. At zerofrequency the value of ˆ Z ( ω ) is equal to the diffusion coeffi-cient. At low magnetization ( β = 0 .
01) and low coupling(Γ = 1) the spectrum in each direction is characterized bya single broad peak at low frequency, associated with thedecay due to Coulomb collisions. The spectrum at lowmagnetization ( β = 0 .
01) and strong coupling (Γ = 100)has a similar broad low frequency component, but alsoan additional sharp peak associated with the plasmondispersion at ω = ω p , and a broader peak at ω ≈ . ω p that is associated with caging of particles by their nearestneighbors. In contrast, multiple peaks are observed spanning abroad frequency range as β increases; see Refs. 49–51 fora thorough discussion of the oscillation spectrum of themagnetized OCP in 3D, and Refs. 52 and 53 for 2D. Theˆ Z ⊥ and ˆ Z ∧ spectra show high-frequency peaks associ-ated with the upper hybrid mode: a longitudinal wavethat propagates perpendicular to the magnetic field at afrequency ω/ω p = (cid:112) β . This high-frequency peakis near the gyrofrequency and sets the timestep require-ment in the simulation. In the parallel direction, a promi-nent peak develops at a frequency a few times lower thanthe plasma frequency (near 0 . ω p for Γ = 1 and near0 . ω p for Γ = 100) as β increases. This peak is likely3 FIG. 11. The modulus of the Fourier transform, | ˆ Z ( w ) | , of the parallel, perpendicular, and transverse velocity correlationfunction. In the legend, β = 0 . β = 1 and 5 are exclusive to Γ = 100. due to the bouncing component of particle motion beingextended by the magnetic field (i.e., the long-range par-allel correlation), as shown in Fig. 7, which occurs moreslowly at Γ = 1 than at Γ = 100. In both cases a broadlow-frequency component also persists, which is likely as-sociated with the continuous drifting aspect from long-range Coulomb collisions. In the perpendicular directionat large coupling (Γ = 100) the correlation peak associ-ated with caging is observed to shift to lower frequency as β increases, as expected from the parallel correlation as-sociated with the formation of a collision cylinder (ratherthan sphere). V. CONCLUSIONS
Molecular dynamics simulations demonstrated thatstrong magnetization ( β (cid:29)
1) gives rise to a long-range correlation associated with interparticle interac-tions. The physical mechanism is a channeling effect, asa magnetic field of sufficient strength confines particlesto move largely in one dimension, scattering in nearly180 ◦ interactions with particles on nearby magnetic fieldlines. The resulting interaction volume is characterizedby an elongated collision cylinder, rather than the spher-ical geometry expected at weak magnetization. Theelongation of the cylinder is reduced when the Coulombcoupling strength increases due to an increased collisionrate. Influences of this correlation on the velocity auto-correlation function were observed to be reminiscent of what is observed at higher coupling strengths (Γ (cid:29) β (cid:29)
1) and ultracold plasmas are beginning to reach thisregime.
One motivation for these experiments isto perform sensitive tests of plasma transport proper-ties. Many of these experiments include on the orderof 10 − charged particles in total. Our work indi-cates that the number of particles required for these ex-periments to represent a plasma transport regime when β (cid:29) ≈ . β ≈ particles in or-der to represent a plasma. VI. DATA AVAILABILITY
The data that support the findings of this study areavailable from the corresponding author upon reasonablerequest.4
ACKNOWLEDGMENTS
The authors thank Dr. Jerome Daligault for supplyingthe MD code used in this work, for training in use of thecode, and for assistance in the data analysis and interpre-tation. This material is based upon work supported bythe Air Force Office of Scientific Research under awardnumber FA9550-16-1-0221 and by the National ScienceFoundation under award number PHY-1453736. It usethe Extreme Science and Engineering Discovery Environ-ment (XSEDE) , which is supported by NSF Grant No.ACI-1548562, under Project Award No. PHYS-150018. D. A. Gurnett and A. Bhattacharjee,
Introduction to PlasmaPhysics: With Space, Laboratory and Astrophysical Applica-tions , 2nd ed. (Cambridge University Press, 2017). B. R. Beck, J. Fajans, and J. H. Malmberg, Phys. Rev. Lett. ,317 (1992). B. R. Beck, J. Fajans, and J. H. Malmberg, Phys. Plasmas ,1250 (1996). D. H. E. Dubin and T. M. O’Neil, Rev. Mod. Phys. , 87 (1999). E. M. Hollmann, F. Anderegg, and C. F. Driscoll, Phys. Rev.Lett. , 4839 (1999). J. M. Kriesel and C. F. Driscoll, Phys. Rev. Lett. , 135003(2001). F. Anderegg, D. H. E. Dubin, M. Affolter, andC. F. Driscoll, Physics of Plasmas , 092118 (2017),https://doi.org/10.1063/1.4999350. M. Affolter, F. Anderegg, D. H. E. Dubin, and C. F. Driscoll,Phys. Rev. Lett. , 155001 (2016). T. Ott and M. Bonitz, Phys. Rev. Lett. (2011). T. Ott, M. Bonitz, and Z. Donk´o, Phys. Rev. E , 063105(2015). T. Ott, M. Bonitz, P. Hartmann, and Z. Donk´o, Phys. Rev. E , 013209 (2017). S. D. Baalrud and J. Daligault, Phys. Rev. E (2017). B. Scheiner and S. D. Baalrud, Phys. Rev. E , 063202 (2020). D. J. Bernstein, T. Lafleur, J. Daligault, and S. D. Baalrud,Phys. Rev. E , 041201 (2020). J. R. Danielson, D. H. E. Dubin, R. G. Greaves, and C. M.Surko, Rev. Mod. Phys. , 247 (2015). J. Roberts, J. Guthrie, and P. Jiang, in (APS, Virtual Meeting,2020). E. Thomas, R. L. Merlino, and M. Rosenberg, Plasma Phys.Control. Fusion , 124034 (2012). M. Bonitz, H. Kalhlert, T. Ott, and H. Laowen, Plasma SourcesSci. Technol. , 015007 (2012). B. Tadsen, F. Greiner, and A. Piel, Phys. Rev. E , 033203(2018). P. Hartmann, J. C. Reyes, E. G. Kostadinova, L. S. Matthews,T. W. Hyde, R. U. Masheyeva, K. N. Dzhumagulova, T. S.Ramazanov, T. Ott, H. K¨ahlert, M. Bonitz, I. Korolov, andZ. Donk´o, Phys. Rev. E , 013203 (2019). O. V. Gotchev, P. Y. Chang, J. P. Knauer, D. D. Meyerhofer,O. Polomarov, J. Frenje, C. K. Li, M. J.-E. Manuel, R. D. Pe-trasso, J. R. Rygg, F. H. S´eguin, and R. Betti, Phys. Rev. Lett. , 215004 (2009). M. Hohenberger, P.-Y. Chang, G. Fiksel, J. P. Knauer, R. Betti,F. J. Marshall, D. D. Meyerhofer, F. H. S´eguin, and R. D. Pe-trasso, Phys. Plasmas , 056306 (2012). M. R. Gomez, S. A. Slutz, A. B. Sefkow, D. B. Sinars, K. D.Hahn, S. B. Hansen, E. C. Harding, P. F. Knapp, P. F. Schmit,C. A. Jennings, T. J. Awe, M. Geissel, D. C. Rovang, G. A. Chan-dler, G. W. Cooper, M. E. Cuneo, A. J. Harvey-Thompson, M. C.Herrmann, M. H. Hess, O. Johns, D. C. Lamppa, M. R. Martin,R. D. McBride, K. J. Peterson, J. L. Porter, G. K. Robertson,G. A. Rochau, C. L. Ruiz, M. E. Savage, I. C. Smith, W. A.Stygar, and R. A. Vesey, Phys. Rev. Lett. , 155003 (2014). M. Shimada, D. Campbell, V. Mukhovatov, M. Fujiwara,N. Kirneva, K. Lackner, M. Nagami, V. Pustovitov, N. Uckan,J. Wesley, N. Asakura, A. Costley, A. Donn´e, E. Doyle, A. Fa-soli, C. Gormezano, Y. Gribov, O. Gruber, T. Hender, W. Houl-berg, S. Ide, Y. Kamada, A. Leonard, B. Lipschultz, A. Loarte,K. Miyamoto, V. Mukhovatov, T. Osborne, A. Polevoi, andA. Sips, Nucl. Fusion , S1 (2007). C. Paz-Soldan, N. W. Eidietis, R. Granetz, E. M. Hollmann,R. A. Moyer, J. C. Wesley, J. Zhang, M. E. Austin, N. A. Crocker,A. Wingen, and Y. Zhu, Phys. Plasmas , 022514 (2014). A. J. Creely, M. J. Greenwald, S. B. Ballinger, D. Brunner,J. Canik, J. Doody, T. F¨ul¨op, D. T. Garnier, R. Granetz, T. K.Gray, and et al., J. Plasma Phys. , 865860502 (2020). K. K. Khurana, M. G. Kivelson, V. M. Vasyliunas, N. Krupp,J. Woch, A. Lagg, B. H. Mauk, and W. S. Kurth, in
Jupiter: ThePlanet, Satellites and Magnetosphere , Vol. 1, edited by F. Bage-nal, T. E. Dowling, and W. B. McKinnon (Cambridge UniversityPress, New York, 2004) Chap. 24, pp. 593–616. A. K. Harding and D. Lai, Rep. Prog. Phys. , 2631 (2006). R. W. Hockney and J. W. Eastwood,
Computer Simulation UsingParticles (McGraw-Hill Inc., 1981). S. K. Tiwari and S. D. Baalrud, Physics of Plasmas , 013511(2018), https://doi.org/10.1063/1.5013320. G. M. Gorman, M. K. Warrens, S. J. Bradshaw, and T. C.Killian, arXiv (2020). I. L. Isaev and A. P. Gavriliuk, Journal of Physics B: Atomic,Molecular and Optical Physics , 025701 (2017). J. P. Hansen and I. R. McDonald,
Theory of Simple Liquids , 4thed. (Academic Press, Oxford, 2013). Y. Feng, J. Goree, B. Liu, T. P. Intrator, and M. S. Murillo,Phys. Rev. E (2014). D. Frenkel and B. Smit,
Understanding Molecular Simulation:From Algorithms to Applications (Academic Press, San Diego,CA, 2002). Q. Spreiter and M. Walter, J. Comput. Phys. , 102 (1999). R. Pathria and P. D. Beale,
Statistical Mechanics , 3rd ed., editedby R. Pathria and P. D. Beale (Academic Press, Boston, 2011). B. J. Alder and T. E. Wainwright, Phys. Rev. A , 18 (1970). L. Landau, Phys. Z. Sowjetunion , 154 (1936). L. Spitzer and R. H¨arm, Phys. Rev. , 977 (1953). S. D. Baalrud and J. Daligault, Phys. Plasmas , 055707 (2014). S. Chapman and T. G. Cowling,
The Mathematical Theory ofNon-Uniform Gases , 3rd ed. (Cambridge University Press, Cam-bridge, UK, 1991). S. I. Braginskii, Rev. Plasma Phys. vol. 1 (1965). D. H. E. Dubin, Physics of Plasmas , 052108 (2014),https://doi.org/10.1063/1.4876749. L. Jose and S. D. Baalrud, Phys. Plasmas , 112101 (2020). A. Stukowski, Model. Simul. Mater. Sci. Eng. (2010). T. Lafleur and S. D. Baalrud, Plasma Phys. Control. Fusion ,125004 (2019). Z. Donk´o, G. J. Kalman, and K. I. Golden, Phys. Rev. Lett. ,225001 (2002). T. Ott, H. K¨ahlert, A. Reynolds, and M. Bonitz, Phys. Rev.Lett. , 255002 (2012). T. Ott, D. A. Baiko, H. K¨ahlert, and M. Bonitz, Phys. Rev. E , 043102 (2013). H. K¨ahlert, T. Ott, A. Reynolds, G. J. Kalman, andM. Bonitz, Physics of Plasmas , 057301 (2013),https://doi.org/10.1063/1.4801522. P. Hartmann, Z. Donk´o, T. Ott, H. K¨ahlert, and M. Bonitz,Phys. Rev. Lett. , 155002 (2013). K. N. Dzhumagulova, R. U. Masheeva, T. S. Ramazanov, andZ. Donk´o, Phys. Rev. E , 033104 (2014). J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither,A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Pe-terson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr, Comput.Sci. Eng.16