The Harris-sheet model provides an elegant solution to the kinetic plasma equation for a steady state 1D current sheet geometry separating regions with oppositely directed magnetic field. However, adding just a small normal magnetic field to the Harris configuration yields thermal streaming of particles into and out of the current sheet, fundamentally changing the form of its kinetic description. The action variable, $ J z$, associated with the oscillatory orbit motion perpendicular to the current sheet is well conserved and can be applied for solving the kinetic equation in the 1D sheet geometry that includes a small normal magnetic field. Revisiting this problem, we develop a new formalism that permits numerical solutions to be readily obtained for general upstream/asymptotic electron and ion distributions. In particular, we consider the case of isotropic ion pressure and anisotropic bi-Maxwellian electrons. The current sheets are then supported by electron pressure anisotropy. Furthermore, the total current across a particular sheet is set by the fire-hose condition based on the electron pressures normalized by the asymptotic magnetic field pressure. Analytical approximations are obtained for the numerical solutions expressed in terms of the asymptotic electron temperature anisotropy and the ion temperature. We discuss a preliminary application of the framework to the electron diffusion region of anti-parallel magnetic reconnection.

## I. INTRODUCTION

Plasma current sheets have been observed in nearly all regions of space accessed *in situ* by spacecraft. For example, the ongoing Parker Solar Probe mission has recently provided direct observations of the heliopsheric current sheet at locations less than 0.2 AU from the Sun,^{1} while the current sheet of the Earth's magnetotail was encountered by the early Imp 1 mission.^{2} Even before the first space observations, current sheets were diagnosed in laboratory pinch experiments^{3} and motivated Dr. Harris to report on the now-renowned Harris model for 1D current sheets.^{4}

^{5}current sheet stability,

^{6}and magnetic reconnection.

^{7}It provides an exact closed-form analytical solution to the full Vlasov–Maxwell system of equations, constituting a kinetic description of a 1D plasma current sheet separating oppositely directed magnetic fields. As an illustration of the geometry, the Harris profiles of the current density and magnetic field are shown in Figs. 1(a) and 1(b) corresponding to the analytical solutions,

*B*

_{0}is the asymptotic magnetic field strength that reverses its direction across the sheet,

*μ*

_{0}is the vacuum permeability, and sech $ ( x ) = 1 / cosh ( x )$. The thickness of the sheet is set by

*δ*.

^{8}), it can be important to include a normal magnetic field,

*B*. Here, the normal magnetic field has been associated with double peaked current sheets, and its presence strongly influences the onset of magnetic reconnection.

_{z}^{9–13}The strong impact of a finite

*B*is perhaps not surprising as it changes fundamental properties of the system. For example, when imposing a non-zero

_{z}*B*onto the Harris solution, an unbalance force appears in the x-direction, $ F x = J y B z$, violating the equilibrium condition. In fact, earlier work

_{z}^{14–16}show that 1D current layers (including a normal magnetic field) can only be in equilibrium if the asymptotic plasma has an anisotropic pressure characterized by the marginal fire-hose condition,

Another complication to the Harris solution is evident in Fig. 1(c), where the black lines represent projections of the magnetic field lines onto the *xz*-plane. A typical trajectory of a thermal electron is illustrated by the blue line. Due to the finite *B _{z}* component, the electrons stream along the magnetic field into the center of the current sheet, where they reflect and are again ejected. This type of trajectory is known as a Speiser orbit,

^{17}with other projections shown in Figs. 1(d) and 1(e). Irrespective of how small $ | B z |$ may be, any finite value of

*B*will result in the described free streaming of particles across the sheet, not respected by the Harris solution.

_{z}*z*-oscillations of the Speiser orbit, the associated action integral $ J z \u221d \u222e v z d z$ is an adiabatic invariant,

^{18,19}provided that the normal magnetic field is small.

^{20}This current sheet invariant has been explored for solving the kinetic equation of the 1D current sheet problem, by writing the phase-space distributions in the form

*U*is the total particle energy, $ \Phi ( z )$ is the electrostatic potential, $ f \sigma 0$ is the asymptotic phase-space distribution, while $ q \sigma $ is the charge for particle species

*σ*. In the present paper, we will also study solutions provided by Eq. (4); these results are complementary and generally in agreement with earlier work

^{21–23}on this problem.

The paper is organized as follows: For fixed profiles of $ B x ( z )$ and $ \Phi ( z )$, in Sec. II, we discuss in more detail the current sheet adiabatic invariant, $ J z$, as well as our general approach for obtaining solutions to Eq. (4). In Sec. III, a numerical scheme is detailed by which self-consistent profiles of $ B x ( z )$ and $ \Phi ( z )$ are determined and analytical approximations for these profiles are derived. Relevant to kinetic simulations of reconnection, Sec. IV includes an analysis of modifications to the numerical solutions caused by non-monotonic electrostatic $ \Phi ( z )$-profiles. The paper is concluded in Sec. V.

## II. SOLVING THE KINETIC EQUATION USING THE $ J z$-ACTION VARIABLE

### A. The $ J z$-action variable, and its relation to the magnetic moment

