Reduction of electron heating by magnetizing ultracold neutral plasma

Ultracold neutral plasma (UCP) experiments provide an excellent testbed for theories of transport properties of moderate to strongly coupled plasmas because they can be precisely diagnosed.¹ Such measurements are interesting from a fundamental physics perspective because they access exotic regimes of plasma physics. Validating theory also contributes to other areas of science, such as inertial confinement fusion^2,3 and astrophysics,^4,5 where moderate and strongly coupled plasmas are encountered in systems that are comparatively difficult to diagnose. In a typical UCP experiment, the coupling strength of electrons and ions are limited to $Γe ≲ 0.1$ and $Γi ≲ 5$⁠, respectively, due to disorder induced heating (DIH) and three-body recombination (TBR) heating processes.^6,7 Here

is the ratio of the average Coulomb potential energy of species “s” at the species interparticle spacing, $as=(3/4πns)1/3$⁠, to the average kinetic energy $kBTs$⁠, and Z_s and n_s are the charge (in electron units) and number density of species “s,” respectively. If either the electron coupling strength or the ion coupling strength, or both, could be increased, UCPs could reach physics regimes that are more interesting from a fundamental physics standpoint and which more directly resemble conditions in many of the other applications of interest.

In this paper, classical molecular dynamics (MD) simulations are used to show that a strong external magnetic field can reduce electron heating, leading to a larger electron coupling strength. The electron temperature is reduced because electron motion is restricted to nearly one-dimension, parallel to the direction of the applied magnetic field $B$⁠. In this paper, electron heating is studied for times up to 10 $ωpe−1$⁠, where $ωpe=e2ne/ϵ0me$ is the electron plasma frequency. Over this duration, neither plasma species acquires thermodynamic equilibrium. The notion of temperature is represented by the “kinetic temperature” defined as

Here, $Ts∥$ and $Ts⊥$ are the kinetic temperatures for species “s” in the parallel (x direction) and the perpendicular direction to the applied magnetic field, respectively. N_s and m_s are the total number of particles and the mass of single particle of species “s,” respectively. Heating occurs predominantly in the parallel direction, so the total kinetic temperature is observed to be approximately one-third of the value obtained in the absence of $B$⁠. A large temperature anisotropy is established as a result. Eventually, this anisotropy relaxes, and heating also occurs in the perpendicular directions. However, if the magnetic field is strong enough, the relaxation can be dramatically suppressed,^8–12 and may be delayed long enough that measurements could be made before heating occurs in the perpendicular directions. This may increase the effective Coulomb coupling strength, Γ, as measured in terms of the kinetic temperature.

An applied magnetic field also introduces interesting transport and collective properties, causing magnetized ultracold plasmas to differ from typical ultracold plasmas to a degree that depends on the strength of the magnetic field.^13,14 For instance, the dynamic structure of a magnetized plasma will include features associated with collective modes such as upper and lower hybrid and shear modes, along with the acoustic modes.¹⁵ The present work focuses only on the increased effective coupling strength of an ultracold plasma due to the applied magnetic field.

In experiments, creation of UCP is followed by DIH on a timescale characterized by $ωpe−1∼$ (nanoseconds)^16–18 and later by TBR heating.¹⁹ Ions are also influenced by DIH, but on a longer timescale characterized by the ion plasma period (⁠$ωpi−1$⁠).^20,21 Previous and ongoing efforts to increase coupling strength largely focus on ions, including laser cooling,^22,23 introducing initial spatial correlation using blockaded Rydberg gases⁷ or Fermi gases,²⁴ or by trapping the gas in an optical lattice.^25,26 We focus on electrons because they are more easily magnetized than ions.

