We are developing a positron–electron plasma trap based on a dipole magnetic field generated by a levitated superconducting magnet to investigate the physics of magnetized plasmas with mass symmetry as well as antimatter components. Such laboratory magnetosphere is deemed essential for the understanding of pair plasmas in astrophysical environments, such as magnetars and blackholes, and represents a novel technology with potential applications in antimatter confinement and the development of coherent gamma-ray lasers. The design of the device requires a preemptive analysis of the achievable self-organized steady states. In this study, we construct a theoretical model describing maximum entropy states of a collisionless positron–electron plasma confined by a dipole magnetic field and demonstrate efficient confinement of both species under a wide range of physical parameters by analyzing the effect of the three adiabatic invariants on the phase space distribution function. The theory is verified by numerical evaluation of spatial density, electrostatic potential, and toroidal rotation velocity for each species in correspondence with the maximum entropy state.

## I. INTRODUCTION

Pair plasmas consist of two species of charged particles having the same mass, but opposite electric charge. Such mass symmetry is in stark contrast with usual plasmas, where ions exhibit a significantly larger mass than the electron component, and it is the reason why pair plasmas are expected to possess peculiar turbulence, stability, and fluctuation properties.^{1} Positron–electron plasmas are one example of naturally occurring pair plasmas and can be found in astrophysical jets and magnetospheres of quasars,^{2–4} pulsars,^{5} and magnetars,^{6–8} where positrons and electrons are formed by pair production. More exotic pair plasmas, such as proton–antiproton plasmas, could occur as well as in suitable astrophysical environments.

Positron–electron plasma experiments based on magnetic mirrors have been considered in the past^{9} and efforts are under way for their realization,^{10} while pair plasmas involving oppositely charged fullerene ions have been produced in the laboratory.^{11,12} Furthermore, positron–electron pair production is expected to occur during the disruptive phase of tokamaks.^{13} At present, positron–electron plasma confinement schemes based on dipole magnetic fields represent an active area of research due to recent technological advancement in trap design,^{14–17} positron accumulation,^{18,19} and positron injection.^{20,21} Furthermore, the inhomogeneity of a dipole magnetic field makes these magnetic configurations suitable to confine both neutral and nonnetutral plasmas, in contrast with magnetic traps relying on a straight homogeneous magnetic field, which are appropriate for nonneutral plasmas.^{22–24}

If realized, a positron–electron plasma based on a dipole magnetic field has several potential applications ranging from the experimental study of astrophysical systems, such as astrophysical jets, quasars, and pulsars, to matter and antimatter confinement, potentially providing confinement to a large amount of charged particles (mainly depending on the properties of the matter/antimatter source) for a long time (300 s and beyond^{25}) without requiring external electric fields and development of coherent gamma-ray sources (gamma-ray lasers). In the laser application envisioned here, the dipole trap would confine a dense low-temperature positron–electron plasma. The coherent gamma-ray source would then be obtained by the process of stimulated annihilation of positron–electron pairs through the action of impinging photons of suitable energy (annihilation laser^{26–28}) rather than by stimulated emission.

The target parameters for a laboratory positron–electron plasma are expected to ensure a high degree of stability not seen in standard ion–electron plasmas.^{29,30} In such a confinement regime, the gyroradius $ r c = m v \u22a5 / | q | B$ is smaller than the Debye length $ \lambda D = ( \u03f5 0 k B T / 2 n e 2 ) 1 / 2$, which is smaller than the system size *R*, i.e., $ r c \u226a \lambda D \u226a R$. Here, *m* denotes the particle mass, $ v \u22a5$ is the velocity perpendicular to the magnetic field ** B**,

*q*is the electric charge,

*B*is the modulus of the magnetic field,

*ϵ*

_{0}is the vacuum permittivity,

*k*is the Boltzmann constant,

_{B}*T*is the temperature,

*n*is the spatial density, and $ e = | q |$. For a positron–electron plasma with $ B \u223c 1 \u2009 T , n \u223c 10 12 \u2009 \u2009 m \u2212 3$, $ k B T \u223c 10 \u2009 \u2009 eV$, and $ R \u223c 1 \u2009 \u2009 m$, one obtains $ r c \u223c 10 \u2212 5 \u2009 m$ and $ \lambda D \u223c 2 \xd7 10 \u2212 2 \u2009 \u2009 m$. Furthermore, the plasma parameter $ \Lambda = 4 \pi n \lambda D 3 = ( 1 / 4 2 \pi ) ( r d / r C ) 3 / 2 \u223c 10 8 \u226b 1$ is large, implying that the typical distance among particles $ r d = n \u2212 1 / 3 \u223c 10 \u2212 4 \u2009 m$ is much larger than the distance $ r C = e 2 / 4 \pi \u03f5 0 k B T \u223c 10 \u2212 10 \u2009 m$ at which the Coulomb energy becomes comparable with the average kinetic energy. In particular, the Coulomb energy at the typical distance

*r*is $ E C = e 2 / 4 \pi \u03f5 0 r d \u2248 2 \xd7 10 \u2212 24 J$. This value is small compared to the kinetic energy since $ E C / k B T \u2248 10 \u2212 6$. This makes the system weakly coupled, and diffusive (entropy maximizing) processes are dominated by collective electromagnetic fluctuations rather than localized Coulomb collisions.

_{d}At present, we are developing a positron–electron plasma trap based on a dipole magnetic field generated by a levitated superconducting coil operating in the plasma regime described above. This requires a preemptive analysis of the achievable plasma confinement. The aim of this study is, thus, to explore the nature of maximum entropy states that are self-organized by the plasma and assess the degree of confinement of both species. In this context, a maximum entropy state is defined as the distribution function of the largest entropy that is compatible with the macroscopic constraints affecting the system, such as conservation of total particle number, total energy, or total magnetic moment. Since a positron–electron system is formally analogous to any two-species plasma, such as an ion-electron plasma, thermal equilibria in a positron–electron plasma can be inferred from those of a two-species plasma. It is, therefore, important to stress that the core issue examined in the present paper is not the study of two-species thermal equilibria, which is an established matter, but the elucidation of the effect of the conservation of adiabatic invariants on self-organized maximum entropy states in two-species plasmas, a problem that is not discussed in the literature.