The present type of kinetic problem involving free streaming of particles along some variable magnetic field can in many cases be solved using the drift-kinetic framework.^{24} For example, for the considered 1D geometry, within the regions where the magnetic moment $ \mu = m v \u22a5 2 / ( 2 B )$ is an adiabatic invariant, a steady state drift-kinetic solution to the free-streaming particles can in certain cases be written as $ f \sigma = f \sigma 0 ( U , \mu )$. However, a solution of this form is not applicable to the center of a current sheet where *μ* breaks as an adiabatic invariant.

*μ*, we recall that its invariance is directly tied to its properties as an action integral. In terms of the gyrophase $\varphi $ and the Larmor radius $ \rho L = m v \u22a5 / ( e B )$, this action integral takes the form

*μ*. In this outside region, the

*v*and

_{y}*v*contributions on the right of Eq. (5) are identical. The breakdown of

_{z}*μ*occurs abruptly as the electrons transition into the Speiser orbits. For the sections of meandering orbit motion across the current layer,

*v*does not change sign and there is no periodic motion in the

_{y}*y*-direction. Therefore, the $ \u222e v y d y$-integral in Eq. (5) cannot be computed, and

*μ*is undefined.

*z*-direction, as it is evident from Figs. 1(c) and 1(d) that orbits in the current layer geometry include well-defined

*z*-oscillations. The action integral $ \u222e v z d z$ is an adiabatic invariant,

^{19,20}but it should be noted that the transition from the regular cyclotron motion into the meandering motion can be understood as a merger at

*z*= 0 between two cyclotron orbits. Thus, during this merger, the $ \u222e v z d z$-integral doubles.

^{20}Based on this observation, we normalize this $ J z$-action by the expression

*ρ*is much larger than the radius of curvature

_{L}*R*of the magnetic field for

_{c}*z*= 0. The curvature parameter $ \kappa = R c / \rho L \u226a 1$ is then small and characterizes the quasi-adiabatic regime where the change in $ J z$ as a particle travels in and out of current sheet is on the order of $ \kappa J z$.

^{20}Furthermore, with our choice for normalization in Eq. (6), it also follows that in the external regions (where $ \u222e v y d y = \u222e v z d z$), we simply have $ J z = \mu $.

^{24}but includes some important differences. Contrary to the drift-kinetic approach where $ \mu = m e v \u22a5 2 / ( 2 B )$ is only the first term in an expansion of the full adiabatic invariant,

^{25}here we will evaluate $ J z ( z , v y , v z )$ without any approximation. As a result, Eq. (4) provides a steady state solution to the full kinetic equation [ $ \u2202 f \sigma / \u2202 t + \u2207 x \xb7 v f \sigma + ( q / m ) \u2207 v \xb7 ( E + v \xd7 B ) f \sigma = 0$].

While our approach is general, we limit our investigations to the case where the ions are characterized by isotropic pressure with an uniform temperature *T _{i}* across the geometry. Meanwhile, the electrons are allowed to be anisotropic, characterized by asymptotic temperatures $ T e | | 0$ and $ T e \u22a5 0$, parallel and perpendicular to the asymptotic magnetic field, respectively. In the asymptotic plasma, the electrons and ions share the common number density

*n*

_{0}.

*B*

_{0}and

*δ*represent the free parameters in the Harris sheet solution. A characteristic electron energy related to the width of the current layer may here be defined as $ E \delta = ( e B 0 \delta ) 2 / ( 2 m e )$. Meanwhile, later it will become clear that the fields of the self-consistent solutions for a particular set of the asymptotic temperatures can be approximated through the Harris profiles in Eqs. (1) and (2) when using $ E \delta $ in Eq. (8) in combination with the expressions for

*B*

_{0}and

*δ*in Eqs. (9) and (10), respectively. Characteristic values of the current density and the $ J z$-action variable are defined in Eq. (11).

### B. Determining $ J z$ as a function of $ ( z , v y , v z )$

The field-reverse-configuration (FRC) considered for magnetic confinement of fusion plasma^{26} turns out to have orbit dynamics very similar to the 1D current sheets considered here. As such, the analysis below for electrons in the 1D current sheet closely follows the approach developed in Ref. 27 for analyzing the fast ion behavior in the FRC. For the present analysis, we choose coordinates such that the infinite current sheet $ J ( z ) e y$ flows in the *y*-direction and is symmetrical about *z* = 0 such that $ J ( \u2212 z ) = J ( z )$. Again, the geometry is illustrated in Fig. 1(a), which by Ampère's law yields the magnetic field reversal in the asymptotic magnetic field *B*_{0} in the *x*-direction illustrated in Fig. 1(b). Throughout the analysis, we include a normal magnetic field such that $ B = [ B x ( z ) , 0 , B z ]$. It is assumed that $ | B z |$ is small, and while the final results of the analysis turn out to be independent of *B _{z}*, all illustrations are obtained with $ B z = \u2212 B 0 / 10$.

In our analysis, we also include a finite electric field, *E _{z}*. Initially, we focus on the case where

*E*is directed toward

_{z}*z*= 0, illustrated in Fig. 2(a). Thus, the electrostatic potential, $ \Phi = \u2212 \u222b 0 z E z d z \u2032$, has the form shown in Fig. 2(b), symmetric about

