Effects of Disorder on Electron Heating in Ultracold Plasmas

Starting from the beginning of their research in the early 2000's, the ultracold plasmas were considered as a promising tool to achieve considerable values of the Coulomb coupling parameter for electrons. Unfortunately, this was found to be precluded by a sharp spontaneous increase of temperature, which was often attributed to the so-called disorder-induced heating (DIH). It is the aim of the present paper to quantify the effect of spontaneous heating as function of the initial ionic disorder and, thereby, to estimate the efficiency of its mitigation, e.g., by the Rydberg blockade. As a result of the performed simulations, we found that the dynamics of electrons exhibited a well-expressed transition from the case of the quasi-regular arrangement of ions to the disordered one; the magnitude of the effect being about 30%. Thereby, we can conclude that the two-step formation of ultracold plasmas - involving the intermediate stage of the blockaded Rydberg gas - can really serve as a tool to increase the degree of Coulomb coupling, but the efficiency of this method is moderate.


I. INTRODUCTION
The so-called ultracold plasmas are neutral systems of charged particles with typical electronic temperatures from a few to several tens of Kelvin, which are obtained by a photoionization of gases cooled in the magneto-optical traps; e.g., reviews [1][2][3][4] .The experimental realization of such plasmas became feasible in the very late 1990's and early 2000's, and they opened a new area of research in the non-ideal plasma physics 5,6 .
It was initially expected that the extremely high values of the Coulomb coupling parameter (where K and U are the kinetic and potential energies of an electron) can be achieved in this way.Really, if energy of the ionizing laser irradiation was chosen to be slightly above the ionization threshold of the cold neutral atoms, the initial kinetic energy of the released electrons would be very low.Therefore, it was expected in the first experiments 7 that very large values of the coupling parameter (1) could be obtained, e.g., tens or hundreds.Unfortunately, it was quickly recognized that the situation is not so simple: In fact, the cold photoelectrons are quickly accelerated by the electric fields of nearby ions, and their temperature spontaneously increase by many times.
A simple pictorial explanation of this effect is the socalled disorder-induced heating (DIH) 8 : The charged particles (firstly, electrons and, at the longer time scale, also the a) Corresponding author's electronic mail: dumin@pks.mpg.de,dumin@yahoo.comb) Electronic mail: lukashenko@dec1.sinp.msu.ru,a_lu@mail.ruions) tend to move to the positions with minimal potential energy.As a result-since the total energy of the system should be conserved-kinetic energy of the particles will increase.Therefore, DIH looks like an unavoidable effect, limiting the temperature from below.
However, it was suggested by the same authors 8 that one can get around the DIH effect by preparation of the system of charged particles in the "correlated" state with a reduced Coulomb energy.An attempt of realization of this approach was undertaken a decade later in the experiment 9 : Namely, cold neutral atoms were initially transferred to the state of Rydberg blockade, where the already excited atoms-due to their strong electric fields-shift energy levels of the nearby atoms, thereby prohibiting their excitation to the same Rydberg state [10][11][12] .As a result, a quasi-regular arrangement of the excited atoms is formed, where their close location to each other is excluded (e.g., Fig. 3 in paper 13 ).Next, after photoionization of Rydberg atoms, the resulting ions will also be well separated from each other.Therefore, one can expect that DIH will no longer take place, because the system of ions is already in the "quasi-crystalline" state with minimal potential energy.
The experiment 9 was really able to trace a plasma formation from the blockaded Rydberg gas, but it remained unclear if the resulting electron temperature was appreciably reduced (and, respectively, the Coulomb coupling parameter increased).This was caused, firstly, by the fact that the ionization proceeded mostly by the spontaneous avalanche process (and, therefore, did not strictly preserve the initial quasicrystalline arrangement of the Rydberg atoms) and, secondly, by the insufficient diagnostic capabilities to measure the electron temperature (S.Whitlock, private communication).Besides, there was no a clear theoretical prediction how much could be the magnitude of the expected reduction in temperature.Surprisingly, while the DIH mechanism was discussed FIG. 1. Examples of the ionic arrangement in the xy-plane (bottom row) and the corresponding histograms of the interparticle separation (top row) for the quasi-regular distributions with different degrees of disorder σ reg .For convenience, a grid of dotted lines illustrates a characteristic space per one particle in the perfect cubic lattice.
for the first time about 20 years ago 8 , it was subsequently studied mostly for ions [14][15][16][17][18] , and there are no reliable calculations of its influence on the electron temperature till now.
So, it is the aim of the present paper to quantify the corresponding effect in electrons as function of the initial disorder.In fact, this problem was partially touched in our previous work about the clusterized plasmas 13 .As a particular case of the nontrivial arrangement of the background ions, we considered there also their quasi-regular distribution and found that the resulting electron temperature was somewhat reduced.However, the degree of such a reduction was comparable to the uncertainty (r.m.s.variation) of that simulations.Besides, it remained unclear how the corresponding effect depends on the degree of disorder and, particularly, how it tends to the purely random case when the disorder becomes sufficiently large.In the present work-due to much better accuracy of the simulations-all these issues will be resolved.

