Kinetic equilibrium of two-dimensional force-free current sheets

Force-free current sheets are local plasma structures with field-aligned electric currents and approximately uniform plasma pressures. Such structures, widely found throughout the heliosphere, are sites for plasma instabilities and magnetic reconnection, the growth rate of which is controlled by the structure's current sheet configuration. Despite the fact that many kinetic equilibrium models have been developed for one-dimensional (1D) force-free current sheets, their two-dimensional (2D) counterparts, which have a magnetic field component normal to the current sheets, have not received sufficient attention to date. Here, using particle-in-cell simulations, we search for such 2D force-free current sheets through relaxation from an initial, magnetohydrodynamic equilibrium. Kinetic equilibria are established toward the end of our simulations, thus demonstrating the existence of kinetic force-free current sheets. Although the system currents in the late equilibrium state remain field aligned as in the initial configuration, the velocity distribution functions of both ions and electrons systematically evolve from their initial drifting Maxwellians to their final time-stationary Vlasov state. The existence of 2D force-free current sheets at kinetic equilibrium necessitates future work in discovering additional integrals of motion of the system, constructing the kinetic distribution functions, and eventually investigating their stability properties.

The properties of magnetic reconnection and the dissipation rate of magnetic energy strongly depend on the current sheet configurations. The simplest magnetic field geometry is the one-dimensional (1D) current sheet, which is either a tangential discontinuity separating two plasmas with different properties or a rotational discontinuity having a finite magnetic field component, B n , normal to the current sheet. The stress balance in 1D tangential discontinuities is established by the gradients of plasma and magnetic field pressures [see Figure 1(a) and Allanson et al. (2015); Neukirch et al. (2020a,b)], whereas the stress balance in 1D rotational discontinuities requires a contribution from the plasma dynamic pressure in order to balance the magnetic field line tension force ∝ B n [see Figure 1(b) and Hudson (1970)].
In two-dimensional (2D) current sheets, B n ̸ = 0 naturally appears (Schindler 2006); thus, the stress balance requires 2D plasma pressure gradients [see Figure 1(c) and Yoon & Lui (2005)], 2D dynamic pressure gradients [see Figure  1(d) and Birn (1992); Nickeler & Wiegelmann (2010); Cicogna & Pegoraro (2015)], or pressure anisotropy (e.g. Sitnov & Merkin 2016;Artemyev et al. 2016b). The latter solution has multiple variants, most of which are devoted to configurations. The coordinate system consists of the l component along the main magnetic field direction, reversing sign at the current sheet neutral plane (B l = 0), the m component along the main current density direction corresponding to variations of B l , and the n (normal) component along the spatial gradient of B l (i.e., 4πjm/c = ∂B l /∂rn). Panel (e) shows the typical parameter regimes of current sheets observed by THEMIS and ARTEMIS (Angelopoulos 2008(Angelopoulos , 2011 in the solar wind (blue), Earth's magnetosheath (shocked solar wind; black), near-Earth magnetotail (red), and lunar-distance magnetotail (green). The black dashed box defines the parameter regime where 2D force-free current sheets are expected. Our dataset includes ∼ 300 solar wind current sheets, ∼ 100 magnetosheath current sheets, ∼ 100 current sheets in the near-Earth magnetotail, and ∼ 500 current sheets in the lunar-distance magnetotail. All current sheets were identified from THEMIS and ARTEMIS fluxgate magnetometer measurements (Auster et al. 2008). The selection procedure and data processing for solar wind current sheets are described in Artemyev et al. (2019a): ion temperatures are calculated using the OMNI dataset (King & Papitashvili 2005;Artemyev et al. 2018), plasma flows in the current sheet reference frame are calculated as the change of the most variable flow component across the current sheet (see details in Artemyev et al. 2019a). The same criteria and data processing methods were applied to the magnetosheath (dayside) current sheets, but without using the OMNI data because THEMIS accurately measures magnetosheath ion temperatures in that region (McFadden et al. 2008). The dataset of the near-Earth magnetotail current sheets (radial distances ∼ 10 − 30 Earth radii) is described in Artemyev et al. (2016a). The same criteria and data processing methods were applied to the lunar-distance current sheets. Details of the calculation of β and MA are described in Artemyev et al. (2019aArtemyev et al. ( , 2017. The class of such 2D force-free current sheets is not limited to observations in the solar wind and distant Earth's magnetotail, but also includes current sheets observed in the cold plasma of the Martian magnetotail (e.g., DiBraccio et al. 2015;Artemyev et al. 2017) and in the low-density Jovian magnetotail (e.g., Behannon et al. 1981;Artemyev et al. 2014). Figures 2(c,d) show examples of force-free current sheets observed in the Martian and Jovian magnetotails. There are not enough known integrals of motion to describe distribution functions of charged particles in such current sheet configurations (see discussion in Lukin et al. 2022). Consequently, the absence of kinetic models of force-free current sheets with B n ̸ = 0 significantly hinders the analysis of their stability, the process responsible for the magnetic reconnection onset and particle acceleration.
In this study, we develop a 2D kinetic force-free current sheet model. In Section 2, we describe a magnetohydrodynamic (MHD) model of 2D force-free current sheets. In Section 3, we initialize self-consistent, particle-in-cell (PIC) simulations using currents and magnetic fields from the MHD equilibrium. We load particle distributions using drifting Maxwellians, which initially do not necessarily satisfy the time-stationary Vlasov equation. In Section 4, we obtain kinetic equilibria by evolving these particle distribution functions using the self-consistent PIC simulations and  Figure 2. Four examples of force-free current sheets in (a) solar wind, (b) lunar-distance magnetotail, (c) Martian magnetotail, and (d) Jovian magnetotail. The top panels show the magnetic field in the local coordinate system Sonnerup & Cahill (1968). Grey curves show the magnetic field magnitudes; B ≈ constant indicates the force-free current sheets. Middle panels show electron and ion β profiles. Bottom panels show profiles of the field-aligned currents with an indication on the dominant ion species and estimations of the current sheet thickness in the ion inertial length. For the solar wind and Earth's lunar-distance magnetotail, we used ARTEMIS magnetic field (Auster et al. 2008) and plasma (McFadden et al. 2008) measurements (see details of data processing procedure in Artemyev et al. 2019a). For the Martian magnetotail, we used MAVEN magnetic field (Connerney et al. 2015) and plasma (McFadden et al. 2015;Halekas et al. 2015) measurements (see details of data processing procedure in Artemyev et al. 2017). For the Jovian magnetotail, we used Juno magnetic field (Connerney et al. 2017a,b) and plasma (McComas et al. 2017;Kim et al. 2020) measurements (see details of data processing procedure in Artemyev et al. 2023). demonstrate the existence of kinetic equilibria for 2D force-free current sheets. We describe the equilibrium distribution functions and the difference between ion-and electron-dominated 2D current sheets. In Section 5, we summarize the results and discuss their applications.

MHD EQUILIBRIA OF FORCE-FREE CURRENT SHEETS
We consider the general case of 2D force-free equilibrium in which all quantities vary with coordinates x and z only (i.e., ∂/∂y = 0), and, as is customary in this definition, the pressure gradient force contribution to the force balance is negligible, i.e., ∇P ∼ 0. The divergence-free magnetic field has the form (e.g., Low & Wolfson 1988) where Ae y is the vector potential in the y direction, and B y is the magnetic field in the y direction. The current density is where c is the speed of light, and ∇ 2 = ∂ 2 /∂x 2 + ∂ 2 /∂z 2 is the Laplacian. The force-free condition, J × B = 0, is equivalent to J = αB, where α = α(x, z) is a scalar function. Comparing Equations (1) and (2), the force-free condition requires that B y is a function of A only, and that B y satisfies and the coefficient α is (c/4π)dB y /dA. Once B y (A) is given, we can solve Equation (3) for A(x, z) with appropriate boundary conditions. As shown in Figure 3, for balancing force-free current sheets with low thermal and dynamic pressures (e.g., low β/low Mach number current sheets in the solar wind, magnetosheath, and lunar-distance magnetotail; see Figure 1), the magnetic pressure B 2 y plays an role analogous to the thermal pressure. Motivated by the functional form of pressure p ∝ exp [2A/(λB 0 )] describing thermal-pressure balanced current sheets (Lembege & Pellat 1982), we adopt for force-free current sheets, so that Equation (3) describing magnetic field lines in the x-z plane (i.e., contours of A) resembles that of the Lembege-Pellat current sheet (e.g., Tassi et al. 2008;Lukin et al. 2018): where λ is the current sheet thickness, and B 0 is the magnetic field at large |z|.
Near-Earth magnetotail Mixed states Solar wind; Lunar-distance magnetotail Figure 3. Force balance of current sheets in different space environments. At one end of the spectrum, representing current sheets in the near-Earth magnetotail, magnetic tension force JyBz/c and magnetic pressure force JyBx/c ≈ ∂(B 2 x /8π)/∂z are balanced by thermal pressure gradients ∂p/∂x and ∂p/∂z, respectively. At the other end of the spectrum, representing current sheets in the solar wind and lunar-distance magnetotail, JyBz/c and JyBx/c are balanced by shears of B 2 y , i.e., ∂(B 2 y /8π)/∂x and ∂(B 2 y /8π)/∂z, respectively. There may be a continuous spectrum of current sheets between these two ends (Yoon et al. 2023), in which JyBz/c and JyBx/c are balanced by a mixture of thermal pressure gradients and shears of B 2 y .
Because the scale length along current sheets is much larger than that across current sheets, we assume A = A(εx, z) to be weakly nonuniform in the x direction (ε being a small parameter). Equation (5) can be approximated as which is accurate up to order ε. The boundary condition is which is equivalent to B where F (x) = exp (−εx/λ). Thus the three components of the magnetic field are where B z ̸ = 0 describes rotational discontinuities. The current density is The direction of magnetic field rotates from pointing in +x at the z > 0 half-space (z/λ ≫ 1) to pointing in +y around the neutral (equatorial) plane (−1 ≲ z/λ ≲ 1), and further to pointing in −x at the z < 0 half-space (z/λ ≪ −1). The magnitude of magnetic field, (B 2 , has a weak dependence on x and is roughly a constant at a given x location. The magnitude of current density, (J 2 , is concentrated in a layer |z| < λF (x). The current density profile [Equation (10)], together with the constant plasma pressure profile, constrains the initial particle distribution functions.

COMPUTATIONAL SETUP
Searching for kinetic equilibria of 2D force-free current sheets, we initialize our simulated current sheets from the results of an MHD equilibrium. Importantly, this may not satisfy the time-stationary Vlasov equation because the initial particle distributions are simply drifting Maxwellians. We simulate their relaxation toward a certain kinetic equilibrium (if any) based on a massively parallel PIC code (Pritchett 2001(Pritchett , 2005. Our simulations have two dimensions (x, z) in configuration space and three dimensions (v x , v y , v z ) in velocity space. The results are presented in normalized units: magnetic fields are normalized to B 0 , lengths to the ion inertial length d i = c/(4πn 0 e 2 /m i ) 1/2 , time to the reciprocal of the ion gyrofrequency ω −1 Here n 0 is the plasma density, m i is the ion mass, and e is the elementary charge. The computational domain is ci . The ion-to-electron mass ratio is m i /m e = 100. The normalized speed of light is c/v A = 20, which gives the ratio of electron plasma frequency to electron gyrofrequency ω pe /ω ce = 2. The reference density n 0 is represented by 1929 particles. The total number of particles is 1.6×10 9 including ions and electrons in each simulation. Figure 4 shows the initial magnetic field configuration and current density, which are determined by Equations (9) and (10) using the input parameters ε = 0.04 and λ = 2d i . There are a few degrees of freedom in choosing the initial particle distribution functions: (1) While the current density is known a priori, the proportion of ion and electron currents is left unspecified; (2) Electron and ion temperatures and plasma β vary across different space environments (e.g., Figure 2); (3) Specific forms of the velocity distribution functions are left undetermined. To this end, we choose initial electron and ion distribution functions as drifting Maxwellians where the subscript s = i, e represents ions and electrons, v ds (x, z) is the position-dependent drift velocity, and T s is the temperature. Note that T i , T e , and n 0 are constants throughout the domain so that the initial plasma pressure n 0 (T i + T e ) is a constant, as required by the force-free condition. The simulations allow either ions or electrons to be the main current carrier, while keeping the relative drift between them commensurate with the current density [Figures 4(c), 4(d), and 4(e)]. Additionally, we vary particle temperatures to search for kinetic equilibrium in different plasma beta regimes, where the plasma beta is β = 2(T i + T e )/m i v 2 A . The detailed parameters for particle distribution functions in the six runs are listed in Table 1. These parameters cover the β and T i /T e ranges of force-free current sheets observed in the solar wind and planetary magnetotails (see Figure 2).
In our simulations, the electromagnetic fields are advanced in time by integrating Maxwell's equations using a leapfrog scheme. These fields are stored on the Yee grid. The relativistic equations of motion for ions and electrons are integrated in time using a leapfrog scheme with a standard Boris push for velocity update. The conservation of charge is ensured by applying a Poisson correction to the electric fields (Marder 1987;Langdon 1992).
For particles crossing the x boundaries, we take advantage of the symmetry between z > 0 and z < 0 [ Figure 4], so that particles exiting the system at a location z with velocity (v x , v y , v z ) are reinjected into the system at the conjugate location −z with velocity (−v x , v y , v z ). This is equivalent to an open boundary condition for particles because the injected particle distribution matches that at one cell interior to the boundary at all simulation times. At the z boundaries, particles striking the boundary are reflected into the system with v z = −v z .
For fields at the x boundaries, the guard values of the tangential magnetic fields are determined by δB n y,g = δB n y,i1 , δB n z,g = δB n−1 z,i1 (2 − ∆x/c∆t) + δB n−2 z,i2 (∆x/c∆t − 1),  where the superscript indicates the time level, and the subscripts g, i1, i2 indicate the guard point, first interior point and second interior point, respectively. This boundary condition for δB z ensures that the magnetic flux can freely cross the x boundaries (Pritchett 2005). The guard values of the normal electric field δE x at the x boundaries are determined by δE n x,g = δE n x,i1 . Similarly, at the z boundaries, the two components of the tangential magnetic fields in the guard point are determined by δB n x,g = δB n x,i1 and δB n y,g = δB n y,i1 ; The normal electric field δE z in the guard point are determined by δE n z,g = δE n z,i1 . Unlike the driven simulations that adds magnetic flux to the simulation box, no external E y is applied at the z boundaries.

RESULTS
We track the evolution of the initially force-free current sheets in the six runs up to t = 180 ω −1 ci , at which time the macroscopic states (e.g., electric and magnetic fields, currents, pressures) of the system are quasi-stationary. Below we examine if the configurations satisfy J × B = 0 at the end of the simulations. We further investigate how the particle velocity distributions deviate from the initial Maxwellians and whether or not the system reaches a kinetic equilibrium.

Initially ion-dominated current sheets
When ions initially carry entirely the (field-aligned) currents (Runs 1A, 2A, and 3A), electrons are accelerated by transient parallel electric fields to form field-aligned currents [e.g., Figure 5(b)], while ions are slightly decelerated by such electric fields [e.g., comparing Figures 5(a) and 4(e)]. These electron currents can be a significant fraction (e.g., ≳ 1/3) of ion currents in the off-equatorial region (1 < |z| < 5). Almost all currents remain field aligned [see J ⊥ ∼ 0 in Figure 5(c)], implying that a force-free configuration may indeed be obtained. During the relaxation of ion-dominated current sheets, electrostatic fields are generated and remain present in the late equilibrium states. These electrostatic fields are perpendicular to the magnetic field: At |z| > 0, they points away from the equatorial plane [ Figure 6 Here the E x component is much weaker than the E z component. The electrostatic fields arise due to the decoupling of the unmagnetized ions and magnetized electrons (Schindler et al. 2012), which can be derived from the ordering among ion thermal gyroradius, current sheet thickness, and electron thermal gyroradius Taking Run 2A for example, we have ρ i : λ : ρ e = 0.32 : 1 : 0.01. As the plasma beta is lowered, the electrostatic field decreases -due to stronger magnetization of ions, which is shown in Figures 7 and 8 below. To show the detailed force balance of these current sheet configurations, we decompose the forces in parallel and perpendicular directions with respect to local magnetic fields. The force balance in the parallel direction is rapidly established because particles move freely along field lines. The unit vectors in the two perpendicular directions are defined as e ⊥1 =ẑ ×b and e ⊥2 =b × (ẑ ×b), whereb andẑ are the unit vectors along the magnetic field and z direction, respectively. It is worth noting that e ⊥1 is roughly along the ±y direction far above/below the equator, and along the −x direction at the equator, while e ⊥2 is roughly along the +z direction all over the domain.
At the equator, the electrostatic fields, E ⊥1 (i.e., along the − , which does not give net currents. All aforementioned perpendicular drift velocities are small, when compared to ion and electron parallel flow velocities that carry currents [e.g., Figure 5].
In all cases of ion-dominated current sheets, the E × B drift of ions and electrons does not give net perpendicular currents. The total electric currents are field-aligned. In addition, the plasma pressure gradients are approximately zero. Thus, the force-free configurations of the ion-dominated current sheets are obtained.

Initially electron-dominated current sheets
When electrons initially carry entirely the (field-aligned) currents (Runs 1B, 2B, and 3B), there is no change in the macroscopic states of the system at the end of simulations [e.g., Figure 9]: The currents remain field aligned, and are carried by electrons only. Compared to the ion-dominated current sheets, the relaxation of electron-dominated current sheets does not alter the macroscopic states: neither the proportion of electron and ion currents is changed, nor do the electrostatic fields develop. In all three cases with different plasma betas, the system remains in a forcefree configuration, with no plasma pressure gradients [Figures 7(d-f), 7(j-l), 7(p-r), 8(d-f), 8(j-l), 8(p-r)]. Despite the system being identical to the initial ones from the fluid point of view, the underlying velocity distribution functions evolve substantially from the initial Maxwellians to reach the Vlasov equilibrium, shown below.

Reaching kinetic equilibrium or not?
Since the force-free equilibria are reestablished toward the end of simulations for both ion-and electron-dominated current sheets from the fluid point of view, a logical next step would be to further examine if such current sheets reach the kinetic equilibrium. Because all the currents are field aligned, we integrate the full distribution function f s (x, z, v ∥ , v ⊥ , ϕ, t) of species s to obtain the reduced distribution g s (x, v ∥ , t) at the equator as a function of x, v ∥ , and t: where v ∥ and v ⊥ are the parallel and perpendicular velocities, respectively, and ϕ is the gyrophase.
, from Run 2A as an example. The ion phase space density is transported toward velocities both above and below the mean ion flow velocity [ Figure 10(a)]. Consequently, the ion velocity distribution becomes wider and less peaked, compared to the initial Maxwellian [ Figure 10(b)]. The mean ion flow velocity is slowed more on the earthward side than the tailward side [see the inset of Figure 10  The macroscopic states of electron-dominated current sheets do not differ from their initial configurations, as seen in Figure 9. The electron reduced distribution function g e (x, v ∥ ), however, shows systematic deviations from the initial Maxwellians, as displayed in Figure 11. Such deviations are similar to those in ion-dominated current sheets. The peak of the electron velocity distributions is reduced and redistributed toward larger velocities of v ∥ > 0 and v ∥ < 0 (but without changes in mean electron flow velocities) [Figures 11(c) and 11(d)]. The deviation of the ion velocity distribution from the initial Maxwellian is small [ Figure 11(b)], but is required to reach the Vlasov equilibrium.
To quantify the convergence of the system toward kinetic equilibrium, we define the metric where ⟨·⟩ x denotes the average over spatial coordinate x. This metric is chosen in such a way that the non-Maxwellian features of g s (x, v ∥ , t) are properly accounted for. During the evolution of the system, it is expected that G s (t) experiences significant changes at the beginning of the simulation and eventually reaches a steady state, if a kinetic equilibrium is established. Figure 12 shows the absolute values of the time derivative |∂G s /∂t| for both ion-and electron-dominated current sheets with different plasma betas. All runs end up with |∂G s /∂t| ≲ 0.01, except for Run 1A (i.e., the ion-dominated current sheet with β = 0.1). For the same plasma beta, electron-dominated current sheets reach the kinetic equilibrium faster than ion-dominated current sheets. In the latter case, the system simply takes more time to redistribute the currents between unmagnetized ions and magnetized electrons. For both types (ion-and electron-dominated) of current sheets, those with higher plasma β evolve toward kinetic equilibrium faster, because the current sheets with higher plasma β can support larger transient electric fields to redistribute particles more rapidly in phase space toward the equilibrium. At present, we cannot simulate further in time the relaxation of the low-β, ion-dominated current sheet in Run 1A toward the final kinetic equilibrium, because the current sheet becomes unstable at a later time (most likely due to boundary conditions). However, the progression of the relaxation for different beta values and the behavior of Run 1A up to t · ω ci ∼ 10 suggests that it also conforms to the explanatory model presented herein. The distribution functions at that time provide a sufficiently good representation of the Vlasov equilibrium to use in future modeling.

CONCLUSION AND DISCUSSION
In summary, we demonstrate that kinetic equilibria of 2D force-free current sheets exist for different proportions of ion and electron currents in various plasma betas. When initial currents are carried purely by ions, field-aligned electron currents are developed by transient parallel electric fields, while perpendicular electrostatic fields are generated due to unmagnetized ions and magnetized electrons. When initial currents are carried exclusively by electrons, the macroscopic state of the system remains unchanged from its initial state, described by the MHD model. In both scenarios, the electron and ion distribution functions at the late equilibrium states show systematic deviations from the initial drifting Maxwellians. These deviations occur in order to satisfy the time-stationary Vlasov equation.
We should mention that the proposed equilibrium must be constructed as a solution of the stationary Vlasov-Maxwell equations, but there is an infinite number of such solutions (Grad 1961). For each system, a specific solution should be determined by the boundary conditions and via relaxation (nonstationary) process. Therefore, our results show the existence of such 2D force-free current sheets, but the obtained solution should not be considered as unique. As an example, in Appendix A, we compare a broad class of theoretical solutions Neukirch et al. 2009;Neukirch et al. 2020b) for 1D force-free current sheets of B z = 0 with PIC simulations. We find that despite the existence of an apriori equilibrium set by an analytical solution for a choice of a vector potential, the simulation converges to a different equilibrium corresponding to a different vector potential, and there is no apparent control on the final solution, that is one of an infinite number of such solutions.
The existence of a kinetic equilibrium for 2D force-free current sheets indicates that there is at least one hidden symmetry in the system, implying the existence of an additional integral of motion (i.e., an additional invariant). Our system has five dimensions N D = 5 in the phase space (2D in the coordinate space and 3D in the velocity space). The Figure 11. Phase space distributions for the electron-dominated force-free current sheet with plasma beta β = 1 in Run 2B. The format is the same as Figure 10. Although the profiles of particle drift velocities are almost the same as their initial [see the insets of (a) and (c)], the particle distribution functions evolve substantially from the initial drift Maxwellians [see panels (b) and (d)].
two known invariants are the total energy H, and the y component of the canonical momentum P y = mv y + eA/c. The degrees of freedom of the system are N f = N D − N I , where N I is the number of invariants. The particle phase space densities at equilibrium are constructed as a function of such invariants of motion. To fully describe such kinetic equilibrium, which we now know it does exist, we must have N I > N f , or equivalently, N I > N D /2. Thus, N I is at least 3 (in our case) and there is at least one hidden symmetry. Starting from particle trajectory data in the equilibrium electromagnetic fields, future research on this topic could make use of machine learning models to help find the number of invariants or even discover the analytic formula of these invariants (e.g., Liu & Tegmark 2021;Liu et al. 2022 Neukirch et al. 2009;Neukirch et al. 2020b) has been obtained for the stationary Vlasov-Maxwell equations of 1D force-free current sheets with B z = 0, i.e., for 1D force-free tangential discontinuities (also called 1D force-free Harris current sheet after Harris (1962) non-force-free solution). These solutions are quite useful, but not unique for the stationary Vlasov-Maxwell equations. Therefore, although there is almost no chance to obtain the same solutions via the relaxation method using PIC simulations (because of the infinite number of possible solutions), it may be informative and useful to demonstrate the applicability of our relaxation method by comparing its 1D results with these analytical solutions.
The magnetic field of 1D force-free Harris current sheet is The current density is The vector potential in Coulomb gauge is given by The theoretical distribution function ) reads where H s is the Hamiltonian, and p s = m s v + qs c A is the canonical momentum. The subscript s = {i, e} indexes ions and electrons. The constant parameters include the density n 0s , the inverse of temperature β s , the thermal velocity v th,s = 1/ √ β s m s , the characteristic drift velocities u xs and u ys , and the dimensionless numbers a s and b s . The distribution function consists of two components: Component 1 is the familiar term carrying current in the y direction as in the Harris current sheet; Component 2 is the new term carrying current in the x direction for the force-free Harris current sheet. There are constraints on a s , b s , u xs , u ys and v th,s to ensure both the positivity of f s and a single maximum of f s in v x and v y (to avoid possible microinstabilities in velocity space; Neukirch et al. 2009). In addition, the quasi-neutrality condition relates the ion parameters with the electron parameters. We summarize the calculation of these parameters from the macroscopic parameters B 0 and λ: 2c 2c In this study, we choose m i /m e = 100, λ = 2d i , T i = 5T e = 5 12 m i v 2 A (equivalently, β i = 0.2β e = 2.4(m i v 2 A ) −1 ), and b = 1 (ensuring f s > 0, as well as a single maximum of f s in v x and v y ). Then the other parameters are determined as u yi = u xi = 5 12 v A , u ye = −u xe = − 1 12 v A , a e = 0.50, a i = 0.76, b e = 1.0, b i = 1.2, n 0e = n 0 , and n 0i = 0.81n 0 . In the PIC simulations, for the initial conditions we initialized a two-component Maxwellian that gives the same density, flow velocity, and temperature as that in Equation (A4): with The boundary condition for both fields and particles is periodic in the x direction, whereas the boundary condition in the z direction is the same as that described in Section 3. With the magnetic field and distribution function given by Equations (A1) and (A11), respectively, we let the system evolve until an equilibrium state is reached. Figure 13 shows the macrostate of the system at the final equilibrium (t = 1000 ω −1 ci ). The system is in a force-free equilibrium as evidence by P zz ≈ constant and J ⊥ ≈ 0. Compared with the initial state, the current density (mainly J y ) is bifurcated, possibly due to the orbit class transitions (Yoon et al. 2021;Yoon et al. 2023) in the relaxation process of the current sheet (see discussion of such biffurcated current sheet models in, e.g., Camporeale & Lapenta 2005;Sitnov et al. 2006). A polarized electric field points from the equator (z = 0) to higher latitudes, which is similar to the force-free Lembege-Pellat current sheet [ Figure 6]. Figure 14 shows the evolution of the reduced distribution function g s (v x ) = ∞ −∞ ∞ −∞ dv y dv z f s (v x , v y , v z ) at four representative z locations. The ion and electron distribution functions asymptotically converge to an equilibrium state, which is different from Equation (A4) , referred to as N2009 in Figure 14). In particular, the equilibrium electron distribution has a substantially larger thermal velocity than the initial one, whereas the equilibrium ion distribution has a similar thermal velocity as the initial one. Because the solution is not unique for the stationary Vlasov-Maxwell equations of 1D force-free current sheets (Grad 1961;Channell 1976), the final equilibrium distribution function is not guaranteed to be the same as Equation (A4).  ci , colored-coded from blue to red. The grey curve in each panel labeled "N2009" represents the reduced distribution from Equation (A4) ).