*z*= 0 with a single global minimum for

*z*= 0. As will be discussed later, the profiles in Fig. 2 are self-consistently obtained for an upstream bi-Maxwellian electron distribution with $ T e \u2225 0 / T e \u22a5 0 = 4$, and with isotropic ions characterized by $ T i = T e | | 0$.

*B*= 0) is fully described by the

_{y}*A*-component of the magnetic vector potential,

_{y}*y*. Given this symmetry in the

*y*-direction, the canonical momentum,

*v*as a function of (

_{y}*x*,

*z*) can then be obtained by analyzing $ p y ( x , z , v y ) = p y 0$. Meanwhile, we notice that

*x*only appears in the expression for

*p*through the linear term

_{y}*xB*in Eq. (12), and we can therefore always pick a coordinate system where $ p y 0 = 0$. To be more specific, with $ x T = x \u2212 p y 0 / ( e B z )$, the constraint $ p y ( x , z , v y ) = p y 0$ can be expressed as

_{z}*x*, the corresponding electron is “tied” to the magnetic field line that passes through $ ( x T , z ) = ( 0 , 0 )$. To illustrate this property, the orbit shown in Fig. 1(c) was computed with $ p y 0 = 0$, such that

_{T}*x*=

*x*.

_{T}*z*-direction, such that this motion can be evaluated at a fixed

*x*and a fixed total perpendicular energy given by

_{T}*x*only depends on

_{T}*z*,

*x*. The profile in Fig. 2(c) is evaluated for $ x T / \delta = \u2212 20$ and displays a local maximum of $ \Phi eff$ for

_{T}*z*= 0. Corresponding to the regular cyclotron motion around the $ A y ( x T , z ) = 0$ field line, electrons with

*U*below the local maximum of $ \Phi eff$ will be trapped in one of the local wells. In contrast, electrons with larger

_{yz}*U*are characterized by meandering orbits across the current sheet. In Fig. 2(d), we show $ \Phi eff ( z )$ evaluated for $ x T / \delta = 7$, where the profile includes a single global minimum at

_{yz}*z*= 0. Generally, for positive values of

*x*, all electron orbits are of the meandering type.

_{T}The color contours in Fig. 2(e) illustrate $ \Phi eff$ as a function of $ ( x T , z )$. The cyan line is again the field line characterized by $ A y ( x T , z ) = 0$, and the cuts considered in panels (c) and (d) are marked by the black dashed lines. The transition points identified in panel (c) fall on the magenta line, marking the bifurcation between cyclotron and meandering orbits. As such, the line encloses the regions of local cyclotron motion; any electron bouncing to the left of the magenta line performs regular cyclotron motion, whereas electrons that reach regions to the right of the line follow the meandering motion.

*z*at a fixed

*x*, let

_{T}*z*represent the turning points where

_{T}*v*= 0. The bounce motion perpendicular to the

_{z}*x*-direction is then fully parameterized by (

*x*,

_{T}*z*), as the perpendicular energy in Eq. (15) is then known $ U y z = e \Phi eff ( x T , z T )$. In turn, for given values of (

_{T}*x*,

_{T}*z*), we may also compute

_{T}*v*as a function of

_{z}*z*,

*z*, and the color contours represent

_{T}*v*evaluated as a function of $ ( x T , z )$.

_{z}With Eq. (17), profiles like the one in Fig. 2(f) are readily evaluated for any *z _{T}*, and all information required for computing $ J z$ in Eq. (6) can therefore be obtained as a function of

*z*for any set of (

*x*,

_{T}*z*). As an example, in Fig. 2(g), the color contours of constant $ J z$ are presented in the (

_{T}*x*,

_{T}*z*)-plane. The white lines highlight a particular contour level of $ J z$ corresponding to the orbit shown by the red line. It is evident that $ J z$ is well conserved for this orbit as all of the orbit bounce points consistently fall on the white lines. We notice from Fig. 2(e) that $ \Phi eff$ increases with the distance away from the

_{T}*A*= 0 contour such that arbitrary large values of

_{y}*U*can be accommodated by considering a sufficiently large domain size.

_{yz}*x*,

_{T}*z*). Meanwhile, for evaluating the solution in Eq. (4), it is desirable to cast $ J z$ as a function of $ ( z , v y , v z )$. Here, from Eq. (12), it is clear that $ A y ( x T , z ) = A y ( 0 , z ) + x T B z$, and using this in Eq. (14), it then follows

_{T}*z*can be determined as the roots of

_{T}*x*and

_{T}*z*to be evaluated as a function of $ ( z , v y , v z )$. Based on these results and observations, we have implemented in MATLAB numerical routines that as a function of $ ( z , v y , v z )$ first solve for (

_{T}*x*,

_{T}*z*) using Eqs. (18) and (19). In turn, $ J z ( z , v y , v z )$ is then readily evaluated through interpolation in the type of $ J z ( x T , z T )$-map shown in Fig. 2(g).

_{T}### C. Solutions to the kinetic equation across a current sheet with $ B z \u2260 0$

