2. Theory of global thermal equilibrium in a dipole field

For a real dipole trap experiment, such as the RT1 experiment (Ogawa et al. Reference Ogawa, Yoshida, Morikawa, Saitoh, Watanabe, Yano, Mizumaki and Tosaka2009) or the planned APEX (A Positron-Electron eXperiment) (Stoneking et al. Reference Stoneking, Pedersen, Helander, Chen, Hergenhahn, Stenson, Fiksel, von der Linden, Saitoh and Surko2020), 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 Reference Dubin and O'Neil1999)

(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 (Reference Goldreich and Julian1969) and Turner (Reference Turner1991). 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 Reference Prasad and O'Neil1979). 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 Reference Jackson1999). 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. Reference Pétri, Heyvaerts and Bonnazzola2002). 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. Reference Spencer, Rasband and Vanfleet1993). 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. Reference Hyatt, Driscoll and Malmberg1987) and lasts for the transport time scale (Driscoll et al. Reference Driscoll, Malmberg and Fine1988). 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 (Reference Pedersen and Boozer2002).

In Sato et al. (Reference Sato, Kasaoka and Yoshida2015) 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. (Reference Sato, Kasaoka and Yoshida2015) can be anisotropic. The anisotropic temperature relaxation rate for an electron plasma in a Penning–Malmberg trap was measured by Hyatt et al. (Reference Hyatt, Driscoll and Malmberg1987) and is in agreement with collisional Fokker–Planck theory. The states obtained by Sato et al. (Reference Sato, Kasaoka and Yoshida2015) 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 (Reference Taylor1964). 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 Reference Bellan2006). 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. (Reference Pétri, Heyvaerts and Bonnazzola2002) 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.