Once built, the positron–electron plasma trap will be used to study both waves and transport phenomena, which are expected to be in the $ 1 \u2009 \u2009 kHz$ and $ 1 \u2009 \u2009 Hz$ frequency range, respectively, based on previous experimental data from the RT-1 device.^{25} Hence, we aim at achieving confinement on time scales *τ _{c}* of the order of $ 1 \u2009 \u2009 s$. Since in the present setting the Coulomb scattering frequency is of the order $ \nu C \u223c 1 \u2009 \u2009 Hz$, the plasma is expected to be quite collisionless, while the main mechanism increasing the entropy of the system is given by collective electromagnetic fluctuations. For example, a fluctuating electric field

**can induce a random walk (diffusion process) across the magnetic field by exciting $ E \xd7 B$ drift motion. Hence, we shall develop the theory by assuming that Coulomb collisions can be neglected.**

*E*It should be emphasized that accumulating 10^{11}–10^{12} positrons in a volume $ \Omega \u223c 1 \u2009 m 3$ requires the development of positron accumulation and injection technologies. In the present trap design, positrons will be produced by a $ 10 \u2009 \u2009 eV$ pulsed positron source located at The National Institute of Advanced Industrial Science^{18} and progressively accumulated into a buffer-gas trap until a total positron number $ N p \u223c 10 11$ is reached. The system above is also designed to minimize the positron energy spread. The positrons will then be released into the dipole trap over a time interval (injection timescale) $ \tau i \u223c 10 \u2009 \u2009 \mu s \u226a \tau c \u223c 1 \u2009 \u2009 s$, reaching a particle density in the target range 10^{11}– $ 10 12 \u2009 m \u2212 3$. Notice here that electrons “wait” inside the dipole field for the positron injection (electrons can be confined inside the dipole trap for time scales of the order $ \u223c 300 \u2009 s$). Efficient injection schemes are also being developed to minimize losses at the injection phase. Once injected, the main source of particle loss is not expected to be the interaction with the confining vessel or the current loop box, but charge exchange with neutral hydrogen impurities, resulting in the formation of positronium. The loss timescale for this process can be roughly estimated to be $ \tau l \u223c 10 \u2009 s$. This value is sufficiently longer than the confinement timescale *τ _{c}*, which is itself longer than the relaxation timescale

*τ*(the time needed by the system to approach the equilibrium state). In particular, experimental evidence from the RT-1 device

_{r}^{31}suggests

*τ*to be a small fraction of

_{r}*τ*(although precise estimates for a positron–electron plasma would need experimental verification). The injection, relaxation, confinement, and loss time scales are, therefore, related according to $ \tau i < \tau r < \tau c < \tau l$.

_{c}In addition to the technological hurdles listed above, other factors, such as particle loss associated with turbulent transport and plasma instability, may affect the quality of confinement. These aspects, which have been investigated in dipole geometry within the frameworks of MHD and gyrokinetics (see, e.g., Refs. 32–34) will not be addressed in the present paper since its scope is limited to establishing the existence of stable maximum entropy equilibria toward which the positron–electron plasma is expected to converge under ideal conditions.

In an inhomogeneous magnetic field such as a dipole magnetic field, the properties of the equilibrium states are strongly dependent on the presence of adiabatic invariants.^{35} When the timescale of collective electromagnetic fluctuations within the plasma is longer than the period of cyclotron gyration, bounce motion, or toroidal drift, the corresponding adiabatic invariant (magnetic moment *μ*, bounce action $ J \u2225$, or magnetic flux $\Psi $) is conserved by the dynamics of a single charged particle. Indeed, the essential feature of adiabatic invariants is that they remain approximately constant as long as perturbations acting on a dynamical system are slow compared with the period of adiabatic motion. This results in a set of constraints on the equilibrium distribution function of the system, which departs from a standard Maxwell–Boltzmann distribution.^{36–38} To understand how this occurs, it is useful to consider the limiting case in which the value of one adiabatic invariant, say *μ*, is preserved exactly. Then, on each phase space submanifold corresponding to a level set of *μ* (a *μ*-leaf), the effective phase space measure is reduced to $ Bdxdydzd v \u2225$, where $ ( x , y , z )$ are Cartesian coordinates and $ v \u2225$ denotes the velocity along the magnetic field. If $ f \mu ( x , y , z , v \u2225 )$ is the distribution function on a *μ*-leaf, entropy maximizing processes such as cross field diffusion lead to a flattening of $ f \mu $, with a corresponding relaxed spatial density $ n = B \u222b f \mu d v \u2225 d \mu \u223c B$, which is inhomogeneous for an inhomogeneous *B* (for further details, see Refs. 38–40). Therefore, a statistical description of maximum entropy states must take into account the nontrivial role played by adiabatic invariants in shaping the distribution function of each particle species.

The present paper is organized as follows. In Sec. II, we derive the maximum entropy distribution functions of positrons and electrons by taking into account the conservation of the first adiabatic invariant and obtain the corresponding form of the Poisson equation for the electrostatic potential. In Secs. III and IV, we study the effect of the second and third adiabatic invariants on positron–electron plasma maximum entropy states. In Sec. V, we report spatial densities and electrostatic potential obtained by numerical solution of the Poisson equation for the electrostatic potential and demonstrate efficient confinement of both species. In Sec. VI, we study the toroidal rotation velocity profile of the positron–electron plasma and find that it departs from the rigid rotation occurring in magnetic traps relying on a straight homogeneous magnetic field. Concluding remarks are given in Sec. VII.