Our analysis centers on the previous and well-established result that $ J z$ is conserved for particles entering and exiting the current sheet. Using the methods developed above for evaluating $ J z ( z , v y , v z )$, solutions for the electron distribution, $ f ( z , v )$ (where for convenience, we have dropped subscript *e*) of the form in Eq. (4) are then readily obtained for specified profiles of *B*(*z*) and $ \Phi ( z )$.

To further elucidate the approach, in Figs. 3(a)–3(c), color contours of constant $ J z$ are shown as a function of (*v _{y}*,

*v*) for $ z / \delta \u2208 { \u2212 7 , \u2212 2 \u2009 , 0}$. Here, the thermal electrons at $ z = \u2212 7 \delta $ are sufficiently far away from the current sheet that they are characterized by the regular perpendicular cyclotron motion. Consequently, in Figs. 3(a), the contours of constant $ J z$ are concentric circles centered on $ ( v y , v z ) = ( 0 , 0 )$. This result is characteristic for locations outside the current sheet where, as discussed above, $ J z = \mu $. For $ z = \u2212 2 \delta $ considered in Fig. 3(b), electrons with $ v y / v 0 \u2272 \u2212 1$ are influenced by the proximity of the current sheet yielding reduced values of $ J z$. This effect is more noticeable at the center of the current sheet (

_{z}*z*= 0) for which values of $ J z$ in Fig. 3(c) are reduced significantly for $ v y < 0$. Consistent with the trajectory shown in Figs. 1(c)–1(e), at the center of the current sheet, large speeds in the negative

*y*-direction will be deflected by the small

*B*into the negative

_{z}*x*-direction. Thus, for meandering orbits (all characterized by $ v y < 0$), a reduced fraction of

*v*is available to be deflected into the

_{y}*z*-direction during the orbit motion. As a consequence, the values of $ J z$ are reduced compared to orbits with similar negative values of

*v*outside the current sheet.

_{y}For the cases where the asymptotic distribution is isotropic, $ f 0 ( U , J z ) = f 0 ( U )$, the full solution $ f ( z , v ) = f 0 ( U , J z )$ remains isotropic for any *z* independent of $ J z$, and will carry zero current. To search for meaningful solutions consistent with a finite current density of the reversing magnetic field, it is required that *f*_{0} be anisotropic. We therefore explore solutions to Eq. (4) where *f*_{0} is a bi-Maxwellian characterized by $ T e | | 0 > T e \u22a5 0$.

*f*

_{0}with $ T e | | 0 / T e \u22a5 0 = 4$, in Figs. 3(d)–3(f), contours of constant $ f ( z , v )$ are evaluated for

*v*= 0 and $ z / \delta \u2208 { \u2212 7 , \u2212 2 , \u2009 0}$. Note that for the present choice of

_{x}*f*

_{0}, from Eq. (4), it follows that

*f*

_{0}. In contrast, for $ z / \delta = \u2212 2$ in panel (e), the region of reduced $ J z$ for $ v y < 0$ yields enhanced values of

*f*. This effect is most pronounced for the center of the current sheet shown in panel (f).

## III. THE SELF-CONSISTENT PROFILES OF *B*_{x} AND $ \Phi ( z )$

_{x}

### A. Profiles of $ B x ( z )$ for Φ = 0

We first discuss how the self-consistent profiles of $ B x ( z )$ can be determined for $ T i \u226a T e \u2225 0$. For this case, the ions respond very strongly to $ \Phi ( z )$, and quasi-neutrality then requires that $ e \Phi \u2243 T i \u226a T e \u2225 0$, such that Φ = 0 is a good approximation when solving for the electron distributions.

The described enhanced levels of *f* in Figs. 3(e) and 3(f) for $ v y < 0$ yield a net electron current in the *y*-direction, and we apply an iterative scheme to determine the self-consistent profile of $ B x ( z )$ that matches the current profiles carried by these electrons. At iteration step *k*, for a given *f*_{0} and profile of $ B x k ( z )$, we apply 200 evenly spaced values of *z* for which distributions $ f k ( z , v )$ of the types in Figs. 3(d)–3(f) are computed. From these distributions, $ J y k ( z ) = \u2212 e \u222b v y f k d 3 v$ is evaluated and the matching magnetic field is then readily computed $ B x J y k ( z ) = \mu 0 \u222b 0 z J y k ( z \u2032 ) d z \u2032$. As illustrated in Fig. 4(a), for the subsequent iterative steps, we use $ B x k + 1 ( z ) = [ B x k ( z ) + B x J y k ( z ) ] / 2$, where the first guess of $ B x 1 ( z )$ is shown by the black line. The solution quickly converges to reach its self-consistent form after just seven iterations. The corresponding profiles of the current density *J _{y}* and number density

*n*are shown in Figs. 4(b) and 4(c). Here, the final profile of

*n*displays a peak at the center of the current sheet, with an amplitude of about 150% of the asymptotic density.

To test the self-consistent properties of the converged solution in Figs. 4(d)–4(f), we display various profiles related to the temperature and pressure components varying across the current sheet. Because *v _{x}* and

*v*only appear in the model through $ v x 2$ and $ v z 2$ [like in Eq. (19)], the distributions are even in these. Thus, the model includes no off diagonal stress and only the elements on the diagonal of the pressure tensor are finite. As usual, these elements are defined by $ P eii = m e \u222b ( v i \u2212 \u27e8 v i \u27e9 ) 2 f e ( v ) d 3 v$, with $ i \u2208 { x , y , z}$.

