On a plasma sheath with a small normal magnetic field separating regions of oppositely directed magnetic field

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,¹ while the current sheet of the Earth's magnetotail was encountered by the early Imp 1 mission.² Even before the first space observations, current sheets were diagnosed in laboratory pinch experiments³ and motivated Dr. Harris to report on the now-renowned Harris model for 1D current sheets.⁴ 

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,¹⁷ with other projections shown in Figs. 1(d) and 1(e). Irrespective of how small $| B z |$ may be, any finite value of B_z will result in the described free streaming of particles across the sheet, not respected by the Harris solution.

The paper is organized as follows: For fixed profiles of $B x ( z )$ and $Φ ( 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 $Φ ( 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 $Φ ( z )$-profiles. The paper is concluded in Sec. V.

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.²⁴ For example, for the considered 1D geometry, within the regions where the magnetic moment $μ = m v ⊥ 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 σ = f σ 0 ( U , μ )$⁠. However, a solution of this form is not applicable to the center of a current sheet where μ breaks as an adiabatic invariant.

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 ⊥ 0$⁠, parallel and perpendicular to the asymptotic magnetic field, respectively. In the asymptotic plasma, the electrons and ions share the common number density n₀.

The field-reverse-configuration (FRC) considered for magnetic confinement of fusion plasma²⁶ 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 ( − 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₀ 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 = − B 0 / 10$⁠.

In our analysis, we also include a finite electric field, E_z. Initially, we focus on the case where E_z is directed toward z = 0, illustrated in Fig. 2(a). Thus, the electrostatic potential, $Φ = − ∫ 0 z E z d z ′$⁠, 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 ∥ 0 / T e ⊥ 0 = 4$⁠, and with isotropic ions characterized by $T i = T e | | 0$⁠.

The color contours in Fig. 2(e) illustrate $Φ 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.

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_T). As an example, in Fig. 2(g), the color contours of constant $J z$ are presented in the (x_T, z_T)-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 $Φ eff$ increases with the distance away from the A_y = 0 contour such that arbitrary large values of U_yz can be accommodated by considering a sufficiently large domain size.

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 $Φ ( 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_z) for $z / δ ∈ { − 7 , − 2 , 0 }$⁠. Here, the thermal electrons at $z = − 7 δ$ 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 = μ$⁠. For $z = − 2 δ$ considered in Fig. 3(b), electrons with $v y / v 0 ≲ − 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 = 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_z into the negative x-direction. Thus, for meandering orbits (all characterized by $v y < 0$⁠), a reduced fraction of v_y is available to be deflected into the 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_y outside the current sheet.

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₀ be anisotropic. We therefore explore solutions to Eq. (4) where f₀ is a bi-Maxwellian characterized by $T e | | 0 > T e ⊥ 0$⁠.

We first discuss how the self-consistent profiles of $B x ( z )$ can be determined for $T i ≪ T e ∥ 0$⁠. For this case, the ions respond very strongly to $Φ ( z )$⁠, and quasi-neutrality then requires that $e Φ ≃ T i ≪ T e ∥ 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₀ 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 ) = − e ∫ v y f k d 3 v$ is evaluated and the matching magnetic field is then readily computed $B x J y k ( z ) = μ 0 ∫ 0 z J y k ( z ′ ) d z ′$⁠. 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_z 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 ∫ ( v i − ⟨ v i ⟩ ) 2 f e ( v ) d 3 v$⁠, with $i ∈ { x , y , 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_ezz shown in Fig. 4(e). Consistent with force balance across the current layer, we observe that $P ezz + B x 2 / ( 2 μ 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 z-direction across the current sheet.

Numerical solutions to the 1D current sheet problem are obtained for a matrix spanned by three separate values of $T e ∥ 0 / T e ⊥ 0$ and three separate values of $T i / T e ∥ 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 ∥ 0 ∈ { 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 ∥ 0 / T e ⊥ 0 ∈ { 3 / 2 , 3 , 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.

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 $Φ ( 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 $Φ ( z ) ≠ 0$⁠.

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_z magnetic field. For this new class of E_z-trapped electrons, they all cross the 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_z-trapped region splits the yellow regions of the cyclotron orbit (which before were continuous across $z = 0 )$⁠.

Similar to the analysis in Sec. II, for the present configuration, the orbit topologies can be understood by analyzing $Φ eff$ introduced in Eq. (16). A 2D profile of $Φ 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 / δ = 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 $Φ eff$ marked by the red and black circles, respectively, can be found as the roots of $∂ Φ eff / ∂ z = 0$⁠. As a function of x_T, we denote the roots with z < 0 as $z l 1$⁠, defining the line marked l₁ in Figs. 7(c) and 7(d). In accordance with Fig. 7(a), the line marked l₂ in Figs. 7(c) and 7(d) is then obtained by identifying the value, $z l 2$⁠, for which $Φ eff ( x T , z l 2 ) = Φ 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_T)-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 π ) ∮ 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 l₁ and l₂, 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₁ and l₂ 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_T)-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 (x_T, z_T)-plane, the contours now do not represent the actual values of $J z$⁠. Rather, the contours provide a mapping between the orbit coordinates (x_T, z_T) 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).

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 $Φ eff$ as a function of x_T for $z T ≃ 0$⁠. The confined electrons are being reflected in this well before reaching the bifurcation line, l₁. As displayed in Fig. 8(b), we may evaluate $Φ eff$ along the l₁ line as a function of $J z , mea$⁠. This profile, $Φ eff l 1 ( J z , mea )$⁠, has a minimum at $J z , mea *$⁠, where for the present configuration, $log 10 ( J z , mea * ) ≃ − 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 Φ eff ( x T , z ) + m e ( v x 2 + v z 2 ) / 2$ must be less than $Φ eff l 1 A ( J z , mea )$⁠, displayed in Fig. 8(b).

To obtain (x_T, z_T) 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 / δ ∈ { − 0.3 , − 0.15 , 0 }$⁠, respectively. The location $z / δ = − 0.3$ of Fig. 9(a) is outside the region of local trapping by E_z, 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 ≃ 0$ and $v z ≃ 0$⁠, large values of $J z$ are observed corresponding to E_z-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 ⊥ 0 = 4$⁠. For the location at the center of the current sheet, the enhanced values of $J z$ for $v z ≃ 0$ causes significant reductions in f, yielding the observed double peaked distributions in the (v_y, v_z)-plane.

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.²⁹ 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 $β e ∞ = 2 − 4$⁠. The run was explored in Ref. 28, where a more complete account of the simulation setup is given.

The full profile of $E z ( x , z )$ from the kinetic simulation is shown in Fig. 10(a), normalized by $E 0 = E δ / ( e δ )$⁠, where $E δ$ is evaluated using Eq. (8). As is typical for simulations of anti-parallel reconnection,³² 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₀ in Eq. (11). A profile cut at x = 0 is shown in Fig. 10(d), where a bifurcated double peaked structure in J_e is clearly observed. Meanwhile, the similar cut across the number density profile in Fig. 10(e) is mostly uniform.

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 ) = ∬ f ( v ) 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 ⊥ 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_z 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.^28,30,31 Future work will explore in more detail the applicability of the 1D model to anti-parallel reconnection.

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).⁴ However, through Eq. (12), we notice how p_y becomes a function of x for any finite B_z. This x-dependency of p_y imposes that any 1D solution independent of x must be also independent of p_y. With these observations, it is then clear that the Harris-type solution is only valid for the case where B_z is exactly zero.

Breaking with the Harris-approach of parameterizing the solution in terms of the canonical momentum, p_y, we instead use that p_y is a constant of motion variable to derive closed form expressions for the $J z ∝ ∮ 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, f₀ 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 $Φ ( 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 $Φ ( 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_z, it is clear that 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₀. 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 ⊥ 0$⁠, the solutions can be parameterized in terms of the temperature ratios $T i 0 / T e | | 0$ and $T e ⊥ 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 ⊥ 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 $Φ ( z )$ are monotonic in $| z |$⁠. However, from reconnection studies, it is known that electron current layers can develop including a non-monotonic profile of $Φ ( | 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 $Φ ( z )$⁠, which yield electrically trapped electron trajectories. By applying profiles of $Φ ( 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 / ω p e$⁠.