Hot electrons between cold walls

A cryogenic pellet en route through a hot plasma ablates under the heat flux from the ambient electrons. The incoming hot electrons bring a negative charge to the pellet, and the resulting electrostatic potential can limit the heat flux unless the pellet releases cold electrons in return. The hot electrons in fusion-grade plasmas are in a multiple kiloelectron volt energy range. Consequently, they tend to build up a multikilovolt potential difference between the pellet and the ambient plasma. These characteristic values of the hot electron energies and the resulting potential are typically much greater than the ionization potential and the work function of the pellet material, which indicates that the pellet emits the cold electrons freely to balance the incoming current of the hot electrons. A very low electric field at the pellet surface is thus sufficient to drive the needed return current. Suppression of the surface electric field due to free emission is similar to the one in the Child–Langmuir problem,¹ and it is very different from the low-temperature regime of plasma–surface interaction, in which there are essential constraints on emission from the surface.^2,3 The electron heat flux to the pellet can still be large in the absence of net current (of the order of $n0TT/m$ where $n0$ is the plasma density, $T$ is the temperature, and $m$ is the electron mass). This flux is much larger than the ion thermal flux, which suggests that the plasma electrons can cool down quickly in a collisionless regime without any significant movement of the ions while the pellet travels across the flux tube. In this paper, we use a simple model to study such rapid cooling and to assess the role of electrostatic shielding.

The characteristic time scale of flux surface crossing by the pellet is $τp=2rp/vp$⁠, where $rp$ is the pellet radius and $vp$ is the pellet velocity. It is instructive to compare this time scale with the electron and ion transit times along the helical flux tube of width $2rp$⁠. The characteristic transit times are $τe=2πN(a)R/vT$ and $τi=2πN(a)R/cS$⁠, where $R$ is the tokamak major radius, $a$ is a radius of a given magnetic flux surface averaged over the poloidal angle, $vT=T/m$ is the electron thermal velocity, $cS=γT/M$ is the ion sound velocity, $M$ is the ion mass, and the number of turns $N(a)$ depends on whether or not the safety factor $q$ is rational in the vicinity of a pellet (see also Ref. 4). Let us estimate the characteristic timescales of the problem for the pellets of interest for disruption mitigation in ITER.^5,6 For $Te=104 eV$⁠, $vp=104 cm/s$⁠, and $rp=1 cm$ we have: $τp=2×10−4 s$⁠, $τe=0.9×10−6N(a) s$⁠, and $τi=0.32×10−4N(a) s$⁠. We note that $N(a)=4a/rp∼102$ on an irrational flux surface, while $N(a)≈ma∼10$ on a rational one with $q=ma/na$⁠. Comparing the timescales, we see that $τi/τp≫1$ everywhere except for rational surfaces with low $ma$⁠. In contrast, we have $τe/τp≪1$ in the plasma core (⁠$a≤50 cm$⁠) and on each of the rational flux surfaces. Consequently, all plasma electrons in the core and near rational flux surfaces have enough time to interact with the pellet, whereas the ions are too heavy to do that. This reiterates the need to analyze electron cooling in the absence of ion movement.

The magnetic field lines of the helical flux tube in a tokamak connect two opposite sides of the pellet. The strong magnetic field suggests a one-dimensional approximation. We, therefore, view two sides of the pellet as two freely emitting “walls” with a plasma slab in between. An incoming hot electron loses its energy to the pellet, whereas an emitted cold electron gains energy in the electrostatic sheath near the wall. As a result, the wall absorbs the perpendicular energy of the incoming electrons and a part of their parallel energy without any particle influx.

An idealized model presented in this paper reveals an interesting counterintuitive subtlety of electron heat losses. We show that the hot electrons contained between two cold walls can retain a significant part of their initial kinetic energy in a self-sustained collisionless phase space structure. This highlights an interesting physics connection between the need to understand the plasma–pellet interaction and the basic features of electron phase space dynamics and steady-states of the Vlasov–Poisson equations.

A hot plasma slab between the solid walls cools down and eventually decays. The energy lifetime of this plasma depends significantly on the properties of the walls. In a conventional scrape-off layer model (see, e.g., Ref. 7), the electron flow to the wall creates a sheath potential that restrains the electron flux and drives the ions to the wall to maintain ambipolarity.⁸ The losses of both, the heat and the particles, are then governed by the ion motion. If, however, the wall does not just absorb the hot incident electrons but also emits cold electrons to maintain charge neutrality, then the electron cooling can occur without any ion motion. Whenever an electron is treated as “cold” in this work, we mean that its energy is much smaller than the initial thermal energy of the plasma electrons.