_{z}The profile of $ T exx = P exx / n$ is constant, whereas $ T eyy = P eyy / n$ displays a substantial increase compared to the external value, imposed by the applied asymptotic distribution. The combination of the more modest increase in *T _{ezz}* and the peaked density profile produce the profile of

*P*shown in Fig. 4(e). Consistent with force balance across the current layer, we observe that $ P ezz + B x 2 / ( 2 \mu 0 )$ is constant across the current layer (see red line marked “sum”). This illustrates how the numerical approach provides kinetic solutions that fulfill the fluid force balance in the

_{ezz}*z*-direction across the current sheet.

*x*-direction. In our numerical scheme, the asymptotic value of

*B*is not prescribed. Rather, as shown in Fig. 4(a), during the course of the iterations, it quickly evolves toward the expected value characterized by

_{x}*B*

_{0}in Eq. (9). Thus, the iterative solution in Fig. 4(f) demonstrates how the profile of $ P e | | \u2212 P e \u22a5 \u2009 \u2009 [ = P exx \u2212 ( P eyy + P ezz ) / 2 ]$ coincides with $ B x 2 / \mu 0$ for $ z / \delta < \u2212 2.5$, and the calculation emphasizes how the pressure anisotropy in the asymptotic region is directly responsible of driving the total current across the current sheet. In fact, with $ K = \u222b \u2212 \u221e \u221e J y d z$ being the total current per unit length of the layer, from Ampère's law and Eq. (3), it follows that in the case of a 1D current sheet with isotropic ions, we generally have

### B. Profiles of $ B x ( z )$ and for $ \Phi ( z )$ parameterized by $ T e \u2225 0 / T e \u22a5 0$ and $ T i / T e \u2225 0$

*T*, we assume that the density of the ions follows the Boltzmann scaling,

_{i}*k*. For the next step,

*k*+ 1, we then use

Numerical solutions to the 1D current sheet problem are obtained for a matrix spanned by three separate values of $ T e \u2225 0 / T e \u22a5 0$ and three separate values of $ T i / T e \u2225 0$. In Fig. 5, the resulting profiles are displayed, where as indicated in panel (f), the blue, green, and red lines correspond to $ T i / T e \u2225 0 \u2208 { 0.1 , 1 , 10}$, respectively. As marked at the top of the columns, panels (a)–(f), (g)–(l), and (m)–(r) correspond to $ T e \u2225 0 / T e \u22a5 0 \u2208 { 3 / 2 , \u2009 3 , \u2009 6}$, respectively. Meanwhile, the black dashed lines are the result of curve fitting to the numerical profiles. In all cases, these approximate profiles are in good agreement with those evaluated numerically. In the following, we describe the approximate form obtained for each of the quantities.

*z*-coordinates are normalized by

*δ*as given in Eq. (10). From panels (a), (g), and (m), however, it is clear that the width of the current layers is also weakly sensitive to $ T i / T | | 0$ as well as $ T e | | 0 / T e \u22a5 0$ in a way not captured by Eq. (10). Through curve fitting, we obtain a more accurate approximation for the current layer width

*B*

_{0}are defined in Eqs. (8) and (9), respectively. The term $ n z = 0 / n 0$ represents the central density normalized by the density of the asymptotic plasma; the term can be estimated directly using the approximations that are given as follows.

*δ*given in Eq. (10), and not $ \delta *$ in Eq. (23).

*n*(

*z*) in Eq. (28), including the pressures of both the electron and ion fluids. From the black dashed curves in Figs. 5(e), 5(k), and 5(q), the validity of this approximation is verified. Finally, given the form for

*n*(

*z*) in Eq. (28), the profile of $\Phi $ is approximated in Eq. (29) using the Boltzmann density scaling imposed through Eq. (21) on the ions,

Again, all of the approximate profiles in Eqs. (23)–(29) are illustrated by the black dashed lines in Fig. 5 against the numerical profiles. In particular, the approximations for *n*(*z*) and $ \Phi ( z )$ were derived by assuming fluid force balance in the *z*-direction. Thus, the validation of these approximations demonstrates that the numerical profiles in fact honor fluid force balance also for the cases where $ \Phi ( z ) \u2260 0$.

## IV. NON-MONOTONIC $ \Phi ( z )$

### A. Globally trapped electrons by *E*_{z}

_{z}

In this section, we explore the effects of a non-monotonic electrostatic potential. As an example, the profile of *E _{z}* in Fig. 6(a) is now permitted to include multiple sign reversals as a function of

*z*, corresponding to the potential in Fig. 6(b). The inner electric fields pointing outward from

*z*= 0 can trap electrons yielding a new class of electron orbits of the type illustrated by the black electron trajectory in Fig. 6(c). In its

*xy*-projection shown in Fig. 6(d), the trajectory is dominated by a circular cyclotron motion around the

*B*magnetic field. For this new class of

_{z}*E*-trapped electrons, they all cross the

_{z}*z*= 0 plane, and we therefore also characterize these as meandering. In what follows, it will be detailed how the orbit topology map in Fig. 6(c) is obtained and how the