We remark that the theory developed in the present paper holds provided that the macroscopic constraints defining each statistical ensemble (total particle number, total energy, total magnetic moment, and so on) are preserved throughout the relaxation of the system. In particular, the period of the slowest adiabatic motion must be considerably smaller than the typical timescale of the fastest turbulent fluctuations in the system. However, in real experiments, particles and energy will be lost via different mechanisms, such as by interaction with the vessel boundary, and the degree of conservation of adiabatic invariants will not be perfect. Therefore, the theory will be quantitatively accurate only if such losses are not too large. Nevertheless, we expect the theory to remain qualitatively consistent even in the presence of large losses since the surviving (trapped component) of the plasma should eventually converge toward the derived maximum entropy states. It is also worth observing that if one considers time scales longer than the typical interval required for a Coulomb collision to occur, the conservation of the adiabatic invariants of each charged particle is progressively broken due to cumulative changes in the particle energies. Nevertheless, we conjecture the theory to remain approximately valid even in such regime since it only relies on the hypothesis that the total value of each adiabatic invariant is preserved, and the total change due to multiple Coulomb scatterings should be, on average, negligible.

Finally, notice that the present theory does not make any predictions on the nature of electromagnetic fluctuations. It only provides information on the maximum entropy state that results from the action of these fluctuations, which act to increase the entropy of the system. Hence, although the model cannot describe transient phenomena such as waves, transport, or diffusion, it is sufficient to describe the “final” state of the system, i.e., the equilibrium toward which the plasma tends to settle over a sufficiently long time interval.

## II. SPATIAL DENSITY AND ELECTROSTATIC POTENTIAL IN MAGNETIZED POSITRON–ELECTRON PLASMA

*p*and

*e*will be used to specify positrons and electrons, respectively. In a static equilibrium, the electrostatic potential $\Phi $ is determined by the Poisson equation,

*ρ*is the electric charge density and

*ϵ*

_{0}is the vacuum permittivity. In a positron–electron plasma, the charge density can be expressed as

*e*is the positron charge, –

*e*is the electron charge,

*n*is the positron number density, and

_{p}*n*is the electron number density. Denoting with

_{e}*f*the probability distribution function of positrons,

_{p}*f*the probability distribution function of electrons,

_{e}*N*the total number of positrons, and

_{p}*N*the total number of electrons, the number densities

_{e}*n*and

_{p}*n*can be evaluated as

_{e}*μ*is a constant of motion of (isolated) charged particle dynamics. Here, the magnetic moment

*μ*

^{41–43}is defined as

*v*,

_{c}*B*is the modulus of the magnetic field

**, and**

*B**m*is the particle mass (which is the same for both positrons and electrons).

*T*of fluctuations in particle energy caused by electromagnetic turbulence must be longer than the timescale of cyclotron motion

_{f}*T*, i.e.,

_{c}*ω*denotes the cyclotron frequency. In the following, we shall, therefore, assume that the condition (5) holds. We also recall that the conservation of the quantity (4), usually referred to as the lowest order magnetic moment, rests on the additional requirement that the wavelength

_{c}*λ*of any fluctuation is longer than the Larmor radius

*ρ*, i.e., $ \rho c \u226a \lambda $. In this setting, the guiding center Hamiltonian functions for positron and electrons have expressions

_{c}*L*is the characteristic scale length of the magnetic field, and

*k*is a characteristic plasma temperature, represents a small perturbation of the particle energy (for additional details on this point, see discussion at the end of Sec. IV or Ref. 42). Then, the Shannon entropies

_{B}T*S*and

_{p}*S*, the total probabilities

_{e}*P*and

_{p}*P*, the total energies

_{e}*E*and

_{p}*E*, and the total magnetic moments

_{e}*M*and

_{p}*M*of the system can be written as

_{e}*P*,

_{p}*P*,

_{e}*E*,

_{p}*E*,

_{e}*M*, and

_{p}*M*(for further details on this approach, see Refs. 37 and 38). Introducing Lagrange multipliers

_{e}*α*,

_{p}*α*,

_{e}*β*,

_{p}*β*,

_{e}*γ*, and

_{p}*γ*, we, therefore, define the target functional

_{e}*f*and

_{p}*f*are

_{e}*β*and

_{p}*β*represent the characteristic inverse temperatures of the two particle species, the Lagrange multipliers

_{e}*γ*and

_{p}*γ*can be interpreted as chemical potentials describing the macroscopic energy changes $ \beta p d E p = \gamma p d M p$ and $ \beta e d E e = \gamma e d M e$ occurring when magnetic moments

_{e}*dM*,

_{p}*dM*are added to the system.

_{e}Next, it is useful to explain why Coulomb collisions have been neglected in the derivation of the collisionless equilibria (9). First, recall that, in the present setting, the (collective) electric potential $\Phi $ changes over a timescale *T _{f}* that is much longer than the cyclotron timescale

*T*. This ensures that the individual magnetic moments of the particles remain constant, unless local Coulomb collisions occur. For the plasma regime under consideration, the frequency

_{c}*ν*of such collisions can be estimated

_{C}^{44}as $ \nu C \u223c 5 \xd7 10 \u2212 11 \u2009 nT \u2212 3 / 2 \u223c 1 \u2009 \u2009 Hz$, where $ n = 10 12 \u2009 \u2009 m \u2212 3$ is the particle density and $ T = 10 \u2009 \u2009 eV$ is the temperature. Hence, over the confinement time scales $ \tau c \u223c 1 / \nu C \u223c 1 \u2009 \u2009 s$ considered in the model they are negligible, their effect being felt only over longer time intervals where particles undergo mostly small deflections with small changes in their kinetic energies. Indeed, recall that the distance $ r C \u223c 10 \u2212 10 \u2009 m$ at which Coulomb interactions are dominant is smaller than the typical particle distance $ r d \u223c 10 \u2212 4 \u2009 \u2009 m$.

We also remark that if a net transfer of total energy or total magnetic moment occurs between positrons and electrons, only the sums $ E p + E e$ and $ M p + M e$ are preserved. Such scenario can be taken into account by setting $ \beta = \beta e = \beta p$ and $ \gamma = \gamma p = \gamma e$.

