Thermal equilibrium of collisional non-neutral plasma in a magnetic dipole trap

Abstract

This paper discusses thermal equilibrium states of single-species plasmas, such as pure electron plasmas and pure positron plasmas, that are confined in a dipole trap. Thermal equilibrium states for such plasmas are routinely realized in the homogeneous magnetic field of Penning–Malmberg traps. We generalize the theory of these states to include inhomogeneous magnetic dipole fields. The approach to thermal equilibrium takes place in two stages with well separated time scales. On the collision time scale, thermal equilibrium is established along each magnetic field line. On the much longer transport time scale, heat conduction and viscosity bring the plasmas on different flux contours into thermal equilibrium, we call this a state of global thermal equilibrium. We present numerical results for local and global thermal equilibria. These results agree with the analytic predictions for charge collections that are large compared with the Debye length. There is, in principle, no limit to the confinement time of a single-species plasma in a global thermal equilibrium state. Experiments with hot electron–ion plasmas performed in the LDX and RT1 devices give us confidence that, in contrast to a Penning–Malmberg trap, a magnetic dipole field can also confine cold quasi-neutral electron–positron pair plasmas on the time scale of the phenomena of interest. Such pair plasmas are assumed to form in the magnetosphere of neutron stars but have so far not been realized in a laboratory. The creation of an electron–positron pair plasma is the main goal of the APEX collaboration.

1. Introduction

A dipole trap is a magnetic confinement configuration in which the magnetic field is produced by a simple circular coil or permanent magnet. Such traps were implemented with either a supported magnet or a levitated coil. The RT1 experiment confined a pure electron plasma in a dipole trap for times of order 300 seconds (Saitoh et al. 2010) as well as a hot electron–ion plasmas for 0.5 seconds (Yoshida et al. 2012). The APEX collaboration will confine quasi-neutral electron–positron plasmas in a dipole trap and test the numerous theoretical predictions for plasmas with perfect mass symmetry (Stenson et al. 2017). The planned injection scheme requires initial confinement of a pure electron or a pure positron plasma in the trap (Stoneking et al. 2020). It is the confinement of these single-species plasmas in a dipole trap that motivates the work presented here. The equilibrium distributions we find for such a non-neutral plasma resemble the magnetosphere of a neutron star (Pétri, Heyvaerts & Bonnazzola 2002). Even though these collections of charged particles are not quasi-neutral we consider them as plasma in the sense that they have a Debye length shorter than the system size and are therefore governed by collective behaviour.

In contrast to the work of Sato, Kasaoka & Yoshida (2015) we consider time scales longer than the collision time. On this time scale the plasma relaxes to thermal equilibrium locally along each magnetic field line (Hyatt, Driscoll & Malmberg 1987), and by azimuthal symmetry the local thermal equilibrium extends to the surface of constant magnetic flux. Such a local thermal equilibrium can be observed in a broader class of configurations (Pedersen & Boozer 2002) compared with a global thermal equilibrium state. On a much longer time scale, transport produces an interchange of particles and heat between flux surfaces and establishes a state of global thermal equilibrium for the whole plasma (Driscoll, Malmberg & Fine 1988).

States of global thermal equilibrium are routinely obtained for single-species plasmas in a Penning–Malmberg trap. For these traps, radial confinement of the plasma is provided by a uniform axial magnetic field and axial confinement by electrostatic fields at the ends of the trap. In the construction of these traps, great care is taken to ensure that the trap configuration is azimuthally symmetric. Consequently, the total canonical angular momentum is a constant of the motion. Both theory (Prasad & O'Neil 1979) and experiment (Driscoll et al. 1988) have shown that plasmas with a single sign of charge may come to global thermal equilibrium in a Penning–Malmberg trap. Such a plasma is in a state of maximum entropy subject to fixed values of the total plasma energy, the total plasma canonical angular momentum and the total particle number (Dubin & O'Neil 1999). A maximum-entropy state is guaranteed to be stable and quiescent and, in principle, to persist indefinitely. In practice asymmetries (Kabantsev et al. 2003) and collisions with residual gas (Chao, Davidson & Paul 1999) apply a torque to the plasma and lead to a slow loss of the plasma to the wall.

The magnetic-field configuration of the dipole trap inherits the azimuthal symmetry of the circular coil, so it is natural to ask if the dipole trap admits confined thermal equilibrium states for a single-species plasma. To understand the confinement, we note that in the purely poloidal magnetic field B the drift motion results in a rotation of the plasma around the central axis of the dipole trap. In addition to the diamagnetic drift, the inherent electric field E of a non-neutral results in a significant $E\times B$ drift. The consequent inward Lorentz force balances the outward forces due to pressure and space charge. Alternatively, one can say that the magnetic field induces an electric potential in the rotating rest frame of the plasma, and when the sum of the induced potential and the electrostatic trap potential possesses a local minimum, the plasma can come to global thermal equilibrium in that potential well. For a global thermal equilibrium state the plasma rotation must be shear free; otherwise viscosity acting on the shear would produce entropy, which is not possible for a state of maximum entropy.

Regarding the confinement of a quasi-neutral plasma one needs to consider that, both the induced and the electric potentials depend on the sign of the charge. What looks like a local minimum for one particle species is therefore a local maximum for a species with opposite sign of charge. In addition, the $E\times B$ drift is absent in a quasi-neutral plasma and the diamagnetic drifts point in opposite directions for positive and negative species. Hence, it is not possible to confine a quasi-neutral plasma in a state of global thermal equilibrium.

For the case of a bare coil, alone in unbounded space, we will see that plasma rotation in the magnetic field does indeed produce a deep potential well. However, the coil is located at the bottom of the well, leading to particle loss onto its surface. Consequently, the coil is expected to charge up until its electric potential prevents further flow of particles to the coil. However, when the charge on the coil is sufficiently large to suppress further charge accumulation, the local minimum in the potential vanishes. In experiments (Saitoh et al. 2010; Stoneking et al. 2020), a levitated coil is surrounded by conducting walls that bound the confinement region and may be electrically biased. As we show below, one can specify electric potentials on sections of the walls that establish a minimum in the effective potential away from material surfaces. These wells tend to be shallower than the well with the coil at the bottom, but they offer experimentally realizable confinement scenarios.

To find specific thermal equilibrium states we solve Poisson's equation for the self-consistent charge density and electric potential, where the functional form of the charge density is given by the thermal equilibrium distribution function. This task is non-trivial since the electric potential enters the distribution function nonlinearly. Using a numerical iteration scheme (Spencer, Rasband & Vanfleet 1993), we solve for several global thermal equilibrium states, including low-temperature states. For these states, considerations involving Debye shielding and zero-pressure-force balance yield simple analytic expressions for the local plasma density in terms of the local magnetic field (Goldreich & Julian 1969; Turner 1991). These expressions are verified by our numerical solutions.