II. NUMERICAL MODEL
A quite general model of ionic background with different degrees of disorder can be formulated by the following way: Let us consider initially a perfectly cubic lattice of ions.The size of its cell l will be used from here on as the unit of length, and all the distances will be normalized accordingly.Next, let each ion be shifted from its original position by a distance given by the normal (Gaussian) distribution with r.m.s.deviation σ reg , as illustrated in the bottom row of Fig. 1.Then, at σ reg ≪ 1 (in dimensionless units) the ionic distribution is quasi-regular; but it becomes more and more disordered when σ reg increases; and finally, at σ reg ∼ 1, we evidently get a com- pletely random distribution.
These properties are well expressed in the histograms of interparticle separation, shown in the top row of Fig. 1.For example, at σ reg ≪ 1 one can see a series of sharp peaks, which represent a set of the preferable interparticle distances in the quasi-regular lattices.Next, when the disorder increases, these peaks are gradually smoothed out; and finally, at σ reg ∼ 1, the histograms take the Gaussian shape, which is typical for a random distribution.
To study dynamics of electrons against the abovementioned kinds of ionic background, we shall numerically integrate their equations of motion: where r i and R i are the electronic and ionic coordinates, respectively; and m e and e are the electron mass and charge.The ions are assumed to be immobile, because we are interested only in the sufficiently short time intervals.The initial electron coordinates r i (0) are given by the uniform statistical distribution (i.e., the electrons are randomly distributed with a uniform average density).The initial electron velocities v i (0) are specified by the normal (Maxwellian) law with the r.m.s.variation σ v : where α is x, y, or z.In the ideal plasmas, the abovementioned parameter is evidently proportional to the square root of electron temperature.However, the situation can be more tricky in the non-ideal case.This is the reason why we prefer to speak here just about the r.m.s.variation.Strictly speaking, the initial values of the electron coordinates and velocities strongly depend on the details of the ionization process.For example, in the case of instantaneous photoionization, the initial electron positions will be strongly correlated with the positions of ions 19 , while distribution of the electron velocities will look like the delta-function.However, if the photoionization takes some time, the released electrons should be somewhat mixed between the ions and the distribution of their velocities smoothed out, resulting in the Maxwellian form.
Next, we shall use the perfectly reflective boundary conditions, and the Coulomb interactions will be calculated within the fixed simulation box.An alternative choice-used in the most of the previous works on ultracold plasmas-is to employ the periodic boundary conditions, when a particle leaving the simulation box through one of its boundaries simultaneously enters it through the opposite boundary.The Coulomb interactions are usually calculated in such a case by the "wrapping" algorithm 19 , when each particle interacts with other particles within a moving cube of size ±L/2, centered at that particle (where L is the total size of the simulation box).
In principle, none of these options is perfect: The reflective boundaries evidently affect bulk properties of the plasmas.However, as was shown in our previous work 13 , averaging over a sufficient number of initial conditions substantially mitigates this problem.On the other hand, the periodic boundary conditions with the above-mentioned wrapping procedure should not distort the bulk properties.Unfortunately, a closer inspection shows that swapping of a charged particle between the boundaries of a moving box results in the abrupt unphysical change in the direction of Coulomb force between the respective pair of particles.From this point of view, the reflective boundaries look better because their effect on the bulk properties has a clear physical meaning 20 , and it will evidently disappear with increasing the number of particles.
Besides, in the context of our simulations, a decisive advantage of the reflective boundary conditions is that they are well consistent with the algorithms of numerical integration with the adaptive stepsize control (ASSC), such as the subroutines odeint, rkck, and rkqs from the Numerical Recipes 21 .These subroutines were already used in our previous work 13 and demonstrated a perfect performance.Particularly, they are able to work with the "true" (singular) Coulomb potentials, without a need for their cut-off or "softening" at small distances.This excludes any artifacts caused by the modified potentials.(Let us mention that dealing with singular potentials without the ASSC algorithms requires the huge computational resources 22 .)Unfortunately, ASSC becomes unreliable in the case of periodic boundary conditions because of the above-mentioned abrupt jumps of the Coulomb forces.
All the results will be presented below in dimensionless units, introduced by the following way: a unit of length is a mean distance between the ions l (the corresponding unitary cells are marked by dotted lines in Fig. 1); a unit of time, up to numerical factor on the order of unity, is the inverse plasma frequency, and a unit of energy is the characteristic Coulomb energy at the interparticle distance, U = e 2 /l.The dimensionless quantities expressed in these units will be denoted by hats.
Let us emphasize that the unit of time (4) has a deep physical meaning: If the initial plasma state is substantially overcooled, i.e., the electron velocities are sufficiently small, then each electron will fall onto the nearest ion.The corresponding Keplerian trajectory will be strongly elliptical, and the ion will be located at the focus of this ellipse opposite to the initial position of the electron.Next, it should be taken into account that a period of revolution along the Keplerian trajectory depends only on its major axis a and equals 23 : e a 3/2 /e.A typical distance from the electron to the nearest ion should be about l/2.Then, taking this quantity as the major axis, we get: T Kep = (π/4)τ, where the numerical factor is very close to unity.Therefore, a characteristic fall time of any electron should be about one half of the unit of time given by formula (4).In this sense, the accelerated electrons behave "synchronously" and-up to numerical factor about unity-are characterized by the inverse plasma (Langmuir) frequency.(This argumentation evidently refers only to the electron dynamics.Treatment of ions, especially when they experience the Debye screening, requires a more complex analysis, which was undertaken, for example, in paper 24 .)Besides, the above-mentioned "synchronous motion" does not imply the applicability of the mean-field approximation, which can be used only for the sufficiently ideal and equilibrium plasmas.
Finally, the electron temperature is assumed to be related to the kinetic energy per electron by the same formula as for an ideal gas, T e = (2/3) K/k B .Although we deal here with a substantially non-ideal case, it was found in paper 19 that such definition reasonably agrees with a more elaborated derivation of T e , based on the approximation of the simulated velocity distributions by the Maxwellian ones.A theoretical explanation of this fact can be found in paper 25 .