*f*and

_{p}*f*can now be used to determine the spatial number densities

_{e}*n*and

_{p}*n*according to Eq. (3). To evaluate the integral in Eq. (3), the phase space measure $ d \Pi $ must be expressed in a more convenient set of magnetic coordinates. To this end, we restrict our attention to magnetic fields of the type

_{e}**, $ v c$ is the cyclotron velocity such that $ m v c 2 = 2 \mu B$ [recall (4)], and $ v d = v d ( x )$ is the particle drift velocity across**

*B***. Notice that Eq. (12) represents the velocity of a charged particle, and not the guiding center velocity. Furthermore, the drift velocity $ v d$ comprising $ E \xd7 B$, gradient, and curvature drifts is treated as a spatial function, which is expected to be a valid approximation in a time-independent setting since in a vacuum field gradient and curvature drifts can be expressed as $ 2 k B T B \xd7 \kappa / q B 2$, where $\kappa $ is the field curvature,**

*B**k*is the temperature, and $ q = \xb1 e$ is the relevant electric charge. Let

_{B}T*r*denote the cylindrical radius, $ \u03d1 c$ the phase of the cyclotron gyration, and $ A = \Psi \u2207 \phi $ the vector potential associated with the magnetic field (10). Then, the particle momentum $ p = m v + q A$ can be decomposed on the orthonormal set of basis vectors

*n*and

_{p}*n*, therefore, have expressions

_{e}*M*and

_{p}*M*is expected to break down. This scenario corresponds to the limit $ \gamma p = \gamma e = 0$, which gives

_{e}*T*and

_{p}*T*associated with the Lagrange multipliers $ \beta p = ( k B T p ) \u2212 1$ and $ \beta e = ( k B T e ) \u2212 1$) are not spatially uniform. Indeed, integrals of the type,

_{e}## III. THE EFFECT OF BOUNCE MOTION

*μ*and the second adiabatic invariant (bounce action) $ J \u2225$, which is defined by

*a*and

*b*denote the bouncing points along a field line with line element

*ds*, while

*E*, magnetic moment

*μ*, and spatial position $ ( \u2113 , \Psi , \phi )$. The conservation of the second adiabatic invariant (27) requires that timescale

*T*of energy fluctuations caused by electromagnetic turbulence is longer than the timescale of bounce motion

_{f}*T*, i.e.,

_{b}*ω*denotes the bounce frequency. The (half) period

_{b}*T*of the bounce oscillation can be written as

_{b}*t*and

_{a}*t*are the instants at which the particle position reaches the bounce points

_{b}*a*and

*b*, while $ \u27e8 \u2009 \u27e9 b$ denotes averaging over a bounce oscillation.

The actual magnetic field strength within the planned positron–electron trap will not be symmetric under vertical reflections $ z \u2192 \u2212 z$ due to a support coil placed at the top of the device to keep the main superconducting coil generating the dipole magnetic field levitated. Therefore, an accurate model of the trap would need to take into account such asymmetry. This problem will not be considered here to simplify the analysis, and a pure dipole magnetic field will be assumed. For completeness, we remark that the maximum field strength $ B supp$ of the support coil is expected to be less than 10% of the local dipole magnetic field *B*. In addition, this error is expected to be significantly smaller in the confinement region close to the superconducting coil. Since the particle density is roughly a linear function of the field strength [recall Eq. (18)], we, therefore, expect the support coil to introduce a relative error for the computed spatial density profiles of less than 10%.

*B*, one has $ a = \u2212 b$ as well as $ \u2202 v \u2225 \u2217 / \u2202 \phi = 0$. Next, suppose that the number of positrons equals the number of electrons,

*N*=

_{p}*N*. Then, we may assume the electric potential energy $ q \Phi $ to be small compared to $ E \u2212 \mu B$ for the part of the bounce orbit close to the equatorial plane that contributes the most to $ J \u2225$. In this case, Eqs. (27) and (30) can be simplified to

_{e}*q*according to Eq. (32), and therefore, they have the same expression for both positrons and electrons. We also remark that, when $ N p \u2260 N e$, Eq. (32) is expected to remain a good approximation provided that the kinetic energy (temperature)

*k*is sufficiently large compared to $ q \Phi $. Now observe that, using the identity (31) for the kinetic energy along the magnetic field, the positron energy and the electron energy can be approximated as

_{B}T*H*,

_{p}*H*with $ H \u0303 p , H \u0303 e$ evidently go to zero when averaged over a bounce oscillation, $ \u27e8 H p \u2212 H \u0303 p \u27e9 b = \u27e8 H e \u2212 H \u0303 e \u27e9 b = 0$. In the following, we shall be concerned with relaxation time scales

_{e}*τ*longer than the bounce period, $ \tau \u226b T b$, and use the approximate expressions (33) for the particles energies. As it will be shown later, this approach simplifies calculations involving the particles distribution functions.

*ζ*and

_{p}*ζ*, the resulting expressions for the distribution functions

_{e}*f*and

_{p}*f*are

_{e}*H*and

_{p}*H*have been replaced with the approximated values $ H \u0303 p$ and $ H \u0303 e$. Our next goal is to evaluate the spatial particle densities (3) so that Poisson's equation (1) can be applied to compute $\Phi $. To this end, the phase space measure $ d \Pi $ must be expressed in a new set of coordinates that is appropriate to carry out integrals in momentum space. First, consider the change of variables $ ( \u2113 , v \u2225 , \phi , \Psi , \theta c , \mu ) \u2192 ( \u2113 , E , \phi , \Psi , \theta c , \mu )$ with $ E = m v \u2225 2 / 2 + \mu B + q \Phi $. Recalling (17), we, thus, have

_{e}*n*and

_{p}*n*can now be computed as follows. For the positron component, we have

_{e}*ω*is a function of $ J \u2225$,