The class of local thermal equilibrium states is much larger than that of global thermal equilibrium states and does not require the existence of a potential well. Global thermal equilibrium states are parameterized by the values of three constants, whereas local thermal equilibrium states are parameterized by an infinite number of constants, that is, the number of particles and the temperature on each magnetic flux contour. Every global thermal equilibrium state is also a local thermal equilibrium state, but the converse is not true.

The theory of local thermal equilibrium states is not as robust as the theory for the global thermal equilibrium states, and there is no guarantee for the stability of these states. The basic assumption underlying the theory of local thermal equilibrium is that the transport time scale is much larger than the collisional time scale. The plasma comes to local thermal equilibrium on the collisional time scale and this state lasts for the transport time scale. This separation in time scales implicitly assumes that each particle is well localized to its own flux surface, that is, that the cyclotron radius is small compared with spatial scale of the magnetic field gradient. Consequently, one expects a guiding-centre drift dynamics, or single-fluid Magnetohydrodynamics (MHD) to provide a valid description of the dynamics. We will characterize the allowed local thermal equilibrium states by requiring that they satisfy force balance using the single-fluid MHD model.

The paper is organized as follows. Section 2 discusses general theoretical considerations regarding global thermal equilibrium states for a non-neutral plasma in a dipole trap. Section 3 describes our strategy to identify promising vacuum-state conditions for realizing a state of confined global thermal equilibrium in the dipole field; this involves looking for configurations with isolated local minima in the effective potential seen by a particle in a rotating reference frame. In § 4 we present numerical solutions for global thermal equilibrium states that make use of tailored electrostatic boundary conditions. Finally, in §§ 5 and 6 we discuss the theory of local thermal equilibrium states and present numerical results for such solutions in a dipole field.

2. Theory of global thermal equilibrium in a dipole field

For a real dipole trap experiment, such as the RT1 experiment (Ogawa et al. 2009) or the planned APEX (A Positron-Electron eXperiment) (Stoneking et al. 2020), the dipole field is produced by a coil bounded by a conductor with a finite cross-section. For analytic simplicity we take the coil to be a circular wire of radius $R$, with current $I$ taken to be in the $-\varTheta$ direction. Also, we neglect the field used to levitate the coil. We let the current in the levitated coil lie in the $z = 0$ plane of a cylindrical ($r, z, \varTheta$) coordinate system with the centre of the coil at $r = 0$, as shown in figure 1.

Figure 1. Simplified model of a levitated dipole trap. The cylindrical ($r, z, \varTheta$) coordinate system is indicated by the black arrows. The red solid lines indicate contours of constant magnetic flux $\alpha$. Due to the symmetry of the trap it is sufficient to consider one quadrant of the trap. The following two-dimensional (2-D) diagrams will only show the upper right quadrant.