*E*-trapped region splits the yellow regions of the cyclotron orbit (which before were continuous across $ z = 0 )$.

_{z}Similar to the analysis in Sec. II, for the present configuration, the orbit topologies can be understood by analyzing $ \Phi eff$ introduced in Eq. (16). A 2D profile of $ \Phi eff ( x T , z T )$ is shown in Fig. 7(c), with the location of two representative cuts marked and displayed in Figs. 7(a) and 7(b), respectively. For $ x T / \delta = 7$ in Fig. 7(a), three wells are observed, where the central well confines the new class of meandering electrons. The two wells on either sides correspond to the regular cyclotron orbits. The local peaks of $ \Phi eff$ marked by the red and black circles, respectively, can be found as the roots of $ \u2202 \Phi eff / \u2202 z = 0$. As a function of *x _{T}*, we denote the roots with

*z*< 0 as $ z l 1$, defining the line marked

*l*

_{1}in Figs. 7(c) and 7(d). In accordance with Fig. 7(a), the line marked

*l*

_{2}in Figs. 7(c) and 7(d) is then obtained by identifying the value, $ z l 2$, for which $ \Phi eff ( x T , z l 2 ) = \Phi eff ( x T , z l 1 )$. By this procedure, the orbit topology map in Fig. 6(c) is readily obtained as a function of the turning point coordinates (

*x*,

_{T}*z*).

_{T}In Sec. II, we applied the normalization in Eq. (6) to ensure that $ J z$ changes smoothly across the bifurcation lines in the (*x _{T}*,

*z*)-plane. However, that normalization does not yield smooth changes for the present scenario including the electrically trapped electrons. Instead, as a first step in Fig. 7(d), we apply the common normalization of $ J z = e / ( 2 \pi ) \u222e v z d z$. In Fig. 7(d), it is observed that the region of electrically trapped meandering orbits (in between the cyclotron regions) are characterized by low values of $ J z$. Along the bifurcation lines

_{T}*l*

_{1}and

*l*

_{2}, we record the separate values, $ J z , mea$ and $ J z , cyc$, corresponding to the two orbit types. As a function of the position along the combined bifurcation line, $ J z , mea$ increases monotonically, whereas $ J z , cyc$ has a global minimum $ ( J z , cyc = 0 )$ where

*l*

_{1}and

*l*

_{2}meet. Given $ J z , mea$ varies monotonically, we may then consider $ J z , cyc$ as a function of $ J z , mea$, as shown in Fig. 8(a). This function is important because it informs how $ J z , mea$ of the meandering orbits maps onto $ J z , cyc$ of the cyclotron orbits where the electrons cross the bifurcation line.

In Fig. 8(c), we again evaluate color contours of constant $ J z$, but now with the modification that for the meandering part of the (*x _{T}*,

*z*)-plane, we instead display contours of $ J z , cyc [ J z , mea ( x T , z T ) ]$, where $ J z , cyc ( J z )$ is the function shown in Fig. 8(a). We hereby obtain a map of $ J z ( x T , z T )$, which is continuous across the bifurcation lines. Again, for the meandering part of the (

_{T}*x*,

_{T}*z*)-plane, the contours now do not represent the actual values of $ J z$. Rather, the contours provide a mapping between the orbit coordinates (

_{T}*x*,

_{T}*z*) to the value of $ J z$ a particular electron had as it first started its journey from the asymptotic distribution into the current sheet. A zoomed in view is provided in Fig. 8(d).

_{T}The red lines in Figs. 8(c) and 8(d) highlight particular contours of constant $ J z$. We observe how the turning points of the displayed orbits follow these contours. The blue orbit [only in Fig. 8(c)] is of the type explored in Sec. III [like the red orbit in Fig. 2(g)]. Meanwhile, the black orbit in Figs. 8(c) and 8(d) provides an example where an electron initially follows a regular cyclotron orbit and then transitions into the electrically confined orbit type. In turn, this is different from the confined cyan orbit, which does not reach the bifurcation lines and remains confined to the central region of the current sheet.

To understand the boundaries of the confined orbits, in Fig. 7(c), we notice the well in $ \Phi eff$ as a function of *x _{T}* for $ z T \u2243 0$. The confined electrons are being reflected in this well before reaching the bifurcation line,

*l*

_{1}. As displayed in Fig. 8(b), we may evaluate $ \Phi eff$ along the

*l*

_{1}line as a function of $ J z , mea$. This profile, $ \Phi eff l 1 ( J z , mea )$, has a minimum at $ J z , mea *$, where for the present configuration, $ \u2009 log 10 ( J z , mea * ) \u2243 \u2212 1$. As shown in Figs. 8(b)–8(d), the part of the bifurcation line with $ J z , mea < J z , mea *$ is labeled $ l 1 A$. For an electron to be globally confined it must have $ J z , mea < J z , mea *$, and its value of $ e \Phi eff ( x T , z ) + m e ( v x 2 + v z 2 ) / 2$ must be less than $ \Phi eff l 1 A ( J z , mea )$, displayed in Fig. 8(b).

### B. Electron distributions and bifurcated current sheets for non-monotonic $ \Phi ( z )$