We consider a one-dimensional slab of the collisionless plasma between the two walls, assuming that the ions are immobile and form a uniform background, the electrons are initially Maxwellian, and the distance between the walls is $L$⁠. The walls are perfectly conducting so that the electric field vanishes at the walls and freely emitting so the emitted current is always equal to the absorbed current.¹ The boundary condition we use is justified for the following reason. The pellet can be considered as a cold dense plasma or, in some cases, a metal. We, therefore, treat it as a perfectly conducting wall. The electric field should then vanish within the pellet volume up to the sharp pellet–plasma interface. To be more precise, the electric field needed to carry the return current of the pellet electrons is, of course, finite, but it is much lower than the electric field within the sheath adjacent to the pellet surface (the latter is of the order of $T/eλD$⁠).

The cooling process we attempt to describe is expected to take place on the electron time scale $L/vT$⁠. This suggests the following rough estimate for the rate of electron temperature decay:

Solving Eq. (1) with the initial temperature $T0$⁠, we obtain

Note that this simple estimate ignores an electric field that can build up during the cooling process. A more accurate approach to the problem requires a self-consistent solution of the kinetic equation for the electron distribution function $f(v,x,t)$ and the Poisson equation for the electric potential $φ(x,t)$⁠. It is instructive to compare Eq. (2) with such a solution. To do that, we use a PIC code. We assume nonrelativistic velocities of all particles, and, therefore, treat the plasma within the electrostatic approximation. In our numerical computations, we set the electric field to vanish on the surface of each wall. Each time an incident electron passes the wall surface, we replace it with a new one located in the vicinity of the wall. We set the initial velocity of the new injected electron to be negligible compared with the initial thermal velocity. We first choose $L=10λD$⁠, where $λD=T0/4πn0e2$ is the Debye length, $e$ is the electron charge, and $n0$ is a constant ion density. In Fig. 1, we plot the normalized values of the electron kinetic energy $wk$ and the electrostatic energy $wE$⁠, both averaged over the system length $L$

As seen in Fig. 1, the red curve follows the analytical estimate (black curve) precisely during the initial phase of cooling at which the space charge is still negligible, but the process of cooling virtually stops after approximately three electron plasma wave periods. The temperature of the wall is chosen to be very low in our simulations. To be precise, the initial velocity of the emitted electrons is set to be less than one percent of the initial thermal velocity. Consequently, the residual electron kinetic energy cannot be attributed to the electron-wall thermalization process. The observed counterintuitive behavior deserves a more detailed discussion that we focus on herein.

In Fig. 2, we compare snapshots of the electron distribution function at the initial phase of electron cooling and at the plateau in Fig. 1. Each point in Fig. 2 represents an individual electron. In blue, we show the emitted cold electrons. In red, we show the confined hot electrons, which do not reach a wall in their motion. As seen in the right panel of Fig. 2, there is a distinct phase space bucket filled with trapped particles during the plateau in Fig. 1. This structure maintains its integrity, albeit some small jitter with the electron plasma frequency time scale. The overall shape of the phase space structure remains unchanged during the plateau. This suggests the existence of a nonlinear steady-state, around which the system oscillates.

At a qualitative level, the following scenario sets the phase space structure in the right panel of Fig. 2. Initially, the hot Maxwellian electrons rapidly transfer their energy to the walls creating a secondary cold electron population in the vicinity of the walls (see the left panel of Fig. 2). Once the cold electron density near the walls becomes large enough for the potential well to form, the primary electrons of the Maxwellian bulk become restricted to the central region where they are trapped (see the right panel of Fig. 2 and Ref. 9). The structure shown in the right panel of Fig. 2 has a two-component electron population. One of the components represents cold beam electrons emitted from the walls with zero kinetic energy (see Refs. 10 and 11). These electrons accelerate within the potential drop near the wall and decelerate near the opposite wall. They do not transfer any energy to the walls if the potential does not evolve in time. The other component represents trapped electrons. Their energy spread can be comparable to the initial thermal energy. The system can thereby retain a significant fraction of its initial thermal energy despite being in contact with the cold walls. The sustainment of this unusual state is, however, constrained by stability requirement, which we will address after considering the steady-state itself.

To construct a steady-state solution, we consider the above-mentioned two-component electron distribution function: the cold beam electrons and the trapped electrons within the potential well. The electron distribution function satisfies the Vlasov kinetic equation

with a self-consistent potential determined by the Poisson equation

where

is the electron density.