_{b}*μ*, and $\Psi $. However, the dependence on the bounce action $ J \u2225$ disappears in the limit $ b / L \u226a 1$ in which the bounce orbit size

*b*is shorter than the characteristic magnetic field line length

*L*. Indeed, in this limit, we may expand the magnetic field in powers of $\u2113$ around the equatorial point $ \u2113 = 0$ of the dipole to obtain a second-order equation for the bounce position,

*B*has a minimum at $ \u2113 = 0$ for each magnetic surface $\Psi $. Setting $ B 0 = B ( 0 , \Psi )$ and $ B 0 \u2033 = ( \u2202 2 B / \u2202 \u2113 2 ) ( 0 , \Psi )$, the positive bounce point

*b*can, therefore, be approximated as

*E*can be decomposed into a perpendicular component $ E \u22a5 0 = \mu B 0$ and a parallel component $ E \u2225 0 = \alpha 0 p E \u22a5 0 = \alpha 0 p \mu B 0$, where $ \alpha 0 p ( \Psi ) > 0$ is the local positron temperature anisotropy at the equator. We may, therefore, estimate

*μ*and $\Psi $ only. Recalling Eq. (41), we, thus, arrive at the following expression for the positron density:

*n*by flipping the sign of the electric charge, and the Poisson equation (1) for the electrostatic potential $\Phi $ becomes

_{e}*ζ*has been neglected, and, thus, the conservation of the total bounce action $ J p$ is not felt by the ensemble, the effect of bounce motion on the particle distribution does not disappear. This is because bounce dynamics is encapsulated in the term $ \omega b J \u2225 / 2$ of the Hamiltonian $ H \u0303 p$ appearing in the distribution function

_{p}*f*. Observing that in the same limit $ \zeta e / \beta e \omega b \u226a 1$ the electron density becomes

_{p}## IV. THE EFFECT OF THE THIRD ADIABATIC INVARIANT

*n*of Eq. (46) is modified by an exponential factor $ exp \u2009 { \u2212 \eta p \Psi }$. A similar modification applies to

_{p}*n*so that the Poisson equation (47) for the electrostatic potential $\Phi $ now reads

_{e}*T*of energy fluctuations must be longer than the characteristic timescale

_{f}*T*of drift dynamics $ v d$ across the magnetic field

_{d}**in the toroidal direction $\phi $, i.e.,**

*B**ω*is the drift frequency. Since

_{d}*ω*is mainly determined by the $ E \xd7 B$, curvature, and gradient drift velocities, this frequency is usually smaller than the cyclotron and bounce frequencies, and (53) poses a rather stringent condition on the allowed turbulent spectrum of electromagnetic fluctuations.

_{d}Finally, we remark that the contributions to the particle energy coming from the guiding center drifts are neglected in the guiding center energies $ H \u0303 p , H \u0303 e$ under the assumption that the corresponding kinetic energy is smaller than the other terms. In particular, this is true if the kinetic energy associated with the $ E \xd7 B$ drift $ v E = E \xd7 B / B 2$ satisfies the ordering $ m | v E | 2 / 2 \u226a H \u0303 p , H \u0303 e$. Such configuration can be achieved, for example, when the electric potential $\Phi $ is a first-order contribution in the ratio $ \u03f5 \u223c \rho c / L$, where *ρ _{c}* is the Larmor radius and

*L*is the characteristic scale length of the magnetic field. For a plasma with inverse temperature

*β*, this implies $ \beta e | \Phi | \u223c \u03f5$ as well as $ \beta m 2 | v E | 2 \u223c \u03f5 2$, and the energies $ H \u0303 p , H \u0303 e$ are accurate at first order in

*ϵ*. For additional details on the expression of the guiding center Hamiltonian, see Ref. 42.

## V. NUMERICAL COMPUTATION OF SPATIAL DENSITIES AND ELECTROSTATIC POTENTIAL

The aim of this section is to numerically study the behavior of the positron–electron plasma as described by the model developed in Secs. II–IV for different values of the physical parameters involved.

For the spatial domain Ω occupied by the plasma, we consider a region $ r \u2208 [ r 0 , r 0 + R ] , \u2009 z \u2208 [ \u2212 R / 2 , R / 2 ] , \u2009 \phi \u2208 [ 0 , 2 \pi )$ mimicking an axially symmetric trap with radial size *R*, height *R*, and a physical axis (center stack) whose external boundary is located at $ r = r 0$. In particular, we choose $ R = 1 \u2009 m$ and $ r 0 = 0.1 \u2009 m$, which are values of the same order of the length parameters of the planned positron–electron trap. The dipole magnetic field used to confine the plasma is generated by a current loop of infinitesimal section and radius $ r l = 0.25 \u2009 \u2009 m$ enclosed in an axially symmetric toroidal box (coil) with squared cross section whose left side is located at $ r box = 0.2 \u2009 m$ (see Ref. 45 for the expression of the magnetic field). The typical magnetic field generated by the current loop is around $ B \u223c 0.1 T$ and reaches $ B \u223c 2 T$ in close proximity of the bounding box. It is also convenient to introduce the characteristic magnetic field $ B \u2217 = \mu 0 I / 2 r l = 1.25 \u2009 T$, where *μ*_{0} is the vacuum permeability and *I* is the electric current flowing within the current loop. The setting described above is shown in Fig. 1.

We will now examine the plasma equilibria corresponding to the statistical ensembles constructed in Secs. II–IV separately. To this end, we numerically solve the Poisson equation (1) for the electrostatic potential $\Phi $ by using the source term (2) obtained from the theoretical model. For example, the conservation of the first adiabatic invariant considered in Sec. II leads to the Poisson equation (20). Then, the distribution functions and the spatial densities can be evaluated according to Eqs. (9) and (18) by using the resulting value of $\Phi $. In order to solve the Poisson equation (1), we use the NDSolveValue solver of Wolfram Mathematica 12.2. The specific parameter values and boundary conditions used in each simulation are described in the corresponding sections.

### A. *μ* equilibrium