To obtain (*x _{T}*,

*z*) as a function of $ ( z , v y , v z )$, we can again use Eqs. (18) and (19). In turn, with the map in Figs. 8(c) and 8(d), we may then evaluate $ J z ( z , v y , v z )$. The color contours $ J z ( v y , v z )$ in Figs. 9(a)–9(c) are obtained for $ z / \delta \u2208 { \u2212 0.3 , \u2212 0.15 , 0}$, respectively. The location $ z / \delta = \u2212 0.3$ of Fig. 9(a) is outside the region of local trapping by

_{T}*E*, and the profile of $ J z ( v y , v z )$ is thus similar to that observed in Fig. 3(c). Meanwhile, as shown in Figs. 9(b) and 9(c), for $ z \u2243 0$ and $ v z \u2243 0$, large values of $ J z$ are observed corresponding to

_{z}*E*-trapped electrons with trajectories similar to the black trajectory in Figs. 8(c) and 8(d). Using Eq. (4), we may again evaluate $ f ( z , v )$. The distributions in Figs. 9(d)–9(f) are obtained for an ambient distribution with $ T e | | 0 / T e \u22a5 0 = 4$. For the location at the center of the current sheet, the enhanced values of $ J z$ for $ v z \u2243 0$ causes significant reductions in

_{z}*f*, yielding the observed double peaked distributions in the (

*v*,

_{y}*v*)-plane.

_{z}### C. Similarity with current layers of anti-parallel magnetic reconnection

The 1D current layers explored in this section have characteristics relevant to magnetic reconnection. During magnetic reconnection, electron pressure anisotropy typically develops within the reconnection inflow regions.^{29} In particular, for anti-parallel reconnection, this upstream pressure anisotropy supports the current layer of meandering electrons that develops within the electron diffusion region (EDR).^{28,30,31} As an example, in Figs. 10(a)–10(f), we display various profiles obtained in a fully kinetic simulation. The data correspond to the case with an ion to electron mass ratio $ m i / m e = 1836$, and an upstream normalized electron pressure $ \beta e \u221e = 2 \u2212 4$. The run was explored in Ref. 28, where a more complete account of the simulation setup is given.

*δ*in (10) and the electron skin depth $ d e = m e / ( \mu 0 n e e 2 )$, used as the characteristic length scale in the kinetic simulation. From Fig. 1(b) of Ref. 28, it is seen that $ T e | | 0 / T e \u22a5 0 = 4$ just upstream of the EDR, which is useful for relating

*δ*and

*d*. In fact, using Eqs. (8)–(10), we find that

_{e}*δ*.

The full profile of $ E z ( x , z )$ from the kinetic simulation is shown in Fig. 10(a), normalized by $ E 0 = E \delta / ( e \delta )$, where $ E \delta $ is evaluated using Eq. (8). As is typical for simulations of anti-parallel reconnection,^{32} the cut of $ E z ( z )$ in Fig. 10(c) for *x* = 0 has a structure similar to that analyzed for the 1D current layer in Fig. 6(a). The total electron current density *J _{e}* is shown in Fig. 10(b), normalized by

*J*

_{0}in Eq. (11). A profile cut at

*x*= 0 is shown in Fig. 10(d), where a bifurcated double peaked structure in

*J*is clearly observed. Meanwhile, the similar cut across the number density profile in Fig. 10(e) is mostly uniform.

_{e}In the above analysis, we observed how applying the simulation *E _{z}*-profiles in the 1D model causes double peaked structures in the model phase space distribution. To help visualize these structures both for the simulation and the model, we compute the reduced phase-space distribution $ g ( v z ) = \u222c f ( v ) \u2009 d v x d v y$. In Fig. 10(f), we display $ g ( z , v z )$ for a cut across the EDR at

*x*= 0. A prominent electron–hole structure centered on

*z*= 0 is observed, similar to those studied in Ref. 32.

Predictions of the 1D current model are evaluated for $ T e | | 0 / T e \u22a5 0 = 4$, while imposing the *E _{z}*-profile shown in Fig. 10(g). Based on the iterative scheme detailed in Sec. III A, we obtain the profiles of $ J e y ( z )$,

*n*(

*z*), and $ g ( z , v z )$ shown in Figs. 10(h)–10(j), respectively. While some differences between the kinetic simulation profiles and those of the 1D current layer are clearly present, we notice the similarity in the profiles of $ J e y ( z )$ as well as $ g ( z , v z )$. It should be emphasized again that in the 1D model the results were obtained by imposing the profile of

*E*causing the bifurcated current layer and electron–hole structure. The self-consistent 1D current layers in Sec. III B did not develop these features, and it is possible that the bifurcated current layer is unique to the reconnecting geometry. However, the similarity between $ J e y ( z )$ of the simulation and the 1D model emphasizes how the reconnection current layer is mainly driven by the electron pressure anisotropy that develops in the reconnection inflow regions.

_{z}^{28,30,31}Future work will explore in more detail the applicability of the 1D model to anti-parallel reconnection.

## V. DISCUSSION AND CONCLUSION