The electric field has to vanish at the surfaces of the emitting walls. Since the potential is defined up to a constant and the walls are identical, we can choose the potential to vanish at the walls. We thus have the following boundary conditions for Eq. (5):

The total electron density in Eq. (5) is

where $nt$ is the density of the trapped electrons and $nb$ is the density of the cold electron beams. In a steady-state, the electron energy

is a constant of motion, and the electron distribution function $f$ depends only on $E$⁠. The trapped electrons are those with $E<0$⁠. Given a distribution function $f(E)$⁠, we rewrite Eq. (6) as

As seen in Fig. 2, the potential well traps predominantly low-energy particles for which the Maxwellian distribution is roughly (albeit not exactly) constant. To interpret our numerical results qualitatively, we assume that the distribution function of the trapped electrons is uniform over their energy range, i.e.,

where $β$ is a positive constant. We then obtain from Eq. (10) that

Note that the phase space bucket traps predominantly slow electrons from around the flattop of the initial Maxwellian distribution. It is then reasonable to choose the value of $β$ to match the flattop, which gives

Equation (13) specifies the value of beta for which the value of the top-hat distribution (11) matches the value of the Maxwellian distribution at $v=0$⁠. The distribution function of the cold beam is

where $α$ is a positive constant. The corresponding expression for the beam density is

We now substitute expressions (12) and (15) into the Poisson equation (5) and transform Eq. (5) to dimensionless variables $ξ=x/λD$ and $ψ=|e|φ/T0$ to obtain

Although the electric field vanishes at the walls due to Eq. (7) and the initial velocity of the emitted electrons is set to zero [see Eq. (14)], the emitted electrons leave the potential hill within a finite time. The diverging electron beam density dominates on the right-hand side of Eq. (16) in the vicinity of the wall. This is reminiscent of the classical problem solved by Langmuir.¹ Using the space charge limited potential obtained in Ref. 1 and the emitted electron energy conservation (9), one can show that $ξ(t)∝t3$ for the emitted electrons. Consequently, the electrons do leave the wall, although their local density diverges at the surface.

We next multiply Eq. (16) by $dψ/dξ$ and integrate it once using Eq. (7), which gives

where

with $Δ≡1−4αβ$ and $ψ01/2≡1/(2β)$⁠. Note that the polynomial $P2$ has two real roots when $Δ>0$

where $ψ−$ is the smallest root and $ψ+$ is unphysical (see Fig. 3).

An electric potential that holds trapped electrons must have a maximum between the two walls, which means that the two branches of Eq. (17) must be matched at the maximum of $ψ$ to construct a continuous nonmonotonic function $ψ(ξ)$ with $ψ=0$ at the walls. The lower panel of Fig. 3 shows a set of candidate single-maximum solutions for different values of $Δ$⁠. The value of $Δ$ can now be chosen to meet the requirement that $dψ/dξ$ vanishes at $ξ=±L/2λD$⁠. The symmetry of the problem suggests that the branches of a single-well electric potential have to merge at the geometrical center (⁠$ξ=0$⁠). Based on these remarks, the value of $Δ$ that guarantees smooth connection between the branches of (17) and ensures that $ψ=0$ at the walls is determined by the following relation:

where $ψ−$ is defined by Eq. (19). The same condition (20) also ensures that the total charge of the system is zero. For convenience, we herein present an analytical solution to the boundary value problem (16) and (7) in terms of the incomplete elliptic integrals of the first $F([ψ/ψ−]1/4|ℵ2)$ and the second kind $E([ψ/ψ−]1/4|ℵ2)$ (see Ref. 12 and  Appendix A)

where

Using Eqs. (21) and (22), we can rewrite Eq. (20) as follows:

where $K(ℵ2)$ and $E(ℵ2)$ are the complete elliptic integrals of the first and the second kind, respectively (see Ref. 10).

Equation (23) gives a one-to-one correspondence between the length of the system $L$ and the value of $Δ$ (or, equivalently, the value of $α$ for a given value of $β$⁠). As seen in Fig. 3, the larger the system length, the smaller the root $Δ$ of Eq. (23). Also, a system of larger length tends to have a larger amplitude of the potential $ψ$⁠. In the limit of $L→∞$⁠, the maximum value of the potential is $ψm=ψ0$⁠.