Maximizing the plasma entropy subject to the constancy of total energy, total canonical angular momentum and total particle number yields the distribution function (Dubin & O'Neil 1999)

(2.1)\begin{equation} f = C \exp\left[-\frac{1}{T}(H-\omega p_{\varTheta})\right], \end{equation}

where $C$ is a normalization constant, $T$ is the plasma temperature and $\omega$ is the plasma rotation frequency. The quantity

(2.2)\begin{equation} H = \frac{m}{2}[\dot{z}^2 + \dot{r}^2 + r^2\dot{\varTheta}^2] + q\phi(r, z), \end{equation}

is the Hamiltonian, or energy, for a charged particle, and the quantity

(2.3)\begin{equation} p_{\varTheta} = m r^2\dot{\varTheta} + \frac{q}{{c}}rA_{\varTheta}(r, z), \end{equation}

is its canonical angular momentum. Here, $q$ is the particle charge, $m$ is the particle mass, c is the speed of light, $\phi (r, z)$ is the electric potential and $A_{\varTheta }(r, z)$ is the azimuthal component of the vector potential. The combination

(2.4)\begin{equation} H - \omega p_{\varTheta} = \frac{m}{2}(\dot{z}^2 + \dot{r}^2 + r^2(\dot{\varTheta}-\omega)^2) + q\phi(r, z) - \omega\frac{q}{{c}}rA_{\varTheta}(r, z) - \frac{m}{2} r^2\omega^2, \end{equation}

is the particle energy in a frame that rotates with angular frequency $\omega$. The last term is the centrifugal potential energy and the penultimate term is the potential energy associated with the electric field that is induced by rotating through the magnetic field. We call this term the induced potential, $\phi_{ind}(r, z) = -\left(\omega r/c \right) A_{\Theta}(r, z)$. Note that a change in the sign of charge changes the sign of the electrostatic potential as well as that of the rotation frequency and thereby the sign of the induced potential.

In the case of an azimuthal symmetric, poloidal magnetic field, the vector potential can be written in terms of the axial magnetic flux per unit angle $\alpha$

(2.5)

where the contour integral is around a circle through some point ($r, z$) with its centre on the axis and lying in a plane of constant $z$. The surface integral is over the disc bounded by the circle. The effective trap potential in the rotating frame is given by the expression

(2.6)\begin{equation} q\phi_{{\rm rot}}(r, z) = q\phi_{{\rm trap}}(r, z) - \frac{\omega q}{{c}} \alpha(r, z) - \frac{m}{2} \omega^2r^2. \end{equation}

The trap potential $\phi _{{\rm trap}}$ in the laboratory frame is due to the electrostatic potential on the coil and the wall. We will be interested in cases where $q\phi _{{\rm rot}}(r, z)$ is a local potential well.

Equations (2.1) and (2.4) imply that the density is described by a Boltzmann distribution

(2.7)\begin{equation} n(r, z) = N \frac{\exp\left[-\dfrac{1}{T}(q\phi_{{\rm rot}} + q\phi_{p})\right]}{\displaystyle \int {\rm d} V \exp\left[-\frac{1}{T}(q\phi_{{\rm rot}} + q\phi_{p})\right]}, \end{equation}

where $N$ is the total number of particles in the trap and $\phi _{p}$ is the plasma potential. The normalization integral in the denominator of (2.7) extends only over the region of the potential well. This form of the density distribution is specified by the three thermodynamic variables: $T$, $\omega$ and $N$. Alternatively, one can replace the number of particles $N$ by the mean density within the well $n_0$

(2.8)\begin{equation} n_0 = \frac{N}{\displaystyle \int {\rm d} V \exp\left[-\frac{1}{T}(q\phi_{{\rm rot}} + q\phi_{p})\right]}. \end{equation}

Both forms are useful, but they describe different classes of equilibria. For the plasma to be confined in a state of global thermal equilibrium, it is necessary that $q\phi _{{\rm rot}}(r, z)$ possess a potential well that is substantially deeper than the temperature $T$ and is only partially filled such that near the edge of the well the density falls off exponentially.

Taking the gradient of the density distribution (2.7) yields the statement of force balance

(2.9)\begin{equation} T \boldsymbol{\nabla} n + n\boldsymbol{\nabla}[q\phi_{{\rm rot}} + q\phi_{p}] = T \boldsymbol{\nabla} n + n\left[q\boldsymbol{\nabla}\phi - \frac{q}{c} \omega r B \boldsymbol{\hat{\alpha}} - m\omega^2r \boldsymbol{\hat{r}}\right] = 0. \end{equation}

Use has been made of the relation $\boldsymbol {\nabla }\alpha = r B \hat {\alpha }$. Here, the unit vector $\hat {\alpha }$ is orthogonal to the flux surface $\alpha$, pointing in the direction of increasing flux. Taking the dot product of (2.9) with the unit vector $\hat {\alpha }$ yields the relation

(2.10)\begin{equation} \omega = \frac{{c}\hat{\boldsymbol{\alpha}}\boldsymbol{\cdot}\boldsymbol{\nabla}\phi}{rB} + \frac{{c}T \hat{\boldsymbol{\alpha}} \boldsymbol{\cdot} \boldsymbol{\nabla} n}{qnrB} - \frac{{c}m\omega^2\hat{\boldsymbol{\alpha}}\boldsymbol{\cdot}\hat{\boldsymbol{r}}}{qB}. \end{equation}

The theory of global thermal equilibrium does not assume that the particle motion is well approximated by guiding-centre drift theory, but when that is the case, (2.10) can be interpreted in terms of drifts. The first term in (2.10) is the $E \times B$ drift, the second is the sum of the diamagnetic, grad $B$ and curvature drifts and the last term is the polarization drift driven by the centrifugal force. This interpretation will be discussed more thoroughly in § 5 where we treat local thermal equilibrium states. For typical experimental conditions, the centrifugal term tends to be small. Quantitatively, the centrifugal term is small compared with the $E \times B$ term when $\omega \ll \omega _{cz} = qB_z/m{c}$ where $\omega _{cz}$ is the cyclotron frequency evaluated for the $z$- component of the magnetic field. This will be apparent in the expression below (2.11) for the zero-temperature limit for the density.

The name plasma is reserved for charge collections that are large compared with the Debye length. We consider this low-temperature/high-density regime from the perspective of local force balance. Let the plasma potential be written as an expansion in temperature $\phi _p \approx \phi _p^{(0)} + \phi _p^{(1)}$. To zeroth order, the first term in (2.9) may be neglected to give $n_{{\rm cold}} \boldsymbol {\nabla } (\phi _p^{(0)} + \phi _{{\rm rot}}) = 0$. Thus, at any point in space either $n_{{\rm cold}} = 0$ or $\boldsymbol {\nabla }(\phi _p^{(0)} + \phi _{{\rm rot}}) = 0$. Assuming that the density is not zero at some point in space, taking the divergence of the gradient and using Poisson's equation yields an expression for the zero-temperature density

(2.11)\begin{equation} 4{\rm \pi} qn_{{\rm cold}}(r, z) ={-}\nabla^2 \phi_p^{(0)} = \nabla^2 \phi_{{\rm rot}} ={-}\frac{2m}{q}\omega( \omega_{cz} + \omega), \end{equation}

which was found by Goldreich & Julian (1969) and Turner (1991). The Laplacian acts on all three terms in (2.6) for $\phi _{{\rm rot}}$. The first term $\nabla ^2 \phi _{{\rm trap}}(r, z)$ vanishes in the vacuum region. The second term yields $\nabla ^2\alpha = 2B_z$ when $\boldsymbol {\nabla } \times \boldsymbol {B} = 0$. For these plasmas, the diamagnetic modification to the magnetic field is negligible, even though the plasma is rotating. Since $qB_z/cm = \omega _{cz}$ is spatially dependent, the cold density distribution is also spatially dependent. The Laplacian acting on the third term is simply $2m\omega ^2$. The particle density $n_{{\rm cold}} (r,z)$ must be positive. According to (2.11), this condition requires that $-\omega (\omega _{cz} + \omega )$ be positive throughout the region occupied by the plasma.

From this perspective, based on force balance, the role of Debye shielding is not apparent. When the Debye length is small compared with the system size and the scale on which the magnetic field varies, the plasma shields out the total effective potential as seen in the rotating frame. The potential $\phi _{{\rm rot}}$ appears to arise from an effective charge density $4{\rm \pi} q^{\prime }n_{{\rm rot}}(r, z) = -\nabla ^2 \phi _{{\rm rot}}$ with $q = -q^{\prime }$. The plasma density adjusts itself to cancel out this effective charge density everywhere inside the plasma. In contrast to the actual plasma density, which is well separated from the coil and the walls, $n_{{\rm rot}}$ extends throughout the entire trap. Hence, a self-consistent boundary condition at the edge of the plasma needs to be satisfied. This boundary condition is given by the fact that, in the presence of the cold plasma, the total effective potential is constant $\phi _{{\rm rot}} + \phi ^{(0)} = {\rm const.}$ up to the plasma edge. Consequently, the plasma has to be enclosed by a contour of constant total effective potential. Within this potential contour the charge density is $n_{{\rm cold}} = n_{{\rm rot}}$, and it is zero outside the contour. Since the plasma potential contributes to the total potential this contour needs to be determined using an iterative scheme as described in § 4.

The second term in the expansion $\phi _p^{(1)}$ balances the pressure gradient in (2.9) to first order in temperature. We thus find that $q\phi _p^{(1)} = -T ln(n_{{\rm cold}}/n_0)$, which depends on position only through the spatial dependence of the magnetic field. At the edge, which extends a Debye length into the plasma, the approximation $\phi _p \approx \phi _p^{(0)} + \phi _p^{(1)}$ fails because $n_{{\rm cold}} \neq n_{{\rm warm}}$. In contrast to the sharp edge of the cold density distribution, the finite-temperature solution falls off on the scale of the Debye length (Prasad & O'Neil 1979). In this region it becomes necessary to determine the exact finite-temperature plasma density $n_{{\rm warm}}$. This density cannot be calculated analytically. In § 4 we present a numerical method to find this solution.

3. Identifying configurations with effective potential minima

We look for situations where the effective trap potential in the rotating frame, as given by (2.6), is a potential well. For simplicity we start by considering a bare circular coil in free space and neglect the thickness of the coil. Let the coil carry a current $I$ and charge $Q$. We scale the magnetic flux function for the dipole field and the trap potential as follows:

(3.1)\begin{equation} -\frac{\omega}{{c}} \alpha(r, z) + \phi_{{\rm trap}}(r, z) ={-}\frac{\omega I R}{{c}^2}[\tilde{\alpha}(\rho, \zeta) + b\tilde{\phi}_{{\rm trap}}(\rho, \zeta))], \end{equation}

where $\rho = r/R$ and $\zeta = z/R$ are scaled by the radius of the coil $R$ and $b = -{Q{c}^2}/{2{\rm \pi} R^2 \omega I}$ is a dimensionless factor. The dimensionless flux function $\tilde {\alpha }(\rho, \zeta )$ and trap potential $\tilde {\phi }_{{\rm trap}}$ can be calculated analytically using complete elliptic integrals (Jackson 1999). They depend only on geometry; the dependence on the current $I$ and charge $Q$ only appears in the scaling factors. As shown in the following section, the centrifugal term in (2.6) can be neglected for the rotation frequencies applicable for APEX.

Figure 2 shows plots of the term in the square brackets in (3.1) as a function of $\rho$ for $\zeta = 0$. The different curves correspond to different values of the parameter $0 \leq b \leq 0.8$. If the radius of the coil and its current as well as the rotation frequency do not change, $b$ is a measure of the charge on the coil. The bottom curve, corresponding to $b = 0$ for an uncharged coil, shows that the induced potential alone produces a deep potential well with the coil at the bottom. The coil in the APEX configuration is floating and would charge up as particles accumulate on its surface. Here, we ask if the charged coil can establish a potential barrier separating the potential well from the surface. For no value of $b$ in figure 2 is there a potential well that is separated from the coil by a potential barrier. For $b \geq 0.6$ a potential barrier starts to form, but at the same time the potential well vanishes. Therefore, the bare coil does not admit global thermal equilibrium states.

Figure 2. Plot of $\tilde {\alpha }(\rho, \zeta ) + b\tilde {\phi }_{{\rm trap}}(\rho, \zeta ))$ versus $\rho$ for the values $b = 0.0, 0.2,\ldots 0.8$. As $b$ increases to 0.6, a potential barrier is formed, preventing further flow of charge to the coil, but confinement at large $\rho$ is lost in the process.