*μ*

Consider the positron–electron plasma equilibrium (9) arising from the conservation of total magnetic moments *M _{p}*,

*M*. We shall refer to such equilibrium as a

_{e}*μ*equilibrium. The condition for a

*μ*equilibrium to hold is that the timescale of the electromagnetic fluctuations driving the system toward the relaxed state is longer than the timescale of cyclotron dynamics as described by Eq. (5). In practice, this means that the term $ q \Phi $ occurring within the Hamiltonians

*H*,

_{p}*H*evolves over long time scales compared to the cyclotron timescale. This ensures that the constancy of the first adiabatic invariant

_{e}*μ*is not broken.

*ω*is

_{c}^{25}the spectrum of electromagnetic fluctuations is well below the $ 10 5 \u2009 Hz$ range. We, therefore, expect the condition (55) to be easily fulfilled by the system under examination. We shall, therefore, assume that (55) holds and apply the equilibrium model developed in Sec. II.

*A*,

_{p}*A*,

_{e}*N*,

_{p}*N*,

_{e}*β*,

_{p}*β*,

_{e}*γ*, and

_{p}*γ*appearing on the right-hand side of equation (20) are determined as follows. First, the target temperature

_{e}*T*for the positron–electron plasma is $ T \u2248 T p \u2248 T e \u2248 10 \u2009 \u2009 eV$, where

*T*and

_{p}*T*denote the temperatures of positrons and electrons, respectively. This fixes the inverse temperatures $ \beta p \u2212 1 = k B T p$ and $ \beta e \u2212 1 = k B T e$, where

_{e}*k*is the Boltzmann constant. The target spatial density

_{B}*n*is $ n \u2248 n p \u2248 n e \u2248 10 11 \u2009 m \u2212 3$ up to $ 10 12 \u2009 \u2009 m \u2212 3$. Since the volume of the positron–electron plasma trap is of the order $ \Omega \u2248 1 \u2009 \u2009 m 3$, we consider a total particle number of the order $ N \u2248 N p \u2248 N e \u2248 10 11$. Next, an estimate of the chemical potentials

*γ*and

_{p}*γ*can be obtained by observing that the changes in the total energies

_{e}*E*,

_{p}*E*caused by the addition of magnetic moments

_{e}*dM*,

_{p}*dM*to the system are given by $ \beta p d E p = \gamma p d M p$ and $ \beta e d E e = \gamma e d M e$. The order of magnitude of $ \gamma p , \gamma e$ is, therefore, expected to be $ \gamma p \u2248 \beta p H p / \mu \u2248 \beta p B$ and $ \gamma e \u2248 \beta e H e / \mu \u2248 \beta e B$. Finally, the parameters

_{e}*A*,

_{p}*A*, which represent measures of the volume Π occupied by the plasma in the phase space, can be estimated by observing that for a Maxwell–Boltzmann distribution

_{e}*R*is the radial size of the device and the value $ n \u2248 10 11 \u2009 \u2009 m \u2212 3$ has been used. In practice, $ \Phi c$ is expected to be much smaller than (58) since the charge separation $ e ( n p \u2212 n e )$ will be sensibly smaller than the characteristic charge density $ e n \u2248 10 \u2212 8 \u2009 \u2009 C \u2009 m \u2212 3$ of each individual species. In summary, the boundary conditions that will be used to evaluate the Poisson equation (20) are

Figure 2 shows the modification of spatial densities *n _{p}*,

*n*and electrostatic potential $\Phi $ obtained from numerical solution of Eq. (20) when varying the coil potential $ \Phi c$ while keeping the other physical parameters fixed. As one may expect, a negatively charged coil attracts the positron component while repelling electrons. The converse occurs for a positively charged coil. The other essential feature of Fig. 2 is the effect of the first adiabatic invariant

_{e}*μ*on the macroscopic equilibrium state of the system: both positrons and electrons tend to accumulate in regions of higher magnetic field strength. This feature combined with the dependence with respect to the coil potential implies that when $ \Phi c > 0$ positrons preferentially populate the interior region close to the center stack and enclosed by the coil, $ r < r box$, while electrons surround the coil surface (first row in Fig. 2). The opposite occurs when $ \Phi c > 0$ (third and fourth rows in Fig. 2). Inspection of Eq. (20) also suggests that if the thermodynamic parameters of the two species are equal, i.e.,

*A*=

_{p}*A*,

_{e}*N*=

_{p}*N*, $ \beta p = \beta e$, and $ \gamma p = \gamma e$, and $ \Phi c = 0$ on $ \u2202 \Omega c$, then, $\Phi $ vanishes throughout Ω while the spatial densities are identical,

_{e}*n*=

_{p}*n*. This case is shown in the second row of Fig. 2.

_{e}Taking the third row in Fig. 2 as reference case, Fig. 3 summarizes how spatial densities *n _{p}*,

*n*and electrostatic potential $\Phi $ change when physical parameters are modified. While equilibrium profiles are not too sensitive to changes in the chemical potentials $ \gamma p , \gamma e$ and in the inverse temperatures $ \beta p , \beta e$, abrupt changes occur when there is asymmetry in the number of protons and electrons. The bottom row of Fig. 3 shows the case $ N e = 2 \xd7 10 11 \u226b N p = 10 9$. Note that positrons are almost completely expelled from the interior region $ r < r box$ and form a radiation belt-like structure on the exterior side of the coil, while electrons are mostly found in proximity of the center stack located at

_{e}*r*=

*r*

_{0}.

### B. *μ*- $ J \u2225$ Equilibrium

*μ*

*T*of electromagnetic fluctuations must be longer than the period

_{f}*T*of a bounce oscillation as described by Eq. (29). An estimate of the bounce frequency

_{b}*ω*(and, thus, of

_{b}*T*) can be obtained with the aid of Eq. (45), which rests on the hypothesis $ \mu B \u226b e | \Phi |$ or $ \beta p \u2212 1 , \beta e \u2212 1 \u226b e | \Phi |$. Therefore, assuming these conditions to hold, the typical bounce frequency is