As the length of the system increases, the electric potential tends to be constant throughout most of the system length (see the $Δ=10−7$ case in the lower panel of Fig. 3). The electric field is then localized near the walls within several Debye lengths. Note that the average kinetic energy density of the phase space structure shown in the right panel of Fig. 2 can be roughly estimated as $wk≈n0T0$⁠. The average electrostatic energy density of a large system can be estimated as $wE≈n0T0λD/L$⁠. As the system length gets larger, the ratio $wE/wk∝λD/L$ decreases. These estimates are consistent with what is seen in Fig. 1 at the plateau phase of cooling.

At this point, we note extra freedom in constructing solutions to Eq. (16). First, we split the system length $L$ into $n$ equal segments, where $n$ is an integer. We can then merge the branches of the electric field (17) at the midpoints of each segment of length $L/n$⁠. This procedure guarantees continuity of the electric potential and its first derivative at each point along the system length. The electric potential constructed this way has $n$ maxima, i.e., $n$ potential wells within which the electrons are trapped. The global maximum of the n-well electric potential depends on the number of wells. As seen in the bottom panel of Fig. 3, the maximum value of the electric potential decreases with the system length. The right panel of Fig. 4 shows the one-well potential and a 4-well potential for a fixed system length $L=20λD$⁠.

As already pointed out, the global maximum of the n-well potential decreases with $n$⁠. As the amplitude of each individual well decreases, the velocity spread of the electrons trapped within each well decreases accordingly, and so does the averaged kinetic energy of the system. It is therefore conceivable that the single-well solution is actually metastable and can break into multiple wells. To examine this conjecture more accurately, we consider a dimensionless energy functional that includes both the electrostatic energy and the kinetic energy of the electrons (see  Appendix B)

where $Wn$ is the normalized energy of the n-well system. In the limit of $L→∞$ for a single-well system, we obtain

For the value of $β$ given by Eq. (13), the corresponding value of $W1(L→∞)=9π/128$ is one and a half times larger than the plateau in Fig. 1.

In Fig. 5, one can see that the system energy decreases with the number of wells. In the limit of $n→∞$⁠, the plasma is effectively cold. It is, however, interesting that the system does not necessarily relax to this lowest energy state.

As shown in Sec. III, the total energy of $n$ identical wells of length $L/n$ is less than the single-well energy. Splitting a single well into multiple wells could, therefore, release the energy difference. Feasibility of such splitting depends on whether the single well is linearly stable. The large length of the system allows us to use a local stability analysis since such a system is uniform throughout most of its length (see the lower panel of Fig. 3). We can, therefore, use a conventional¹³ dispersion relation for electron plasma

where $k$ is a wave number, $ω$—wave frequency, and $ωp=4πn0|e|2/m$—electron plasma frequency. Substitution of the steady-state electron distribution function [see Eqs. (11) and (14)] into Eq. (26) reduces Eq. (26) to

Equation (27) is quadratic in $ω2$ and has the following roots:

Note that both roots are positive so that the system is linearly stable.

The boundary conditions (7) and Eq. (13) determine the density ratio of the trapped electrons and the beam electrons. They also define the beam velocity, prescribing the entire shape of the electron distribution function (see the right panel of Fig. 6). It is instructive to generalize Eq. (27) to the case of an arbitrary density ratio of the trapped electrons and the beam electrons $nt/nb≡r/(1−r)$⁠, where $0≤r≤1$⁠, and an arbitrary value of the beam velocity $v0$

The case of $r=0$ is clearly unstable because it represents a two-stream distribution function. The opposite case of $r=1$ is stable. This reveals the stabilizing nature of the trapped electrons. The trapped electrons stabilize the two-stream instability when their density overcomes the beam electron density (⁠$nt≥nb$⁠) or, equivalently, $r≥1/2$ (see also the analysis of Ref. 10). In particular, all roots of Eq. (27) are stable because we have $nt=3nb$ and $r=3/4$ in the problem of our interest.

The linear stability of the single-well solution and the lower energies of the multiwell solutions show that the single-well solution is only metastable. The system will not necessarily relax to the single-well state if the initial distribution function differs significantly from the steady-state ansatz (11) and (14) (see also the left panel of Fig. 2). To illustrate that, we use a PIC simulation of the system with $L=20λD$⁠. The ions are immobile, and the initial electron distribution is Maxwellian.

In Fig. 7, we plot the electron distribution function after one, five, fifteen, and thirty electron plasma wave periods. In the top right panel of Fig. 7, one can see that the system splits into a chain of three buckets. The further system behavior comes down to rapid cooling of the large buckets. After 30 electron plasma wave periods, the system becomes essentially cold (see the bottom right panel of Fig. 7). In Fig. 8, as in Fig. 1, one can see that the red curve (averaged electron kinetic energy) initially follows the black curve [analytical estimate of temperature loss (2)]. The difference between the system of length $L=10λD$ and the system of length $L=20λD$ becomes evident after five electron plasma wave periods. After 30 electron plasma wave periods, the plasma is essentially cold, although it still contains some residual energy in the form of plasma oscillations.