However, real experimental configurations do not consist of a bare coil in infinite space; the confinement region is bounded by azimuthal symmetric conducting walls. To create a potential well that is separated from the coil by a potential barrier, we need to provide an additional inward radial force at large $\rho$. This can either be achieved with tailored electric potentials or through a distortion of the magnetic dipole field. One example is displayed in figure 3. It relies on grounded walls at $\rho = 3$ and $\zeta = \pm 1$. The walls are represented by a numerical solution to Laplace's equation that cancels the value of the coil's potential at the boundary. We assume that the conductors are not ferromagnetic and do not significantly modify the magnetic flux function. Implementing these boundary conditions yields the grey curve in figure 3 plotted for the value $b=0.8$. The proximity of the grounded top and bottom surfaces effectively shorts out the radial electric field at large $\rho$ while the field near the coil is relatively unchanged. In this case, there is a well with the bottom near $\rho =2$. The 2-D contour plot in figure 4 verifies that the minimum in the effective potential along the midplane is, in fact, a minimum in 3-D space rather than a saddle point. We note that particles with the opposite sign of charge could be confined at the centre of the coil. This region resembles a traditional Penning–Malmberg trap. The simultaneous confinement of two species with opposite signs of charge in the respective potential well agrees with the model of the magnetosphere of a neutron star (Pétri et al. 2002). In this model the pair plasma that originates from the surface of the neutron star is separated into a disc in the equatorial plane and a dome beyond the poles. The potential well on the outboard side is shallow compared with the one in the centre but it covers a larger volume. We find, that the central well can confine 2.5 times the number of particles compared with the outboard well. Details of the respective plasma are given in § 4. For a given trap geometry and fixed coil current the number of particles that can be confined scales linear with the rotation frequency as long as the charge on the coil is scaled by the same factor. Only at very high rotation frequencies (compared with what is relevant for APEX) does the quadratic term in (2.6) due to the mechanical part of the canonical angular momentum become significant and eventually destroy the potential well. For a Penning–Malmberg trap this effect is known as the Brillouin density limit. At the end of § 4 we provide parameters relevant for APEX as well as at the density limit.

Figure 3. The scaled effective potential, $\tilde {\alpha }(\rho, \zeta ) + b\tilde {\phi }_{{\rm trap}}(\rho, \zeta ))$ along the midplane for a value of the dimensionless charge on the coil, $b = 0.8$. The red line represents the coil in free space while the grey line includes the conducting boundary and indicates a shallow potential well.

Figure 4. The solid contour lines represent the effective vacuum potential in the rotating frame. The dotted grey lines indicate magnetic flux surfaces.

4. Numerical solutions for global thermal equilibrium states

In this section we find solutions to Poisson's equation

(4.1)\begin{equation} \left(\frac{1}{r}\frac{\partial}{\partial r} r \frac{\partial}{\partial r}+ \frac{\partial^2}{\partial z^2}\right) \phi_{p}(r, z) ={-}4{\rm \pi} q n_0 \exp \left[-\frac{1}{T}(q\phi_{{\rm rot}} + q\phi_{p})\right], \end{equation}

for global thermal equilibrium plasmas confined in the potential well of figure 4. We introduce the scaled potential $\psi = q\phi /T$ and express the density in terms of the scaled Debye length ${\lambda _{D0}^2}/{R^2} = {T}/{4{\rm \pi} q^2 n_0 R^2}$. This allows us to rewrite Poisson's equation in a dimensionless fashion

(4.2)\begin{equation} \left(\frac{1}{\rho}\frac{\partial}{\partial \rho} \rho \frac{\partial}{\partial \rho}+ \frac{\partial^2}{\partial \zeta^2}\right) \psi_{p}(\rho, \zeta) ={-}\frac{R^2}{\lambda_{D0}^2} \exp [-\psi_{{\rm rot}} - \psi_{p}]. \end{equation}

Solving (4.2) for the plasma potential $\psi _p(\rho, \zeta )$ is non-trivial since the potential enters the density distribution nonlinearly which requires a numerical method. The calculations are carried out on a rectangular grid that covers one quadrant of the trap with 3000 radial and 1000 longitudinal grid points. The resolution is increased by two orders of magnitude in the confinement region, from $10^{-2}R_{{\rm coil}}$ to $10^{-4}R_{{\rm coil}}$. We use an iterative approach similar to that used to find equilibrium states in Penning–Malmberg traps (Spencer et al. 1993). In each iteration, Poisson's equation is solved using a finite difference method with respect to the density distribution from the previous step, starting with the vacuum case. Due to the symmetry of the trap we can use $\partial \phi(\rho, \zeta=0)/\partial \zeta=\partial \phi(\rho=0, \zeta)/\partial \rho = 0$ as boundary condition at the inner edges of the quadrant. At the outer edges, which align with the grounded vacuum chamber, the potential is set to zero. The plasma potential is updated with the weighting factor $\gamma$ yielding the updated potential $\psi _{i+1} = \gamma \psi _p(n_i) + (1 - \gamma )\psi _i$ according to (4.2). The updated density is calculated using the updated potential. The density distribution is evaluated differently for a finite-temperature plasma and the zero-temperature limit. In the finite-temperature case, the density distribution is calculated from (2.7) and is set to zero everywhere outside the potential well. The boundary of the potential well is identified as the last closed contour for the total effective potential that does not enclose or intersect a material object. This contour needs to be evaluated for the respective potential of each iteration. For a finite-temperature solution in a deep enough potential well, the density at the boundary is exponentially small but not zero. For the zero-temperature limit, the density according to (2.11) does not explicitly depend on the electrostatic potential. We still need to employ an iterative scheme to find a self-consistent potential contour that determines the edge of the plasma. The potential contour defined for the finite-temperature solution would now result in the maximum amount of particles that can be confined. One might also choose a potential contour that encloses a subdomain of the potential well. That is equivalent to partially filling the trap with a cold plasma. The convergence of the solution is validated by taking the difference between the left and the right-hand side of Poisson's equation. The volume integral over the deviation from Poisson's equation can be related to a number of particles.