Experiments have shown that a magnetic field of a modest magnitude (⁠$∼50$ G) can substantially reduce transverse expansion of an ultracold plasma.²⁷ Previous theoretical studies have also suggested that a strong magnetic field can substantially reduce the rate of TBR in a weakly coupled plasma, and therefore the heating associated with it.^28,29 The magnetization strength of a species can be quantified by the parameter $βs=ωcs/ωps$⁠, which is the ratio of the gyrofrequency (⁠$ωcs=qsB/ms$⁠) to the plasma frequency. Here, we focus on the regime where electrons are strongly magnetized $βe≫1$⁠. In this regime, electrons are constrained to small gyro-orbits, which reduces their collision rate significantly in the transverse direction. This suppresses DIH in the transverse direction. The weakest magnetic field at which any influence is observed is $βe≃1$⁠. At the density $n≃108$ cm⁻³ relevant to UCPs, $βe=1$ if B = 32 G. A much more substantial reduction in transverse DIH is observed for $βe≳10$ (corresponding to $B≳320$ G). Previously, Glinsky et al.,⁸ Dubin et al.,^10,11 and Ott et al.¹² have shown that a long lived temperature anisotropy can be maintained for $β>1$ in weakly, moderately, and strongly coupled plasmas, respectively. Recently, Baalrud and Daligault have provided the temperature anisotropy relaxation rates ranging from weak to strong coupling³⁰ and for magnetic field strength β ranging from 0 to 100 using MD simulations.³¹ Based on this data, it is expected that a magnetization strength of $βe≳100$ (corresponding to $B≳3200$ G) may be sufficient to delay the relaxation of the electron temperature anisotropy long enough (⁠$∼106ωpe−1$⁠) to make measurements over the duration of a typical experiment (⁠$∼100μ$s).^6,21

In the strong field regime, we also observe a reduced rate of heating due to TBR. Like DIH, the reduction is observed to be a factor of three in the strong field regime (⁠$βe≳20$⁠). Previous theoretical work has suggested that the rate of TBR may be suppressed by a factor of 10 in a strong magnetic field at weak coupling. Our work addresses the moderately coupled regime.

The paper is organized as follows: Section II provides details of the simulation model. Section III discusses the heating processes in an unmagnetized ultracold plasma, and this is extended to include strong magnetization in Sec. IV. Potential implications for experiments are discussed in Sec. V, and a summary is provided in Sec. VI.

Non-equilibrium classical MD simulations³² were carried out using the open source code LAMMPS.³³ The simulation model involved a cubic box with periodic boundary conditions in each of the three dimensions. Equal numbers of electrons and ions (10⁴ particles of each species) were introduced at random positions, so as to avoid any initial spatial pre-correlation amongst the charged particles that would have a suppressing effect on DIH.^7,26 Each simulation was carried out with three different initial random configurations, and the presented results are a mean of these independent runs. Each species was initialized with zero temperature (i.e., all particles had zero initial velocity). Simulations were also carried out with a very small but finite initial temperature (corresponding to $Γ≈200$⁠) and the results were confirmed to be the same as those with zero initial temperature. The results show a focus on plasma with singly charged strontium ions (⁠$mi=1.62×105me$⁠) to model a common experimental choice. Cases with different ion to electron mass ratio will be stated explicitly. At such a high mass ratio, ions were approximately stationary over the timescale of several $ωpe−1$ that is the focus of the present study. It was confirmed that similar results were obtained over this time interval if ions were kept stationary at their initial positions.

An electron-ion plasma was modeled taking like-charged particles to interact through the repulsive Coulomb potential and particles with opposite charge to interact through the attractive Coulomb potential along with a repulsive core

Here, r is the separation between two charged particles, e is the elementary charge, and α is an adjustable parameter that sets the e-i repulsion length scale. The repulsive core was included to avoid the Coulomb collapse between oppositely charged point particles. For the purpose of this work, α was chosen to be sufficiently small that it did not influence the results in any way. It is a purely numerical necessity to avoid the rare occurrence of the very close interactions between electrons and ions, which is necessary to prevent because arbitrarily close interactions cannot be resolved by a finite time step integration, leading to energy conservation issues.

Recently, Tiwari et al. used this model to calculate thermodynamic state variables in recombining ultracold plasmas.³⁴ It was found that the magnitude of α controls the concentration of free charges compared with classical Rydberg states at equilibrium because it determines the depth of the potential well formed between oppositely charged particles. “Bound states” refers to the classical analog of the true TBR recombination processes, which would include discrete quantum energy levels. By distinguishing the free and classically bound populations using both an analytic ansatz and by counting them directly from the particle trajectories, it was shown that the excess pressure and excess internal energy of the free charges becomes independent of the choice of repulsive core scale length when it is sufficiently short-ranged. This enabled simulations of thermodynamic properties of an ultracold plasma at fixed conditions. Here, this model is applied to study the dynamic properties on short timescales. This is a simpler situation because few bound states are formed on the short timescales of interest.