On the other hand, a single-well solution should form when the initial condition lies within its convergence range. In order to demonstrate that, we choose the initial state that is qualitatively similar (see the left panel of Fig. 9) to the steady-state solution (11), (14) and (21), and let the system evolve. The initial ratio of the electron beam density to the trapped electron density in the geometrical center of the system is $nt/nb=3$⁠. We keep the system length fixed $L=20λD$⁠, as in the previous case.

In the left panel of Fig. 9, we plot the initial electron distribution that consists of the trapped electron population and two counterpropagating beams. We note that, after thirty electron plasma wave periods, the electron distribution function still preserves its initial structure, except for a small jitter. To have a better picture of how the system reaches the steady-state, we plot the time traces of the averaged electron kinetic energy $wk$ and the averaged electrostatic energy $wE$ (see Fig. 10).

As seen in Fig. 10, it takes approximately 10 electron plasma wave periods for the system to reach the steady-state. A violent initial jitter of the electron distribution function eventually settles down. The electron kinetic energy saturates and the electrons preserve a significant (>60%) fraction of the initial thermal energy. The “roughness” of the initial setup in this example leads to the electrostatic activity during the transition period within the relaxation time scale. This roughness suggests that the radius of convergence of the steady-state is large enough, so even a significantly perturbed initial state converges toward the steady-state.

Concluding our analysis we note that the steady-state we have constructed may have a finite lifetime due to the electron–electron collisions, the ion motion, or the 3D effects. An individual electron colliding with another electron can gain the kinetic energy sufficient enough to escape from the potential well and eventually penetrate into the wall. This process decreases the average kinetic energy of the system thus forces the plasma to become cold. The ion motion changes the shape of the potential well and, due to the conventional scrape-off layer model of Ref. 7, also cools down the plasma. Finally, the 3D effects may provide an instability, that the 1D model is not capable of taking into account, that possibly destroys the steady-state.

The presented model shows that the initially hot, collisionless plasma slab contained between two cold and freely emitting walls demonstrates rather counterintuitive behavior. The electron temperature does not decay as predicted by the simple arguments [see Eq. (2)] but rather saturates and stays constant on the electron time scale. We have found steady-state solutions for the electron distribution function that consists of the phase space buckets filled with trapped electrons, and two electron beams edging the trapped electron population. The potential, in this case, has an n-well structure. The linear stability analysis shows that an infinite-length phase space bucket is stable throughout its almost entire length. Our numerical simulations have demonstrated that the regions near the walls do not affect significantly the one-bucket system stability. It was also found that the initial conditions determine whether or not the initial electron distribution function will reach the one- or n-bucket steady-state, or whether it will eventually decay into a cold plasma.

The authors are grateful to the organizers and participants of the Festival de Théorie (2019) for their support and stimulating discussions. This work was supported by the U.S. Department of Energy Contract Nos. DEFG02-04ER54742 and DESC0016283.

We first integrate (17) from the left wall to $ξ$⁠, where $ξ≤0$

where Eq. (19) defines $ψ∓$⁠. Let $x/ψ−=z4$⁠, which transforms (A1) to

or, equivalently

where

Following,¹² we define the incomplete elliptic integrals of the first $F([ψ/ψ−]1/4|ℵ2)$ and the second $E([ψ/ψ−]1/4|ℵ2)$ kind:

With (A5), we have

where $−L/2λD≤ξ≤0$⁠. Going through the same steps but integrating (17) from the right wall to $ξ$⁠, we generalize (A6) to

In a steady-state, Eq. (4) reduces to

We multiply (B1) by $mv/2$ and integrate over velocities to obtain

where $w̃k≡1n0T0∫−∞+∞mv22fdv$ is a normalized electron kinetic energy density. We next use Eqs. (12) and (15) to express $ne$ in terms of $ψ$ and then integrate (B2) to find

This expression takes into account that $w̃k$ vanishes at the walls where $ψ=0$⁠.

We now add $w̃k(ξ)$ to the normalized electrostatic energy density $w̃E=0.5(dψ/dξ)2$ and average the result over the system length to find an average normalized value of the energy density in the n-bucket system

Note that the first term under the integral can be written as

where

[see Eq. (16)]. We also note that

because of the boundary conditions. One can, therefore, transform (B4) into