(4.3)\begin{equation} N_{{\rm error}} = \int {\rm d} V \left|n + \frac{1}{4{\rm \pi} q} \nabla^2 \phi\right|. \end{equation}

We use a dynamic weighting factor for faster convergence. If $N_{{\rm error}}$ decreases compared with the previous iteration the weighting factor increases and vice versa. The distribution with finite temperature displayed below converged at $N_{{\rm error}}/N_{{\rm tot}} \approx 10^{-9}$. For solutions with low temperatures or high densities it can be useful to start with a lower density or hotter plasma and gradually adjust the input parameters. In order to avoid rounding errors when evaluating the exponential we subtract the average of the total potential within the potential well from the argument in the exponent

(4.4)\begin{equation} n(\rho, \zeta) = n_0^\prime \exp [-\psi_{{\rm rot}} - \psi_{p} - \langle \psi_{{\rm rot}} + \psi_{p} \rangle_{{\rm well}}]. \end{equation}

The normalization factor is adjusted accordingly to $n_0^\prime = n_0 \exp [\langle \psi _{{\rm rot}} + \psi _{p} \rangle _{{\rm well}} ]$.

The example of a zero-temperature solution given here converged at $N_{{\rm error}}/N_{{\rm tot}} = 3 * 10^{-3}$. Since this solution is subject to finding a self-consistent plasma edge we assume that this method is more sensitive to the resolution of the grid compared with the finite-temperature case. Accordingly, the error of the zero-temperature solution doubles if we decrease the resolution by one half (1500 radial and 500 longitudinal grid points), while it does not change significantly for the finite-temperature solution.

Illustrative results of these calculations are shown in figure 5. It shows colour density contours for four plasmas of decreasing temperature from top to bottom and left to right. In each panel, the white contour shows the boundary of the potential well. For finite temperatures a solution for (4.2) was found for a given temperature, number of particles and rotation frequency. In contrast to the finite-temperature results, the lower right panel displays the density according to the zero temperature limit (2.11). For this case, density is cut off at a contour of constant total potential. The potential contour is chosen such that the well in the presence of the plasma has a depth of 10 % compared with the vacuum potential well. This constraint determines the total number of particles. The density distribution of the warm solutions are normalized to match the number of particles of the cold solution. Introducing a cold plasma in the same way but in the centre of the coil, which resembles a Penning–Malmberg trap, results in 2.5 times the number of particles compared with the plasma on the outboard side.

Figure 5. Density distribution for global thermal equilibrium solutions with different temperatures. The temperature is expressed in terms of the Debye length $\lambda _{D0}$ with respect to the mean density within the potential well $n_0$ and the coil radius $R$.

Next, we confirm that these results comply with the characteristics of a global thermal equilibrium state, force balance (2.9) and rigid rotation (2.10). In order to verify the statement of force balance we distinguish between the zero-temperature limit and a finite-temperature plasma. The blue curve in figure 6(a) shows the sum of the effective vacuum potential in the rotating frame $\phi _{{\rm rot}}$ and the plasma potential in the zero-temperature limit $\phi _p^{(0)}$. The red dashed curve in figure 6(a) is the sum of the effective vacuum potential, the plasma potential $\phi _p \approx \phi _p^{(0)}+\phi _p^{(1)}$ for the finite temperature solution with $\lambda _{D0}/R = 0.04$ and the pressure term to first order in temperature $T\ln (n_{{\rm cold}}/n_0)$ (only defined for $\ln (n_{{\rm cold}}>0)$). There is no gradient visible in the presents of the plasma for both cases. We can therefore say that the total potential is constant in the zero-temperature limit. For finite temperatures the variation in the plasma potential balances the pressure term. The density profile for different temperatures, including the zero-temperature limit is given in figure 6(b). In the zero-temperature limit the density is cut off at an equipotential contour (blue solid curve). The entire effective charge density $n_{{\rm rot}}$ is represented by the dashed blue line. As the temperature increases (blue to red curves) the edge of the plasma becomes smoother and falls off on the scale of a Debye length.

Figure 6. Values for the potential (a) and the density (b) in the equatorial plane ($\zeta = 0$). The sum of the effective potential in the rotating frame and the plasma potential in the zero-temperature limit (blue solid line) as well as the sum of the effective potential, the plasma potential of the finite temperatures and the pressure term (red solid lines) is displayed in (a). The density for different temperatures (solid line) as well as the complete solution of (2.11) (dashed line) are displayed in (b).

Maximum entropy implies the absence of rotational shear. More explicitly, the sum of all drifts matches the rotation frequency of the frame of reference $\omega _0$. In figure 7, the blue curve is the rotation frequency due to the $E\times B$ drift and the red curve is the rotation frequency due to the sum of the diamagnetic drift, grad $B$ drift and curvature drift in the equatorial plane. The purple curve is the sum of the two, which does in fact provide the unsheared rotation expected for global thermal equilibrium and equals the rotation of the frame of reference. The $E\times B$ drift is calculated from the potential that is obtained by solving Poisson's equation. Therefore, the constancy of the rotation frequency verifies that the solution to Poisson's equation yields the same potential as enters the density distribution. This is a self-consistency check of the solution.

Figure 7. The $E\times B$ rotation (blue line), the diamagnetic, grad $B$ and curvature rotation (red line) as well as the sum of the two (purple line) are evaluated in the equatorial plane ($\zeta = 0$). All frequencies are scaled by the rotation frequency of the frame of reference $\omega _0$.

The minimum target value of $10^{10}$ particles in the APEX trap would require a rotation frequency of the order of $10^6$ rad s$^{-1}$. For the warmest solution in figure 5 with $\lambda _{D0}/R = 0.08$ this would corresponds to a temperature of $T = 0.5$ eV. The ratio of the plasma rotation frequency and the cyclotron frequency associated with the $z$-component of the magnetic field is $\omega /\omega _{cz} < 4*10^{-7}$, meaning that the centrifugal term is negligible. The ratio between the total plasma charge and the charge on the coil is $Q_{{\rm plasma}}/Q_{{\rm coil}} = 0.74$ which corresponds to a potential of $\phi _{{\rm coil}} = 1.5$ kV. For this configuration we find the maximum number of particles at a rotation frequency of $3 * 10^8$ rad s$^{-1}$ with $2 * 10^{12}$ particles. However, this implies a coil potential of $\phi _{{\rm coil}} = 300$ kV. Increasing the rotation frequency further leads to a deformation of the potential well due to the effective centrifugal potential and thereby reduces the number of particles that can be confined. We chose this configuration because of its simplicity and did not optimize it for the maximum number of particles it can confine.