The system was evolved according to the equation of motion

where $j$ is the index of any particle in the system (⁠$j=1…N$⁠), $Fj=−∑i=1N−1∇vij$ is the total force on jth particle due to interaction with all other particles, and $B=Bex$ is the magnetic field. A modified Velocity-Verlet scheme (as described by Spreiter and Walter³⁵) was used to evolve the equation of motion.³⁶ The time step was chosen to resolve the gyromotion, electron plasma frequency, and the orbits of electron–ion interactions. For $β≤10$⁠, it was chosen to be $Δt=10−3ωpe−1$⁠, and for $β>10$ to be $Δt=10−2ωce−1$⁠. The total energy was conserved to within 0.01% for a time span of 10 $ωpe−1$⁠.

The ultracold plasma evolution has been modeled using classical molecular dynamics with a soft repulsive core potential by Kuzmin and O'Neil^16,37 with a hybrid molecular dynamics model by Pohl et al.³⁸ and recently by Isaev and Gavriliuk to study the heating and diffusion processes in the presence of a homogeneous magnetic field.³⁹ It will be shown below that the results are insensitive to the particular form of repulsive core potential.

In an unmagnetized ultracold plasma, electron-DIH saturates on a timescale of 1–2 $ωpe−1$⁠, when the kinetic energy of the electron species becomes comparable to their potential energy at the average interparticle separation, $kBTDIH∼e2/a$⁠.²⁴ Figure 1 shows the kinetic temperature evolution of electrons (black line) and ions (blue dashed line). A rapid increase in the electron temperature due to DIH occurs immediately, and is saturated within 1–2 $ωpe−1$⁠. A monotonic but slower increase in the electron temperature at later times is attributed to TBR heating. The rate of each processes can be determined from the slope of a linear fit in each time interval (dashed line for DIH and dashed-dotted line for TBR). As expected, ion heating is negligible on this short timescale.

During DIH, electrons gain kinetic energy largely via the ballistic motion associated with their electrostatic attraction to the nearest ion. Since the timescale associated with DIH is so short, electrons are not thermalized during this period, and the electron velocity distribution is not expected to be Maxwellian. Figure 2 shows the electron velocity distribution in the x-direction. Black and red solid lines show the distribution of electron velocities at 0.2 $ωpe−1$ and 1 $ωpe−1$⁠, respectively. Fits to a Maxwellian, obtained by matching the height of the velocity distribution, are shown as dotted lines. It is clear from the plot that the actual velocity distribution is not Maxwellian. Temperature is not thermodynamic, but rather interpreted in terms of the kinetic temperature of Eqs. (2a) and (2b).

To quantify the electron kinetic energy gained by accelerating toward the nearest ion, we provide a simple model consisting of an immobile ion and an electron with zero initial velocity. The initial distance between the electron and ion was chosen from the initial nearest neighbor electron distribution around ions from a distribution where both charged species were generated in a random homogeneous manner, as in the MD simulations. The equation of motion $mer̈=−∇vei$ was solved for each individual nearest neighbor pair (picked one pair at a time). The kinetic energy was calculated as the electron moved directly toward the ion over a time interval $Δt=ωpe−1$⁠. Figure 3(a) shows an example electron trajectory, and Fig. 3(b) the corresponding electron kinetic energy. In this example, the initial separation was $r=0.7a$⁠. After a representative sample of initial nearest neighbor distance configurations was explored, the distribution of electron velocities was obtained.

Figures 4(a) and 4(c) show the averaged electron velocity distribution dN/N (probability of finding electrons in velocity range of v_x to $vx+50$ m/s) obtained from the model and from MD simulations, respectively. The velocity distributions of the model and simulations qualitatively agree. The average kinetic temperature over three separate time intervals (⁠$0.1–0.2,0.5–0.6$⁠, and $0.9–1.0ωpe−1$⁠), calculated using Eq. (2a), is also indicated in the figure. While temperatures obtained from the model were 0.0784 K, 0.315 K, and 0.565 K (i.e., $1/Γe=0.0628,0.2522$⁠, and 0.4523, respectively), the temperatures obtained from the MD simulations were approximately twice as high with values 0.183 K, 0.802 K, and 1.225 K (i.e., $1/Γe=0.1465,0.6421$⁠, and 0.9807, respectively). The additional heating observed in the MD simulations, compared with the binary interaction model, must be associated with many-body interactions. To test this, we carried out two particle MD simulations for several pairs of nearest neighbors and compared the kinetic energy gained by the electron with its value from the full MD simulation in the presence of all other charges. In each case, we found that the kinetic energy obtained by electron in full MD is greater than what it attains in two particle simulation.