_{b}*μ*is also preserved. Again, experimental evidence from the RT-1 device suggests that the requirement (61) is usually fulfilled since $ \omega f < 10 5 \u2009 \u2009 Hz$ there.

^{25}

In addition to the physical parameters *A _{p}*,

*A*,

_{e}*N*,

_{p}*N*,

_{e}*β*,

_{p}*β*,

_{e}*γ*, and

_{p}*γ*already encountered for the

_{e}*μ*equilibrium, the

*μ*- $ J \u2225$ equilibrium described by the Poisson equation (47) includes the temperature anisotropies $ \alpha 0 p$ and $ \alpha 0 e$ and the Lagrange multipliers

*ζ*and

_{p}*ζ*associated with the conservation of the total bounce actions $ J p$ and $ J e$. In the following, we study a plasma with $ \alpha 0 p = \alpha 0 e = 1$. Notice that if the plasma is isotropic $ \alpha 0 p = \alpha 0 e = 1 / 2$ since perpendicular motion carries twice the degrees of freedom of parallel dynamics. Nevertheless, as long as $ \alpha 0 p$ and $ \alpha 0 e$ are constant, the spatial densities profiles remain unchanged, but they are only scaled by factors $ 1 / \alpha 0 p$ and $ 1 / \alpha 0 e$, respectively [recall Eq. (46)]. We shall also consider the limiting case $ \zeta p / \beta p \omega b \u226a 1$ and $ \zeta e / \beta e \omega b \u226a 1$ to simplify the evaluation of the integrals on the right-hand side of Eq. (47). Then, the relevant Poisson equation is given by Eq. (50), where $ \zeta p , \zeta e$ do not appear explicitly. Notice that physically this implies that the total bounce actions $ J p , J e$ are not preserved during the relaxation of the system. However, the effect of bounce dynamics is still felt by the ensemble through the term $ \omega b J \u2225 / 2$ appearing in the energies $ H \u0303 p , H \u0303 e$ and the bounce averaged phase space measure $ d \Pi \u0303$ of Eq. (40).

_{e}*μ*- $ J \u2225$ equilibria obtained by numerical solution of the Poisson equation (50) for different values of the coil potential $ \Phi c$ when the other physical parameters are kept constant. As in the case of

*μ*equilibria, a charged coil tends to push the species with charge of the same sign in the internal region $ r < r box$ while attracting the species with opposite charge. Next, observe that, as in the previous case, the conservation of the magnetic moment

*μ*results in the tendency of particles to accumulate in regions of high magnetic field strength. However, there is a peculiar feature of

*μ*- $ J \u2225$ equilibria: bounce dynamics squeezes spatial density profiles along the equatorial line

*z*= 0. As a result, radiation belt-like structures are formed on both sides of the coil. This behavior is a consequence of the term $ 1 / B 0$ appearing in the expression of the spatial densities

*n*,

_{p}*n*[recall Eq. (48)], which is related to the characteristic bounce velocity

_{e}*v*given by

_{b}*B*

_{0}, which is located at the radial position of the current loop $ r = r l$. Physically, the formation of radiation belt-like structures can be understood as follows. Particles with a high characteristic bounce velocity are statistically less likely to occur, and the net result is the balance between the tendency of particles to populate regions of high magnetic field strength

*B*as a result of the first adiabatic invariant, and the depletion effect caused by the second adiabatic invariant at those radial positions where

*B*

_{0}, and thus

*v*, are higher.

_{b}### C. *μ*– $ J \u2225$– $\Psi $ equilibrium

*μ*and $ J \u2225$. For the system under examination, the drift frequency is determined by the frequency of the toroidal rotation around the vertical axis. At equilibrium, the toroidal drift velocity $ v \phi = r v d \xb7 \u2207 \phi $ is given by the toroidal component of the guiding center drift velocity $ v d = v E + v k$, which is the sum of $ E \xd7 B$ drift

*β*and $ \alpha = E \u2225 / E \u22a5$ are the inverse temperature and the temperature anisotropy of the particle species under consideration. In deriving (64), we used the fact that in a vacuum magnetic field the term $ \u2202 \u2113 \xd7 \u2207 B$ and the curvature $ \u2202 \u2113 2$ are related by $ \u2202 \u2113 \xd7 ( \u2207 \xd7 B ) = \u2202 \u2113 \xd7 ( \u2207 B \xd7 \u2202 \u2113 ) \u2212 B \u2202 \u2113 2 = 0$, which implies $ \u2202 \u2113 \xd7 \u2207 B = B \u2202 \u2113 \xd7 \u2202 \u2113 2$. When the plasma temperature $ \beta \u2212 1$ and the magnetic field curvature $ \u2202 \u2113 2$ are sufficiently small, the $ E \xd7 B$ drift velocity becomes the dominant contribution to $ v d$. Then, the order of the drift frequency can be evaluated as

^{25}suggests that the frequency of electromagnetic fluctuations within the trap will be comparable to Eq. (65), implying that the criterion for the existence of a third adiabatic invariant $\Psi $,

*n*,

_{p}*n*and the electrostatic potential $\Phi $ obtained by numerical solution of Eq. (52) for different values of the chemical potentials $ \eta p , \eta e$ associated with the third adiabatic invariant $\Psi $. The effect of the preservation of total magnetic fluxes $ \Psi p , \Psi e$ is felt by the system only for large values of $ \eta p , \eta e$, corresponding to large changes in energy $ \beta p d E p = \eta p d \Psi p \u2248 10 , \u2009 \beta e d E e = \eta e d \Psi e \u2248 10$ when magnetic fluxes $ d \Psi p , d \Psi e \u2248 R 2 B \u2217 \u2248 1 \u2009 m 2 T$ are added to the system. The spatial densities are initially flattened at the equator (compare the spatial densities of first and second rows in Fig. 6), eventually splitting into upper and lower lobes separated by the equatorial line