5. Theory of local thermal equilibrium in a dipole field

We can identify local thermal equilibrium states as separate entities only when the cross-field transport time scale is much longer than the collisional time scale. The state of local thermal equilibrium is created on the collision time scale (Hyatt et al. 1987) and lasts for the transport time scale (Driscoll et al. 1988). Collisions drive the velocity distribution along each magnetic field line towards an isotropic Maxwellian with a drift velocity. The drift velocity and the density distribution along the field line is determined by force balance. Due to the azimuthal symmetry this state of local thermal equilibrium extends to a 2-D surface of constant magnetic flux.

The separation of time scales implicitly assumes that each particle is well localized on its flux contour, that is, that the characteristic cyclotron radius is small compared with the scale length of the magnetic field gradient. For such a plasma single-fluid theory provides a good description. A discussion of the theory for a non-neutral plasma confined on the magnetic surfaces of a stellarator configuration was conducted by Pedersen & Boozer (2002).

In Sato et al. (2015) the authors consider equilibrium states for time intervals shorter than the collision time scale. They maximize the entropy while holding adiabatic invariants, such as the total cyclotron action and bounce action, constant. This strikes us as problematic since entropy is maximized dynamically by collisions but the adiabatic invariants are not conserved by collisions. In contrast to our work, the velocity distribution in Sato et al. (2015) can be anisotropic. The anisotropic temperature relaxation rate for an electron plasma in a Penning–Malmberg trap was measured by Hyatt et al. (1987) and is in agreement with collisional Fokker–Planck theory. The states obtained by Sato et al. (2015) are interesting equilibria, but they are not singled out as special by the collisional dynamics.

It is useful to introduce a set of orthogonal coordinates based on the magnetic field lines. The first of these coordinates is the magnetic flux function $\alpha (r, z)$ that was introduced earlier, and the second is the magnetic potential $\chi (r, z)$, defined through the relation $\boldsymbol {B} = \boldsymbol {\nabla } \chi$. The existence of the magnetic potential assumes that $\boldsymbol {\nabla } \times \boldsymbol {B} = 0$, which is consistent with a negligible diamagnetic current. The third coordinate is the azimuthal angle $\varTheta$. These orthogonal coordinates ($\alpha, \chi, \varTheta$) were introduced by Taylor (1964). With the corresponding infinitesimal line elements ${\rm d} s_{\chi } = {\rm d}\chi /B$, ${\rm d} s_{\alpha } = {\rm d} \alpha /rB$ and ${\rm d} s_{\varTheta } = r\,{\rm d} \varTheta$ we find the Laplacian

(5.1)\begin{equation} \nabla^2 = B^2(\alpha, \chi) \left(\frac{\partial}{\partial \alpha} r^2(\alpha, \chi) \frac{\partial}{\partial \alpha} + \frac{\partial^2}{\partial \chi^2}\right) + \frac{1}{r^2(\alpha, \chi)}\frac{\partial^2}{\partial \varTheta^2}. \end{equation}

For a plasma in local thermal equilibrium, the particle distribution is a local Maxwellian characterized by a local density $n(\alpha, \chi )$, local azimuthal drift velocity $v_{\varTheta } = \omega (\alpha, \chi )r(\alpha, \chi )$ and local temperature $T(\alpha )$. There is no dependence in these functions on the azimuthal angle $\varTheta$ since the system is assumed to have azimuthal symmetry. The temperature $T(\alpha )$ is independent of $\chi$ because collisions have established local thermal equilibrium along each field line. Associated with the density distribution is the self-consistent electric potential $\phi (\alpha, \chi )$.

Relationships between the functions $n(\alpha, \chi )$, $\omega (\alpha, \chi )$, $T(\alpha )$ and $\phi (\alpha, \chi )$ are established by the requirement that the equilibrium satisfy the time-independent single-species fluid equations. The continuity equation is satisfied automatically since the flow is solely azimuthal around the axis of symmetry of the equilibrium. The fluid statement of force balance is given by

(5.2)\begin{equation} \boldsymbol{J}\times\frac{\boldsymbol{B}}{c} - \boldsymbol{\nabla}(nT) - qn\boldsymbol{\nabla}\phi =0, \end{equation}

where $\boldsymbol {J} = qnr\omega \hat {\boldsymbol {\varTheta }}$ is the electric current density. In the context of the local thermal equilibrium, the potential refers to the electrostatic potential of the trap and the plasma $\phi = \phi _{{\rm trap}} + \phi _p$. We only consider cases where the centrifugal term is small. For a locally Maxwellian velocity distribution the pressure tensor is a diagonal and can be represented by the scalar $p(\alpha, \chi ) = n(\alpha, \chi )T(\alpha )$. Hence, (5.2) states that the Lorentz force balances the Coulomb and the pressure force.

Using the relations $\boldsymbol{J}\times\boldsymbol{B} = qnr\omega B\hat{\boldsymbol{\Theta}} \times \hat{\boldsymbol{\chi}}$ and $\hat {\boldsymbol {\varTheta }} \times \hat {\boldsymbol {\chi }} = \hat {\boldsymbol {\alpha }}$ and taking the scalar product of (5.2) with $\hat {\alpha }$ results in a force-balance relation transverse to the magnetic flux contour. This relation yields the rotation frequency

(5.3)\begin{equation} \omega(\alpha, \chi) = {c}\frac{\partial \phi(\alpha, \chi)}{\partial \alpha} + \frac{{c}}{qn(\alpha, \chi)} \frac{\partial (n(\alpha, \chi)T(\alpha))}{\partial \alpha}, \end{equation}

where $rB \partial / \partial \alpha = \hat {\boldsymbol {\alpha }} \boldsymbol {\cdot } \boldsymbol {\nabla }$ is a derivative transverse to the magnetic flux contour in the direction of increasing flux. Result (5.3) was obtained from force balance, but can also be obtained from guiding-centre drift theory. Indeed, (5.2) can be obtained by expressing the current $J$ as the sum of currents from the $E\times B$ drift, the curvature drift, the grad $B$ drift and the diamagnetic drift (Bellan 2006). This is a consequence of the equivalence between double-adiabatic MHD theory and drift theory. Here, the temperature is isotropic so double-adiabatic MHD reduces to ordinary MHD. The first term in (5.3) is the $E\times B$ drift and the second term is the sum of the curvature drift, the grad $B$ drift and the diamagnetic drift. This latter result may be surprising to readers with experience in local thermal equilibrium states for plasmas in Penning–Malmberg traps, where the magnetic field is homogeneous and there is no curvature drift and grad $B$ drift. For such a trap, the second term in (5.3) would be just the diamagnetic drift. In the curved dipole field, the same term represents the sum of the curvature drift, grad $B$ drift and diamagnetic drift. There are subtle cancellations in the sum.