Figure 5 presents the dependence of the electron temperature evolution profile on (a) the repulsive core parameter α, (b) the ion to electron mass ratio $mi/me$⁠, (c) the choice of simulation models, and (d) the choice of pseudo-potential for electron–ion interaction. Figure 5(a) suggests that the T_e evolution profile does not appreciably change for values of the repulsive core parameter $α<0.01$ over the time interval of interest. Here, close interactions $r<αa$ are sufficiently rare that further decreasing the value of α causes negligible changes to the heating. For the remainder of this work, we choose the value $α=0.01$⁠. Figure 5(b) shows that the T_e profiles are indistinguishable for the mass ratio of 1000 and $1.62×105$⁠. This suggests that for any realistic mass ratio $mi/me≥1000$⁠, electron and ion dynamical timescales are sufficiently well separated that ion dynamics does not affect the electron heating process. For the remainder of this work, we choose $mi/me=1.62×105$⁠, which corresponds to strontium plasma.

Figure 5(c) provides a comparison of T_e evolution profiles obtained from the one-component plasma (OCP) model (blue dashed line) with those obtained from two-component plasma (TCP) models. The black solid line represents the T_e profile when both (electron and ion) species were dynamic during the simulation. The red dashed line represents T_e profile when only electrons were dynamic and ions were held stationary. Both TCP models provide the same T_e evolution profile, suggesting that ion motion does not influence electron heating over this time scale. This agrees with earlier simulations made by Kuzmin and O'Neil.³⁷ An OCP model predicts a much lower $TDIH$ and suggests that the electrons could remain strongly coupled (⁠$Γe>1$⁠) after the saturation of DIH. The OCP model does not agree with the TCP models or with experimental observations where electrons are observed to be weakly coupled.⁴⁰ An expected reason is that the early stage (⁠$Δt<ωpe−1$⁠) heating of electrons is dominated by kinetic energy gain by attraction toward nearest neighbor ions, which is absent in the OCP DIH models.²⁴ The Yukawa-OCP model has been found to accurately explain ion temperature evolution, in agreement with experiments.²¹ This is likely because at late times (⁠$Δt∼ωpe−1$⁠), electrons have already moved toward ions, so that ions move in the presence of a polarized (screened) charge cloud created by electrons (i.e., a Yukawa one-component plasma, YOCP).