_{e}*z*= 0 (compare the spatial densities of first and third rows in Fig. 6). Physically, this behavior can be understood as follows. The kinetic energy associated with toroidal drift motion can be written as

This term does not appear in the guiding center Hamiltonians (6) because it is usually smaller than the kinetic energies associated with parallel and cyclotron dynamics. In particular, this implies that $ p \phi \u2248 q \Psi $, with $ p \phi $ the canonical momentum of a charged particle in the $\phi $ direction and where we used the fact that in a dipole magnetic field the vector potential is $ A = \Psi \u2207 \phi $. Since states with large deviations in the canonical momentum $ p \phi \u2248 q \Psi $ would break the conservation of the total magnetic flux $ \Psi tot$, regions with large $\Psi $ are penalized in the distribution functions through the factors $ e \u2212 \eta p \Psi $ and $ e \u2212 \eta e \Psi $. If the chemical potentials $ \eta p , \eta e$ are sufficiently large, this effect prevails on the tendency caused by the first adiabatic invariant *μ* to concentrate particles in regions with strong magnetic field strength *B*, resulting in a preferential depletion of the equatorial region outside the coil (recall that along magnetic field lines, which correspond to level sets of $\Psi $, the magnetic field strength increases when approaching the coil from the outside).

## VI. TOROIDAL ROTATION

*μ*does not affect the profile of the spatial density [recall Eq. (18)]. The only relevant constraint is, thus, that given by the third adiabatic invariant $\Psi $. Indeed, denoting with $ v = v \u2225 + v E$ the velocity of a charged particle (other drifts are absent due to the homogeneity of the magnetic field), one has

*f*is the particle distribution function and

*N*is the particle number, then leads to a spatial density distribution

*A*is a positive real constant and

*η*is a Lagrange multiplier associated with the preservation of total magnetic flux. The Poisson equation for the electrostatic potential $\Phi $ in an infinite vertically symmetric $ \u2202 \Phi / \u2202 z = 0$ plasma column now reads

*z*-axis with velocity

^{23,46}We remark that, however, we are not aware of positron–electron experiments in this context.

*μ*, which forces particles in regions of high magnetic field strength, and not by the magnetic flux $\Psi $ as in a straight magnetic field. As shown in Sec. II, this also implies that, in contrast with a straight magnetic field, a dipole magnetic field is suitable to trap both neutral and nonneutral plasmas. Notice also that the toroidal rotation velocity (68) now depends on both the electrostatic potential $\Phi $ and the underlying magnetic field geometry through $ \u2207 \Psi , \u2009 \u2207 \u2113$, and $ \u2202 \u2113 2$. Figure 7 shows the toroidal rotation velocity obtained by numerical evaluation of Eq. (68) for the case reported in the second row of Fig. 4. Both the positron toroidal velocity $ v \phi p$ and the electron toroidal velocity $ v \phi e$ increase with

*r*but have opposite directions. This implies that, for the case considered, the charge-dependent drift $ v k$ is dominant with respect to $ v E$. Furthermore, there is a net toroidal current density $ J \phi = e ( n p v \phi p \u2212 n e v \phi e )$. Nevertheless, the resulting magnetic field $ B \u2032$ is negligible with respect to the dipole magnetic field. Indeed,

*μ*

_{0}is the vacuum permeability and the characteristic values $ R = 1 \u2009 m , \u2009 n \u2248 10 11 \u2009 m \u2212 3$, and $ v \phi \u2248 10 5 \u2009 ms \u2212 1$ were used.

## VII. CONCLUDING REMARKS

In this work, we studied the maximum entropy states of a collisionless positron–electron plasma trapped by a dipole magnetic field with the aim of elucidating the confinement properties of the system for different ranges of physical parameters. Such dipole magnetic field trap has several potential applications, including containment of pair and antimatter plasmas, experimental investigation of exotic and astrophysical plasmas, as well as technology development such as realization of coherent gamma ray lasers.

In a dipole magnetic field, the nature of plasma equilibria depends on the presence of adiabatic invariants. For such conserved quantities to hold, the timescale of electromagnetic fluctuations affecting the energy of a charged particle must be longer than the timescale of the periodic motion associated with each adiabatic invariant. Each adiabatic invariant constrains the maximum entropy state, resulting in a departure from standard Maxwell-Boltzmann statistics of an ideal gas.

Compared to a plasma trap with a straight magnetic field (such as a Penning–Malmberg trap) where radial confinement is provided by the conservation of the fragile third adiabatic invariant (the canonical momentum $ p \phi \u2248 q \Psi $), in a dipole magnetic field containment is realized through the first adiabatic invariant (the magnetic moment *μ*), which results in a tendency of each charged species to move toward regions of high magnetic field strength *B*. For this reason, a dipole magnetic field is suitable to confine both neutral and nonneutral plasmas. The effect of the second adiabatic invariant (the bounce action $ J \u2225$) is to squeeze the spatial densities of both positrons and electrons along the equatorial line, with the formation of characteristic radiation belt-like structures around the coil.

By solving Poisson's equation for the electrostatic potential with the charge density obtained from the maximum entropy states as source term, we put the theoretical model to the test and showed efficient confinement of both species for a wide range of physical parameters. The equilibrium profiles appear to be mostly sensible to asymmetries in the number densities of the two species and to significant changes in the coil potential. This latter fact suggests that the capability to control the coil potential would be a desirable property of any experimental dipole magnetic field trap design.

## ACKNOWLEDGMENTS

The research of NS was partially supported by JSPS KAKENHI Grant Nos. 21K13851 and 22H04936. The author is grateful to H. Saitoh, who made helpful suggestions and critiqued a preliminary draft of the paper.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Naoki Sato:** Conceptualization (lead); Formal analysis (lead); Writing – original draft (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## REFERENCES

*Mechanics*

^{3}He fusion reactor based on a dipole magnetic field

*The Framework of Plasma Physics*

*Introduction to Plasma Physics*