Next, we consider force balance along a field line. Taking the scalar product of (5.3) with $\hat {\boldsymbol {\chi }}$ yields

(5.4)\begin{equation} \frac{T(\alpha)}{n(\alpha, \chi)}\frac{\partial n(\alpha, \chi)}{\partial \chi} +q\frac{\partial \phi(\alpha, \chi)}{\partial \chi} = 0, \end{equation}

which implies the Boltzmann factor

(5.5)\begin{equation} n(\alpha, \chi) = h(\alpha)\exp\left[-\frac{q\phi(\alpha, \chi)}{T(\alpha)}\right]. \end{equation}

Depending on the provided input parameters, $h(\alpha )$ can be expressed differently. For a known number of particles on each flux surface

(5.6)\begin{equation} N(\alpha) = \int {\rm d} \chi \frac{2{\rm \pi} n(\alpha, \chi)}{B^2(\alpha, \chi)}, \end{equation}

we can write $h(\alpha )$ as

(5.7)\begin{equation} h(\alpha) = \frac{N(\alpha)}{\displaystyle \int {\rm d} \chi \frac{2{\rm \pi}}{B^2(\alpha, \chi)}\exp\left[-\frac{q\phi(\alpha, \chi)}{T(\alpha)}\right]}. \end{equation}

Thus, force balance plus the two functions $N(\alpha )$ and $T(\alpha )$ determine the density distribution of a local thermal equilibrium. Alternatively, for a given density profile $\hat {n}(\alpha, \chi _0)$ along some contour of constant $\chi = \chi _0 = {\rm const.}$, we can write

(5.8)\begin{equation} h(\alpha) = \hat{n}(\alpha, \chi_0) \exp\left[\frac{q\phi(\alpha, \chi_0)}{T(\alpha)}\right]. \end{equation}

Plugging this expression for $h(\alpha )$ into (5.5) yields the desired density profile $n(\alpha, \chi ) = \hat {n}(\alpha, \chi _0)$ for $\chi = \chi _0$. Finally, some physical insight can be obtained when expressing $h(\alpha )$ in the form

(5.9)\begin{equation} h(\alpha) = n_0 \exp\left[-\frac{q\tilde{\phi}(\alpha)}{T(\alpha)}\right], \end{equation}

with the constant $n_0$. In this form it is easy to see that the local thermal equilibrium is also a global thermal equilibrium when $T(\alpha ) = T$ and $\omega (\alpha ) = \omega$ are independent of $\alpha$ and $q\tilde {\phi } = -q\omega \alpha /{c}$. Note here that, within guiding-centre drift dynamics, the canonical momentum reduces to $p_{\varTheta } = mr^2\dot {\varTheta } + q\alpha /{c} \approx q\alpha /{c}$.

Using the density expressed in terms of $\tilde {\phi }(\alpha )$ to evaluate (5.3) for the rotation frequency yields the alternate expression

(5.10)\begin{equation} \omega ={-}{c}\frac{\partial \tilde{\phi}(\alpha)}{\partial \alpha} + \frac{{c}}{q} \frac{\partial T(\alpha)}{\partial \alpha} \left[ 1 - \ln \left(\frac{n(\alpha, \chi)}{n_0} \right) \right]. \end{equation}

If the temperature is independent of $\alpha$, the rotation frequency reduces to $\omega (\alpha ) = -{c}\partial \tilde {\phi }/\partial \alpha$. Likewise, the variation of the rotation frequency along a magnetic field line is given by the expression

(5.11)\begin{equation} \frac{\partial \omega}{\partial \chi} ={-}\frac{{c}}{qn}\frac{\partial T}{\partial \alpha}\frac{\partial n}{\partial \chi}. \end{equation}

To explore the low-temperature or small-Debye-length limit, we again write the plasma potential as an expansion $\phi _p \approx \phi _p^{(0)} + \phi _p^{(1)}$, where $\phi _p^{(0)}$ is zero order in temperature and $\phi _p^{(1)}$ is first order in temperature. Setting the temperature to zero in (5.4) results in $n\partial \phi _p^{(0)} / \partial \chi = 0$, a constant electric potential in the presence of the cold plasma. Comparing the rotation frequency according to (5.3) and (5.10) yields $\omega /{c} = \partial \phi _p^{(0)} / \partial \alpha = - \partial \tilde {\phi } / \partial \alpha$ and hence $\phi _p^{(0)}(\alpha ) + \tilde {\phi }(\alpha ) = {\rm const.}$ in the cold limit. Using these results in Poisson's equation gives us the zero-temperature density

(5.12)\begin{equation} n_{{\rm cold}} ={-}\frac{1}{4{\rm \pi} q}\nabla^2 \phi_p^{(0)} ={-}\frac{B^2}{4{\rm \pi} {c} q}(\alpha, \chi)\frac{\partial}{\partial \alpha}(r^2(\alpha, \chi) \omega ). \end{equation}

We can identify $\partial r^2 / \partial \alpha = 2B_z/B^2$ to retrieve a generalization of (2.11) for the case of a local thermal equilibrium

(5.13)\begin{equation} n_{{\rm cold}} ={-}\frac{1}{4{\rm \pi} q}\left(2m\omega\omega_{cz} + \frac{r^2B^2}{{c}}\frac{\partial \omega}{\partial \alpha}\right). \end{equation}

This expression is identical to the one found by Pétri et al. (2002) for the magnetosphere of a neutron star.

From force balance (5.4) we obtain the variation of the plasma potential to first order in temperature along a magnetic field line

(5.14)\begin{equation} \frac{\partial}{\partial \chi} \left[ T(\alpha) ln\left(\frac{n_{{\rm cold}}(\alpha, \chi)}{n_0}\right) + q\phi_p^{(1)}\right] = 0. \end{equation}

The following numerically determined local thermal equilibrium states will corroborate these theoretical findings.

6. Numerical solutions for local thermal equilibrium states

The basic method to calculate solutions for local thermal equilibrium states is the same as for global thermal equilibrium states. Poisson's equation is still solved in cylindrical coordinates, but the density distribution given by (5.5) has to be treated in ($\alpha, \chi$)-coordinates. We use the density profile of the global thermal equilibrium result from figure 5 along the equatorial plane with the constant temperature $\lambda _{D0}/R = 0.08$ as input parameter. This choice allows us to illustrate that a global thermal equilibrium is a special case of a local thermal equilibrium. However, it is an arbitrary choice. The function $h(\alpha )$ can be evaluated according to (5.8). Since the magnetic flux $\alpha$ can be calculated on each grid point in cylindrical coordinates according to (2.5), we can interpolate $h(\alpha )$ onto the rectangular grid without the need to introduce a grid in ($\alpha, \chi$)-coordinates.

For the alternative method that uses the number of particles per field line as input parameter such a grid in magnetic field line coordinates becomes necessary to evaluate the integral in (5.7). We then interpolate the potential onto the magnetic field lines, determine the density with the desired number of particles per field line and interpolate the density back onto the cylindrical grid to solve Poisson's equation.