To investigate the role of pseudo-potential choice on the electron heating mechanism, we modeled the electron–ion interaction using three different forms: $vei=e2/r2+(αa)2$ (O' Neil's form³⁷) $vei=(e2/r)(1−exp (−r/αa))$ (Hansen' form⁴¹) and the Kelbg form [Eq. (3b)].⁴² Figure 5(d) shows that the T_e profiles obtained from the simulations using all three different forms of pseudo-potential are in good agreement for a small value of the repulsive core parameter $α=0.01$⁠.

When an external magnetic field is applied, electrons gyrate around the field lines. At high field strength, the gyroradius is so small that the electrons move primarily in one-dimension, parallel to the magnetic field, and consequently, gain kinetic energy in only this direction. Because electron kinetic energy gain is exceptionally slow in the cross field direction, this leads to an effective drop in the total kinetic T_e in comparison to an unmagnetized plasma.

Figure 6 shows the effect of increasing magnetic field strength in the MD simulations. At a sufficiently strong magnetic field, the total kinetic temperature due to the DIH saturates (in $∼ωpe−1$⁠) at approximately one-third of the value in the case of no magnetic field. We further observe a drop in TBR heating rates at later times with the increase in magnetic field strength. Figure 7 shows that the TBR heating rate saturates at one-third of its value in the unmagnetized case. The TBR heating rate is calculated as the slope of the T_e profile in the time interval $1.5ωpe−1<t<10ωpe−1$ (see Fig. 1).

Due to the strong magnetic field, a strong temperature anisotropy develops in the plasma.¹² Figure 8 shows the evolution of the kinetic electron temperature in the presence of a strong magnetic field (⁠$βe=200$⁠). It is evident that the parallel kinetic temperature (blue dashed line) rises as it does in an unmagnetized plasma, and saturates once the kinetic energy of the electrons in this direction becomes comparable to the potential energy of particles at an average interparticle separation, i.e., $kBT||∼e2/a$⁠.

However, because of the strong magnetic field, the heating time in the perpendicular direction to the field becomes so long that there is essentially no increase in the kinetic temperature in this direction. The magenta dashed-dotted line in the Fig. 8 shows the evolution of the perpendicular kinetic temperature. The black line in Fig. 8 shows the total kinetic temperature from Eq. (2c). Out of three directions, the temperature remains negligible in two of them due to the strong magnetic field. This leads to a total electron temperature which is close to one-third of electron temperature in the unmagnetized case. Typically, the ratio $T∥/T⊥≈30$ at 10 $ωpe−1$ with the magnetic field strength of $βe$ = 200.

The temperature anisotropy can also be observed in the velocity distribution function; see Fig. 9. Subplots (a) to (d) show the averaged electron velocity distribution over four time intervals in the range from 0 to 1 $ωpe−1$⁠. A spherically symmetric (circularly symmetric in the 2D plot) electron velocity distribution is observed when there is no magnetic field. In contrast, the electron velocity distribution is highly asymmetric in the presence of a strong applied magnetic field [subplots (e)–(h)].

The spontaneously generated temperature anisotropy will relax at a rate that depends on the coupling strength Γ and the magnetic field strength β.^10–12 Recent results³¹ (for the one component plasma) suggest that for a coupling strength of $Γ=0.1$ and magnetic field strength of β  =  100, the temperature isotropization time is $∼106ωpe−1$⁠. This would be a sufficiently long delay that the electron temperature anisotropy would last for the complete ultracold plasma lifetime, which can last up to 250 μs (or few hundreds of $ωpi−1$⁠).

Similar to the unmagnetized case, the early time (⁠$Δt<ωpe−1$⁠) kinetic energy gained by electrons is primarily due to their ballistic motion in the presence of a nearest-neighbor ion (with approximately half of the contribution being due to the multi-body interactions). At high field strength, the ballistic motion is not directly toward the ion, but rather is restricted to the direction parallel to the magnetic field. This effectively reduces the average electrostatic potential energy accessible to the electrons in comparison with an unmagnetized plasma.

To demonstrate this, we again introduce a simple model. The equation of motion $medve/dt=−∇vei−e(ve×B)$ was solved for each individual nearest neighbor pair (picked one at a time) from the nearest neighbor electron distribution (as described in the unmagnetized case) and the electron kinetic energy was calculated. An example trajectory is shown in Fig. 3(c), and the corresponding kinetic energy in Fig. 3(d). Similar to the unmagnetized case, we use this model for all possible nearest neighbor distances between electrons and ions obtained from a random initial distribution. Figure 4(b) shows that the electron velocity distribution in the parallel direction obtained from this model qualitatively matches the one obtained from the MD simulations with $βe=200$ [Fig. 4(d)]. The time averaged (over $0.1ωpe−1$⁠) velocity distributions obtained from both ways do not follow the Maxwellian. The average parallel kinetic temperatures in the time intervals $0.1–0.2, 0.5–0.6$⁠, and $0.9–1.0ωpe−1$ from the model are 0.0784 K, 0.315 K, and 0.565 K (i.e., $1/Γe=0.0628,0.2522$⁠, and 0.4523), respectively (the same as the unmagnetized case because there is no force due to the magnetic field in the parallel direction). The average parallel temperatures obtained using MD simulations during these same intervals are 0.145 K, 0.672 K, and 1.095 K (i.e., $1/Γe=0.1161,0.5380$⁠, and 0.8766).

Finally, as a measure of increased coupling strength, we calculate the radial distribution function (RDF) for the electrons. The RDF, or pair correlation function, is a structural statistical property of the medium which provides the probability of finding a particle at some distance r from a reference particle normalized to the background. In a homogeneous and isotropic plasma, this depends only on the radial distance from the reference particle. In molecular dynamics, the RDF is calculated by binning distances between all the particles to determine how many particles lie within distance r to r + dr away from the reference particle. Further, these values are normalized by the particle histogram for an ideal gas which has no correlations. In Fig. 10, we plot RDFs of electrons for an unmagnetized plasma (dashed black line) and for a plasma with magnetic field of strength $βe=200$ (blue solid line). Each RDF is averaged over duration $3−5ωpe−1$⁠, a timescale at which electron DIH is saturated. The fall of the RDF toward $r→0$ is steeper in the case with a magnetic than in the case without. This steepness is a signature of a higher coupling strength. This fact can be verified by RDFs obtained from the OCP for Γ = 1 and Γ = 3; see the inset of Fig. 10.

Experiments focused on electron dynamics make measurements on 1–200 ns timescales.^43,44 Such experiments are designed in such a way that the electrons possess a very small initial kinetic energy due to the photoionization of ultracold atoms.⁴³ Here, we suggest that an external magnetic field of approximately one-tenth of a Tesla or higher could reduce the subsequent DIH and TBR heating to one-third of the value in the absence of a magnetic field. This will increase the effective electron coupling strength by a factor of three. Recent experiments (at electron time scales (ns)) have reported an electron coupling strength $Γe∼0.35$ (⁠$Te∼1.6$ K),⁴⁴ so the presence of an external magnetic field creating $βe∼50$ (or higher) will be sufficient to observe the $Γe∼1$ on this timescale.

Experiments focused on ion dynamics make measurements at μs timescales. Over this period (⁠$>103ωpe−1$⁠), electrons continuously gain kinetic energy due to the TBR heating process. A magnetic field of $βe∼50$ (0.16 T for $ne=1×1014$ m⁻³) is sufficient to reduce the DIH and TBR heating by a factor of three over a timescale of several $ωpe−1$⁠, but to extend any relaxation of the temperature anisotropy to μs timescales experiments will need a strong magnetic field satisfying $βe>100$⁠. The sustained temperature anisotropy may extend the increased electron coupling strength (compared with unmagnetized case) to the μs timescale.

Though UCPs develop a strong temperature anisotropy due to the applied magnetic field, such plasmas are not expected to be susceptible to temperature or pressure anisotropy driven instabilities (such as fire-hose^45,46 or Weibel⁴⁷) These instabilities are known to exist at high values of plasma beta (⁠$βpl=nkBT/(B2/2μ0)≳1$⁠). The β_pl is related to the magnetic field strength as $βpl=(ṽ/c)2/β2$⁠. For UCPs, temperature is a very small quantity, $ṽ/c≪1$ (⁠$ṽ$ is the mean thermal velocity). Also, we are studying UCPs under strong magnetic field strengths (⁠$β>100$⁠). Under these conditions, β_pl becomes a very small value.

In summary, we addressed the issue of the electron heating and its reduction through the application of an external magnetic field in an electron-ion ultracold plasma. Using classical MD simulations, we showed that in the presence of a strong magnetic field (⁠$βe≥50$⁠), electron DIH is reduced by a factor of three compared to its value in an unmagnetized plasma. We also showed that the electron TBR heating rate was reduced by the factor of three at a similarly strong magnetic field strength. The reduction in heating occurs due to electron motion being constrained to one-dimension (parallel to magnetic field) resulting in little heating in the perpendicular direction. The ultracold plasma parameters suggested that a magnetic field of approximately a Tesla would be sufficient to see this effect in experiments at timescales of the lifetime of UCP. This suggests the possibility to observe an increased effective coupling strength of electrons along with the moderately coupled ions in UCP experiments. We also see that a strong external magnetic field causes a strong temperature anisotropy $T∥/T⊥∼30$⁠. Using a simple model consisting of individual electron-ion pairs, it was shown that approximately one-half of the electron disorder induced heating occurs due to electron motion directly toward its nearest neighbor ion at the early times (⁠$0<t<1ωpe−1$⁠).

We thank N. Shaffer for useful discussions. This material is based upon work supported by the Air Force Office of Scientific Research under Award No. FA9550-16-1-0221. It used the Extreme Science and Engineering Discovery Environment (XSEDE),⁴⁸ which is supported by NSF Grant No. ACI-1053575, under Project Award No. PHY-150018. This research was supported in part through computational resources provided by The University of Iowa.