III. RESULTS OF THE SIMULATIONS
Our simulations were performed for a system of 1000 particles of each kind (electrons and ions), i.e., the simulated volume was composed of 10 × 10 × 10 unitary cells, as depicted in Fig. 1.This is 8 times greater than in the previous simulations of the clusterized plasmas 13 : since influence of the quasi-regular ionic arrangement is a finer effect, we had to increase the number of simulated particles.
In the particular case presented below, we used σv = 0.3, which implies that the initial kinetic energy of electrons was about an order of magnitude less than their potential (Coulomb) energy, i.e., the plasma was substantially overcooled.It should be expected that final results will be insensitive to the particular degree of overcooling, because the major FIG. 2. Individual profiles of the electron kinetic energy as function of time (thin curves) and their average behavior (thick curves) at σ reg = 0.01 and 1 (top and bottom panels, respectively).effect comes from the subsequent spontaneous heating.(Dependence of the simulations on σv is discussed in more detail in Appendix A.) In other words, the plasma was already in the non-ideal (strongly-coupled) state, and our aim was to check how this state will survive against the subsequent heating.
The simulations were performed for the following degrees of disorder: σ reg = 0.01, 0.03, 0.1, 0, 3, 1 and 3, as well as for the purely random ionic distribution.To get the statistically significant results, five versions of initial conditionsrandomly generated both for the ions and electrons-were used for each of these values.Despite dealing with singular interparticle potentials, the algorithm of the adaptive stepsize control enabled us to get a sufficiently high accuracy of integration.It was estimated, as usual, by a conservation of the total energy of the system.In the worst case, such error was 1.4%, but usually by one or two orders of magnitude better.
Examples of the simulated temporal behavior of the dimensionless kinetic energy of electrons are presented in Fig. 2. To avoid obscuring the plots with a lot of sharp peaks, caused by the close interparticle collisions, they were smoothed out over a running window of width ∆t = 0.1.In principle, a general shape of these curves-a quick initial jump at the timescale t ≈ 0.5 followed by a much longer gradual increase-was well known already from the pioneering works by Kuzmin and O'Neil 26,27 .It is the aim of our study to reveal how they depend on the arrangement of the ionic background.Figure 3 shows the average temporal profiles of the electron kinetic energy for the entire variety of the disorder parameters σ reg , ranging from an almost perfect ionic lattice to the completely random distribution.One can see here three types of the curves: Type I corresponds to the small values of σ reg (0.01, 0.03, and 0.1), i.e., the weakly distorted lattices.Type II with σ reg = 0.3 is the intermediate case, correspond- ing to the moderately distorted lattice.At last, the curves of Type III belong either to the cases of strongly distorted lattices, σ reg = 1 and 3, or to a completely random distribution.As could be naturally expected, Type II indicates a noticeable change in the temporal behavior of the electron kinetic energy K .Namely, when σ reg increases, the initial jump (at the timescale 0 ≤ t 0.5) becomes more pronounced.The subsequent gradual increase in K at t 0.5, in principle, also changes but insignificantly.
Referring to the histograms of interparticle separation, presented in the top row of Fig. 1, one can see that at small values of σ reg (e.g., 0.03 and 0.1) they exhibit a series of sharp peaks, which are typical for the crystalline-like structures.On the other hand, at the large values of σ reg (e.g., 1) their shape closely resembles the purely-random Gaussian distribution.It is interesting that in the intermediate case σ reg = 0.3 the his- togram is almost Gaussian, but the curve K (t) in Fig. 3 exhibits a clearly distinct behavior.In fact, a closer inspection of the patterns of ionic arrangement in the bottom row of Fig. 1 shows that the case σ reg = 0.3 preserves some features of the quasi-regularity: namely, there are no occasional "voids", typical for the purely random distributions.
To quantify the three above-mentioned types of the temporal behavior K (t), we calculated average values of the electron kinetic energy at two time intervals: t ∈ [0.5, 1] and t ∈ [9, 10].The first of these quantities characterizes a magnitude of the initial jump of K, which was often attributed just to the disorder-induced heating; while the second quantity represents the total increase in K after a sufficiently long period of evolution, which includes also a subsequent heating of plasma FIG. 4. Average values of the electron kinetic energy established at the time intervals [0.5, 1] (top panel) and [9, 10] (bottom panel), as well as their r.m.s.deviations, as functions of the disorder parameter σ reg .due to recombination.First of all, this is the three-body recombination, which was discussed by various authors starting from the first experiments with ultracold plasmas.However, as follows from our recent simulations of the clusterized plasmas 13 , the multi-body processes (involving two electrons and more than one ion) should be also important in the non-ideal case: Really, since the rate of recombination substantially depended on the degree of clusterization, a simultaneous scattering of two electrons by two or more ions came into play there.The above-mentioned quantities are presented in Fig. 4 as functions of σ reg .
estimate a statistical significance of the simulations, we plotted in this figure also the r.m.s.variations of the average energy with respect to various versions of the initial conditions σver and with respect to time σt over the corresponding intervals (for mathematical details, see Appendix B).These quantities are shown by the dashed and dotted lines, respectively.As is seen, σver is appreciably smaller than σt .In other words, five versions of the initial conditions are quite sufficient for the reliable averaging.On the other hand, the temporal r.m.s.variation σt is evidently unavoidable, and its effect can be reduced only by increasing the number of particles in the simulated system.