If the trap configuration and the plasma temperature is the same as for the global thermal equilibrium from which we recovered our input parameter, the resulting local thermal equilibrium is identical to that global thermal equilibrium. This is the case in the upper left panel in figure 8. For the remaining results in figure 8 the charge on the coil is lowered and the density profile is scaled to keep the total number of particles unchanged. With less charge on the coil, the potential well vanishes and the plasma is in a local, but not a global, thermal equilibrium state. Note that only the trap configuration has changed in the transition from a global to the local thermal equilibrium states. The input parameters for the plasma itself were left unchanged.

Figure 8. The upper left panel corresponds to the hottest global thermal equilibrium in figure 5. In the other panels the charge of the coil was lowered while the temperature and total number of particles remains unchanged. The solid white line in the lower right panel indicates a magnetic field line along which the potential and rotation frequency is evaluated.

All the solutions in figure 8 obey the statement of force balance (5.2). This statement is validated along the magnetic flux contour indicated by the white solid line in figure 8. According to (5.4), force balance along the magnetic flux contour implies that the sum of the electrostatic potential and the pressure term is constant. In figure 9 the blue line for the electrostatic potential and the red line for the pressure term indeed add up to the constant magenta line. This is also true if the temperature varies across the flux contours.

Figure 9. Electric potential (blue), pressure term (red) and the sum of the two (magenta) along the field line indicated in figure 8. The magnetic potential $\chi$ is normalized by the magnetic field in the centre of the coil and the coil radius.

Force balance across magnetic flux contours yields (5.10) for the rotation frequency. figure 10 shows the rotation frequency along the magnetic field line indicated by the white solid line in figure 8. In figure 10(a) the temperature is constant and in figure 10(b) we introduced a temperature gradient across the magnetic flux contours. The sum of the drifts result in a constant rotation frequency along a field line, only if the temperature is constant across field lines. The variation of the rotation frequency in figure 10(b) agrees with (5.10). Hence, the statement of force balance is satisfied.

Figure 10. The $E \times B$ rotation (blue), diamagnetic rotation (red) and the sum of the two (magenta) along the field line indicated in figure 8. Panel (a) shows the case with no temperature gradient and (b) with a temperature gradient. The black dashed line in (b) is the deviation from constant rotation along the magnetic flux surface according to (5.10). The frequency $\omega _0$ refers to the rotation to the corresponding global thermal equilibrium.

Finally, we explore the zero-temperature limit of the local thermal equilibrium according to (5.13). Similar to the case of a cold global equilibrium, the challenge in finding a self-consistent solution for the cold local equilibrium is to determine the location of the edge of the plasma. The zero-temperature limit in figure 11 is calculated with respect to the mean potential along a field line $\langle \phi \rangle _{\alpha }$ of the coldest displayed finite-temperature solution. The density according to (5.12) is calculated in the regions where $|\phi (\alpha ) - \langle \phi \rangle _{\alpha } |/\phi (\alpha ) < \delta \phi$. Regions where this condition is not satisfied are governed by Laplace's equation. Introducing the cold solution will change the potential. Consequently, the edge of the cold solution needs to be adjusted in an iterative manner. We start with the potential of the finite temperature solution and adjust it to the cold density with a weighting factor $\gamma$ as it was described for the global thermal equilibrium. The allowed deviation from the input potential $\delta \phi$ is lowered as the solution converges. This is the least robust of the methods described in this paper and requires an initial guess that is close to the cold solution. The displayed cold solution has an accuracy of $\delta \phi = 0.1\,\%$ and converged at $N_{{\rm error}}/N_{{\rm tot}} = 10^{-4}$.

Figure 11. Density distribution for local thermal equilibria with different temperatures. The hottest solution in the upper left panel is identical to the one with zero charge on the coil in figure 8. The lower right panel shows the corresponding zero temperature limit. The white line in the upper left panel indicates a magnetic field line.

The density distributions in figure 11 are local thermal equilibria with decreasing temperatures (upper left to lower right). The temperature is expressed in terms of the Debye length $\lambda _{D0}$ for the mean density $n_0$ divided by the coil radius $R$. The density distribution in the lower right panel was obtained with the method described above for the zero-temperature limit. As the temperature decreases we expect the pressure term in (5.4) to decrease. As the pressure term vanishes in the cold limit the electrostatic potential along a magnetic flux contour is constant. The electrostatic potential for different temperatures along the magnetic flux contour indicated by the white solid line is displayed figure 12. The variation of the electrostatic potential decreases with decreasing temperature as expected. For the zero-temperature limit we obtain a deviation from the desired constant potential $\langle \phi \rangle _{\alpha }$ within the error margin $\delta \phi = 0.1\,\%$ of our method.

Figure 12. Electrostatic potential along the magnetic field line indicated in figure 11 for decreasing temperatures (red to blue) including the zero-temperature limit. The magnetic potential $\chi$ is normalized by the magnetic field in the centre of the coil and the coil radius.

7. Conclusion

We have discussed both local thermal equilibrium states on magnetic flux surfaces and global thermal equilibrium states for single species plasmas in a azimuthally symmetric dipole trap. Confinement of global thermal equilibrium states requires that the effective trap potential in the rotating frame forms a potential well. For the simplest dipole trap, a circular current carrying coil in unbounded space, the coil is at the bottom of the potential well. A charged coil can form a potential barrier but also destroys the potential well. This issue was resolved when we introduced a tailored conducting vacuum chamber around a levitated coil to provide a potential well that does not contain the coil. We considered global thermal equilibrium states with a finite temperature as well as the zero-temperature limit in such a configuration. To zeroth order in temperature, the plasma potential balances the vacuum potential. An analytic expression for the corresponding density can be found. However, this expression only provides a self-consistent solution within a contour of constant effective potential, resulting in an abrupt edge. Finding this potential contour requires a numerical approach. To first order in temperature, the pressure is balanced. For finite temperatures the plasma density falls off on the scale of the Debye length. This density distribution is governed by the Poisson–Boltzmann equation, which again requires a numerical solver. Results for plasmas in a global thermal equilibrium state with finite temperature as well as in the zero-temperature limit were obtained. The characteristic features of force balance and rigid rotation were validated. The global thermal equilibrium was shown to be a special case of a local thermal equilibrium. Local thermal equilibrium states are characterized by an isotropic, Maxwellian velocity distribution with an azimuthal drift velocity and force balance on magnetic flux surfaces. Force balance as well as the consequent drift velocity are retrieved in the numerical results. The rotation frequency along a field lined turned out to be constant only if the gradient of the temperature across field lines or the gradient of the density along field lines vanishes. The treatment of the cold local equilibrium revealed a generalization of the known analytic expression for a cold global equilibrium. The transition to the zero-temperature limit was demonstrated for which the electrostatic potential of the numerical result is constant along a field line.