In the present paper, we have developed a framework for understanding and modeling 1D current layers including a small normal magnetic field. First, it should be noted that the familiar Harris sheet solution is not compatible with any normal magnetic field, as such a field causes thermal streaming of particles across the current sheet. Mathematically, the underlying kinetic solution of the Harris sheet is parameterized in terms of *p _{y}* in Eq. (13).

^{4}However, through Eq. (12), we notice how

*p*becomes a function of

_{y}*x*for any finite

*B*. This

_{z}*x*-dependency of

*p*imposes that any 1D solution independent of

_{y}*x*must be also independent of

*p*. With these observations, it is then clear that the Harris-type solution is only valid for the case where

_{y}*B*is exactly zero.

_{z}Breaking with the Harris-approach of parameterizing the solution in terms of the canonical momentum, *p _{y}*, we instead use that

*p*is a constant of motion variable to derive closed form expressions for the $ J z \u221d \u222e v z d z$ action integral. In contrast to the magnetic moment, $ J z$ is a well-conserved adiabatic invariant. Based on this observation and with the development of numerical routines for evaluating $ J z ( z , v y , v z )$, we explore the class of 1D solutions simply expressed by the form $ f ( x , v ) = f 0 ( U , J z )$. Here,

_{y}*f*

_{0}is the asymptotic electron distribution, and

*U*is the total energy of the electrons. It should be noted that at bifurcation points where the orbit topology changes from the cyclotron type into the meandering type, the values of $ J z$ change abruptly. However, these changes are readily characterized and accounted for in the model.

Based on the framework outlined above, for given profiles of $ B x ( z )$ and $ \Phi ( z )$, the full solutions for $ f ( z , v )$ are readily obtained. This, in turn, facilitates a numerical scheme to determine the profiles of $ B x ( z )$ and $ \Phi ( z )$, consistent with Ampère's law and the condition of quasi-neutrality. While not imposed numerically, we find that the converged kinetic solutions are in fluid force balance across and along the current sheet, validating the numerical approach. Given the finite values of *J _{y}* and

*B*, it is clear that

_{z}*J*×

*B*-forces are present in the

*x*-direction along the current layer. In the self-consistent solution, these forces are balanced by thermal forces associated with the upstream pressure anisotropy. Concisely, the overall force balance in the

*x*-direction is expressed through the marginal fire-hose condition in Eq. (3). This equation is fundamental to the 1D current layer, as it parameterizes how the pressure anisotropy regulates

*B*

_{0}. In other words, it is the upstream pressure anisotropy that controls the value of the integrated current across the current layer. Our approach is similar to previous work in the literature,

^{22,23,33}but is perhaps more general. For example, in contrast to Ref. 22, we do not include any fluid closure assumptions for evaluating the electron current profile.

For the case of isotropic ions with temperature $ T i 0$ and bi-Maxwellian electrons with asymptotic parallel and perpendicular temperatures $ T e | | 0$ and $ T e \u22a5 0$, the solutions can be parameterized in terms of the temperature ratios $ T i 0 / T e | | 0$ and $ T e \u22a5 0 / T e | | 0$. No other parameters influence the mathematical form of the solutions. As a function of $ T i 0 / T e | | 0$ and $ T e \u22a5 0 / T e | | 0$, closed form numerical approximations are obtained across the current sheets of all essential fluid quantities. These approximations provide quantitative insight into how 1D current layers are affected by the asymptotic boundary conditions and may become useful to the analysis of spacecraft data or future theoretical studies of the present type of 1D current layer.

The numerically identified self-consistent profiles of $ \Phi ( z )$ are monotonic in $ | z |$. However, from reconnection studies, it is known that electron current layers can develop including a non-monotonic profile of $ \Phi ( | z | )$. We therefore generalize the 1D model to include such more general electrostatic potentials. In particular, kinetic simulations of anti-parallel magnetic reconnection include profiles of $ \Phi ( z )$, which yield electrically trapped electron trajectories. By applying profiles of $ \Phi ( z )$ of this nature in the 1D model, a range of features of the reconnection geometry are reproduced. Among these features is the creation of bifurcated double-peaked layers with a width on the order of the electron skin depth, $ c / \omega p e$.

Future studies will clarify the applicability of the 1D model to the structure of the current layers in anti-parallel reconnection. Preliminary results (not included here) suggest that although the inductive *E _{y}* reconnection electric field is not considered in our geometry, the approximations derived for $ J z$ remain valid. In fact, the new framework in combination with the electron distribution model in Ref. 29 appears to provide an accurate account of the

*J*current profile as well as the on-diagonal electron pressure components within the reconnection inflow and electron diffusion region. Furthermore, the model may be generalized to accommodate finite values of

_{y}*E*by imposing the previous observed rotation of the central reconnection current layer into the outflow directions. This, in turn, can account for the main off-diagonal electron stress components discussed in Ref. 28. The details of these findings will be reported in separate publications.

_{y}## ACKNOWLEDGMENTS

The author acknowledges discussions, comments, and suggestions by Young Dae Yoon, Samuel Greess, Abhishek Mhatre, and William Daughton. This work was supported in part by H. I. Romnes Faculty Fellowship by the UW-Madison Office of the Vice Chancellor for Research and Graduate Education.

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts to disclose.

### Author Contributions

**Jan Egedal:** Writing – original draft (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*MHD Description of Plasma*