IV. DISCUSSION AND CONCLUSIONS
1.The main result of our study is the identification of the clear transition of the electron dynamics from the case of quasi-regular ionic background (e.g., caused by the Rydberg blockade, as discussed in the Introduction) to the random one.Let us mention that some reduction of the electron temperature in the regularized ionic arrangement was found already in our previous work 13 , but this effect in K(t) was comparable to the uncertainties σver and σt ; see Table 1 and Fig. 8 in that paper.Besides, it was rather surprising why there was no a well-expressed transition of K(t) to the purely random case when σ reg tended to unity.In the present study-due to the enhanced accuracy of simulations-this puzzle was resolved, and the clear transition was identified.
Unfortunately, the effect of reduction of the electron temperature is not so large: as is seen in Fig. 4, it is approximately 30% both immediately after the sharp jump, occurring at the timescale of one half of the dimensionless unit, and at the longer time interval, when a heat release due to recombination came into play.
It is interesting to mention that yet another method to reduce the electron temperature (and, thereby, to increase the Coulomb coupling parameter) was suggested in papers 28,29 .This is adding the Rydberg atoms with binding energies |E b | (2 − 3) k B T e into the ultracold plasmas.As a result, their inelastic collisions with free electrons will lead to further excitation of the atoms and cooling of the electrons.It was found in the above-cited papers that the overall efficiency of such a process should be about 20-30%, i.e., actually the same as in the method based on the Rydberg blockade 9 .
At last, reduction of the electron heating rate by a factor of 3 was obtained by simulating the strongly-magnetized ultracold plasmas 30 , when the electron motion was effectively constrained to a single dimension.However, since such plasmas involve a strong and long-lasting temperature anisotropy, they represent substantially different physical systems as compared to the ones considered in the present paper.
2. Yet another unexpected finding in our simulations is that there was a rather strong spontaneous heating of the ultracold plasma even in the case of almost regular lattice, i.e., when there was essentially no disorder.So, we believe that the concept of disorder-induced heating (DIH) may have a limited scope of applicability.An alternative explanation can be based, for example, on the concept of "virialization" 25 , which was suggested even before DIH but did not attract attention till now.In that case, a sharp increase in temperature is attributed to the establishment of virial distribution between the kinetic and potential energies and, therefore, the rate of initial heating should not depend appreciably on the degree of disorder.Of course, a more detailed analysis of the behavior of different kinds of energy should be performed to discriminate more definitively between the concepts of DIH and virialization; this is planned to be done in a separate paper.
3. A rather restrictive assumption in our simulations was imposition of the reflective boundary conditions, which was necessary to ensure a smooth operation of the numerical integration algorithms with the adaptive stepsize control (ASSC).In principle, the reflective walls might affect the bulk properties of the simulated plasmas (especially, when the number of particles in the model is not so large, and a considerable fraction of them is located near the walls).However, as was shown in our previous paper 13 , this problem can be substantially mitigated by averaging over a sufficiently large set of initial conditions.Then, the results become almost independent on the number of particles in the simulation cell.
Of course, it would be desirable to use in future simulations the periodic boundary conditions, which are more relevant from the physical point of view.A promising approach to reconcile the periodic boundary conditions with ASSC might be the so-called Ewald summation of the Coulomb interactions 31,32 .Then, there should be no sharp jumps of forces when a particle is transferred from one boundary of the cell to another, and ASSC should work much better.But the implementation of such approach evidently requires a lot of additional work.red curves, respectively.
As could be expected, in the case of strongly overcooled initial state ( σv ≪ 1), a subsequent evolution of the electron kinetic energy is almost independent of σv : Really, as is seen in both panels, the blue ( σv = 0.1) and cyan ( σv = 0.3) curves are almost indistinguishable from each other.Moreover, the cyan curve is often invisible at all, since it is obscured by the blue one.
However, when the initial plasma state was taken to be moderately coupled (non-ideal), e.g.σv = 1.0, there was a noticeable change in the temperature evolution: it is seen that a green curve lies above the blue and cyan ones.This is not surprising because the electrons possessed an additional energy already in the initial state.Besides, despite a sharper initial jump, a subsequent increase of the temperature proceeds a bit slower.This is also rather natural, because it is well-known that the rate of three-body recombination considerably drops with increase in the electron temperature 33 (and the same should be expected for the multi-body recombination).
At last, when the plasma was initially almost ideal (for example, σv = 3.0, i.e., its kinetic energy exceeded the potential one by an order of magnitude), the corresponding red curve lies well above other curves.Moreover, it is almost horizontal, i.e., a heating due to recombination is negligible.This case is evidently irrelevant to the situation discussed in the main body of the paper.
In principle, the fact that dynamics of the overcooled plasmas is almost independent of the particular values of initial temperature was mentioned already in paper 30 .However, that simulations were performed with the "soft-core" potentials (i.e., Coulomb interactions between the electrons and ions were artificially cut off at small distances).Since the dynamics of overcooled electrons involves a lot of very close encounters with ions, it is unclear in advance if the abovementioned independence will survive in the case of more realistic Coulomb interactions.So, as follows from the present simulations, this really takes place.

FIG. 3 .
FIG.3.Average (over 5 realizations) temporal profiles of the electron kinetic energy K for the entire set of the disorder parameters σ reg .