The equilibrium potential structure in the scrape-off layer (SOL) of the field-reversed configuration (FRC) can be affected by the penetration of edge biasing applied at the divertor ends. The primary focus of the paper is to establish a formulation that accurately captures both parallel and radial variations of the two-dimensional (2D) potential in SOL. The formulation mainly describes a quasi-neutral plasma with a logical sheath boundary. A full-f gyrokinetic ion model and a massless electron model are implemented in the GTC-X code to solve for the self-consistent equilibrium potential, given fixed radial potential profiles at the boundaries. The first essential point of this 2D model lies in its ability to couple radial and parallel dynamics stemming from resistive currents and drag force on ions. The model successfully recovers the fluid force balance and continuity equations. These collisional effects on 2D potential mainly appear through the density profile changes, modifying the potential through electron pressure gradient. This means an accurate prescription of electron density and temperature profiles is important in predicting the potential structure in the FRC SOL. The Debye sheath potential and the potential profiles applied at the boundaries can be additional factors contributing to the 2D variations in SOL. This comprehensive full-f scheme holds promise for future investigations into turbulent transport in the presence of the self-consistent 2D potential together with the non-Maxwellian distributions and open boundary conditions in the FRC SOL.

## I. INTRODUCTION

Field-reversed configuration (FRC) is an alternative approach to realizing magnetic fusion, distinguishing itself from the more prevalent designs such as tokamaks and stellarators. Characterized as an elongated prolate compact toroid (CT) with a solely poloidal magnetic field, the FRC is lauded for its engineering simplicity.^{1} Due to its easy construction in engineering aspects, FRC continues to draw research interest, especially in light of recent advancements at the TAE Technologies, Inc. The C2-W experiments at TAE successfully extend steady FRC plasmas for more than 30 ms by using a neutral beam injection (NBI) system and mirror plugs.^{2–5} Notably, the performance is reported to be mainly limited by the NBI duration. Energetic ion population from the NBI helps to stabilize the FRC plasmas because of large Larmor radius (FLR) effects^{6} and strong radial electric field shear.

The C2-W experiments also use an edge-biasing system to produce a negative radial electric field, resulting in an E × B toroidal rotation opposite to the ion diamagnetic flow. The induced rotation plays a crucial role in mitigating macroscopic wobble and rotational instabilities, enhancing the stability of the system.^{7} Additionally, the $E\xd7B$ shearing effects have been demonstrated to be beneficial for reducing the turbulent transport in FRC.^{8,9}

Microscopic drift-wave turbulence on the scale of thermal ion gyroradius is typically suppressed in the FRC core, but in the scrape-off layer (SOL), ion-to-electron scale drift-wave turbulence has been observed from the experiments and simulations.^{10–13} This observation opens up the potential for utilizing the biasing system in the SOL as a means to reduce turbulent transport and improve the FRC performance. Indeed, reduction of turbulence correlation length was reported in previous FRC measurements.^{12} In tokamaks, similar equilibrium $E\xd7B$ sheared flows have been conclusively demonstrated to suppress the drift-wave turbulence by reducing linear growth rate, nonlinear eddy size, fluctuation intensity, and turbulent transport.^{14–17} This motivates our research to investigate the effects of equilibrium $E\xd7B$ sheared flow under edge-biasing conditions on turbulent transport within the FRC SOL.

Using the gyrokinetic toroidal code (GTC),^{18} we have conducted several studies on the turbulence behavior in the FRC. Lau *et al.*^{13} demonstrated that both electron and ion-scale drift waves could become unstable in the SOL, exhibiting critical pressure gradients in line with experimentally observed thresholds. Subsequent nonlinear simulations using a global particle code ANC^{19} find that linear drift-wave instabilities first grow in the SOL, then nonlinearly spreads from the SOL to core, which exhibits a toroidal wavenumber spectrum comparable to the experimental measurements.^{20,21} A modified GTC version for the FRC geometry, GTC-X, is developed for the global simulation of nonlinear turbulent transport in the whole device.^{22} It is reported that the ion temperature gradient (ITG) mode is globally connected and axially varying across central FRC region, mirror throat area, and formation exit area. The self-generated zonal flows can suppress such ITG instability in the nonlinear simulation of the FRC SOL.^{23}

In our preliminary effort to integrate an equilibrium potential into FRC turbulence simulations using GTC-X, we discovered that combining the additional $E\xd7B$ flow with the diamagnetic flow yields a total shearing rate. This total $E\xd7B$ shearing rate can effectively suppress turbulent transport in the FRC SOL region.^{8} However, it is important to note that the equilibrium potential introduced was an unrealistic one-dimensional (1D) function of the flux surface by ignoring the presheath potentials that vary along the magnetic field-line.

When the biasing system is applied in the FRC SOL, it is important to understand how potential boundary profiles can modify the potential structure in the SOL. This understanding is crucial as the $E\xd7B$ shearing effects, introduced by the biased potential, have a direct impact on the behavior of turbulence within the region. Considering that current FRC experiments can well sustain for over 30 ms, we can reasonably anticipate the establishment of a steady-state equilibrium in the presence of a biasing system. In this stable scenario, the potential structure within the SOL can be treated as a time-invariant background, facilitating a more straightforward study of turbulence phenomena. Previous research has indicated that both the linear growth rate and the saturation timescale of turbulent transport in the FRC SOL typically occur within a 1 ms timeframe.^{8,22} This result supports the hypothesis that the background equilibrium and the turbulence phenomena can be studied independently, with the former serving as an equilibrium to the latter.

To study the effects of biasing potential on turbulence, a self-consistent model to calculate the equilibrium potential is needed. Given that particles go from the core region to the SOL, eventually exiting the device via the divertors, our sought-after equilibrium model inherently requires a two-dimensional (2D) approach, accompanied by open boundary conditions in the parallel direction. A transport model, a Quasi-1D (Q1D) code with fluid ions and electrons, together with a kinetic neutral beam species,^{24} has simulated the C2 plasma evolution with a self-consistent potential. Later, this model was extended to a 2D version, incorporating parallel transport in the SOL region. The simulation captured the density profile evolution which can be compared to the experimental observations.^{25}

To study the potential structure in SOL, a KSOL code with kinetic electrons has been developed by the TAE team to study the 1D parallel variation of presheath potential between the outer mirrors and the divertors.^{26} The simple Boltzmann response of electrons is extended by using the Vlasov–Fokker–Planck equation. This approach captures the kinetic effects on potential structure in an expanding magnetic field with anisotropic electron pressures. By comparing with the Q2D fluid model and the Boltzmann relation, it was concluded that the 1D parallel potential profile was insensitive to the different electron models, while the inclusion of ion acceleration in the model could be more important to the profile.^{25,26} This finding motivates us to consider a fluid electron model and focus on the 2D effects of poloidal flow in SOL.

In this paper, we take a different approach from the 1D KSOL simulations at TAE and focus on the 2D potential structure in the central SOL near the core region. By developing a fluid electron model in GTC-X code, our model integrates a resistive current from electron–ion collisions or other enhanced transport, as well as a drag force accounting for ion interactions with impurities and neutrals. These elements are crucial for establishing a connection between radial and parallel transport and creating a 2D potential structure in the FRC SOL. In addition, to provide the kinetic equilibrium for turbulence study, gyrokinetic (GK) ions are incorporated to obtain the kinetic effects on ions. The primary purpose is to obtain a steady-state equilibrium, ensuring its applicability as a reliable equilibrium for the upcoming turbulence simulations in the central vessel.

In the FRC SOL, electrostatic potential varies on both microscopic and macroscopic scales. On the microscopic scale, there exists a parallel potential drop over a Debye length (Debye sheath) in front of the divertors where the magnetic field-lines intercept the conducting surfaces. There is also a perpendicular potential variation over a thin layer within a width of several ion gyroradii (magnetic presheath)^{27} at the outer radius in front of the cylindrical wall in the presence of an oblique magnetic field-line. These microscopic regions are far away from the FRC core and are not simulated when we investigate the macroscopic transport of the FRC SOL in this paper. We discuss a simplified Debye sheath model at the end of the paper and highlight the need for a more accurate sheath model in our simulations. Future upgrades to the simulation boundary conditions of our hybrid gyrokinetic-fluid model could benefit from implementations found in fluid transport codes like SOLPS-ITER for tokamak SOLs.^{28,29}

The focus of the present work is on the parallel potential drop over the macroscopic region of the SOL where ions are accelerated to the ion sound speed before entering the Debye sheath. The equilibrium electrostatic potential in this macroscopic presheath^{30} depends on the magnetic structures, collisions, and atomic processes in the SOL, which affects the penetration of the divertor biasing into the confinement vessel to improve the FRC confinement.^{5,7} The gyrokinetic simulation model is valid for the quasi-neutral plasma in this macroscopic SOL region as documented in our previous turbulence simulation papers.^{8,19,22,23}

Please note that there are inherent limitations in our proposed model described in this paper, which will necessitate further improvements in future research. The model addresses only the quasi-neutral plasma before the sheath layer, employing a logical sheath^{31} boundary condition to connect our simulations with the wall potential at the divertors. Appropriate sheath models are required for accurate comparisons with experiments. Additionally, the radial coupling in the fluid equations is limited to resistive currents and a phenomenological drag force on ions. Other physics which may cause the potential structure changes due to the modification of particle distributions, such as the secondary electron emission, ionization of neutrals, and other factors, are not covered in this paper.

The rest of the paper is organized as follows. In Sec. II, we formulate the 2D equilibrium potential model with gyrokinetic ions and fluid electrons. The model describes a quasi-neutral plasma with a logical sheath boundary in FRC SOL. Both the radial and parallel dynamics are considered by including the resistive currents and drag force on ions. We verify this simulation model in the GTC-X for a simplified 1D case by neglecting the resistive and drag force in Sec. III. The ion force balance and the ion continuity equation are verified for this model. In Sec. IV, the 2D model is discussed in detail for the GTC-X simulation. Parameter scans for the resistivity and the drag force are carried out to understand how the potential structure in the FRC SOL region is established. The current density distribution in the FRC SOL is also discussed with this new model. Section V discussed the influence of the potential boundary profiles on the 2D structure of equilibrium potential. Section VI discussed the Debye sheath potential model that was required to connect our simulations to the actual wall potential at the divertor ends. The demand for a better sheath potential model was discussed in this section. Section VII is a summary of the simulation results and potential issues and future plan for this 2D FRC equilibrium model with edge biasing.

## II. FORMULATION OF 2D EQUILIBRIUM POTENTIAL MODEL IN FRC SOL

The primary objective of this 2D equilibrium potential model including the presheath region is to calculate a self-consistent potential, given a preassigned boundary. This equilibrium can subsequently be linked to the edge-biasing applied at the divertor walls by incorporating an appropriate Debye sheath model. The 2D model described here lays the foundation for future transport studies by providing a self-consistent gyrokinetic particle distribution function as an equilibrium. As illustrated in Fig. 1, the gyrokinetic (GK) ion simulation enables us to derive the ion guiding center distribution. From this, we can calculate essential parameters such as particle density, ion parallel flow velocities, and pressures. With these quantities, we employ the Hall-MHD electron model and a quasi-neutral condition to solve for the electrostatic potential, thereby advancing our understanding of the system's behavior and laying the groundwork for future explorations of turbulence phenomena.

### A. Gyrokinetic model for ion species in equilibrium simulation

The gyrokinetic formulation is applicable when dealing with physical phenomena characterized by a frequency much smaller than the cyclotron frequency and a gyroradius much smaller than the system size. In regions like the FRC core and near magnetic null points in the FRC, where the gyroradius is relatively large, a fully kinetic approach is more appropriate. However, in the SOL, the gyrokinetic assumptions remain valid due to the stronger magnetic field and significant temperature drop at the edge, which reduces the gyroradius. While high $\beta $ values (the ratio of kinetic pressure to magnetic pressure) may be expected near the FRC core, our simulations typically yield a $\beta $ of around 0.01 in the SOL. As a result, we employ the electrostatic gyrokinetic equation to calculate the ion distribution function, which can be later utilized in turbulence studies.

^{32}

^{,}

^{33}$Cfi$ represents a general collision operator as discussed in Ref. 34, and the electrostatic potential, $\varphi $, is gyro-averaged through the operator $\varphi =\u222bd\alpha 2\pi \u222bdx\varphi \delta x\u2212R\u2212\rho $.

^{35–38}In this method, the distribution function $fi$ is divided into two components: an equilibrium part, denoted as $fi0R,v\u2225,\mu $, and a time-varying perturbation part, $\delta fiR,v\u2225,\mu ,t$. The equilibrium component satisfies the zeroth-order equation

^{8,22,23,39}

However, when computing a self-consistent equilibrium, the local Maxwellian solution becomes inappropriate. This is particularly evident in the SOL region, where particles can escape the machine through open field-lines. To address this issue, we must assume open boundaries at the divertor ends, which can result in a non-Maxwellian distribution of particles with a macroscopic flow velocity. Consequently, constructing an analytical form of $fi0$ for Eq. (2) in this open boundary scenario becomes challenging.

To construct the complicated equilibrium, we now pursue a full-f simulation approach. In this method, we no longer separate the $fi0$ from $\delta fi$. Instead, we evolve the complete gyrokinetic equation with the full-f weight, $wf=fig$, where $g$ is the marker distribution. The full-f weight should satisfy $Lwf=Cfi$, given the marker distribution $g$ as also a solution for the Lagrangian operator. In this case, the full-f weight remains constant throughout the simulation if we ignore the effects of collision on $g$.^{40}

^{33}To obtain the particle density, $Ni=\u222bdx\u2032dv\u2009Fi\delta (x\u2032\u2212x)$, the first term gives us the guiding center density

^{41}

When calculating the pressure, it is important to consider that the pressure is defined by the averaged kinetic energy in the reference frame moving with the average flow velocity of each species. Particularly in the FRC SOL, a notable parallel flow is directed toward the divertors. Therefore, when computing pressure in Eq. (6), it is essential to subtract the fluid flow velocity $Ui0$. The contribution of perpendicular flow velocity, being 1–2 orders of magnitude smaller than the parallel flow, is omitted in the formulation. This observation is confirmed through simulations discussed in subsequent sections.

To compute polarization quantities, obtaining the gyro-averaged potential $\varphi $ is crucial. Presently, the GTC-X code performs gyro-averaging for a specific toroidal mode number only,^{8} so we here use a long-wavelength approximation for the gyro-average.^{41} Future developments may focus on implementing a more generalized gyro-averaging algorithm within the GTC-X framework.

Another challenge in analytically integrating the polarization quantities in Eqs. (3)–(6) is to get the guiding center density distribution $fi$. Since the actual density distribution can deviate from the local Maxwellian, we adopt a shifted Maxwellian as a simplified model, $fi\u2248fSM=ni0mi2\pi Ti03/2exp\u2212miv\u2225\u2212Ui02+2\mu B2Ti0$. This heuristic model offers valuable insights into the approximation of polarization terms in Eqs. (3)–(6) in the long-wavelength limit. Evaluating integrals with the shifted Maxwellian is straightforward, which yields the guiding center contributions: $n\u0131\xaf=ni0$, $ui\u2225\xaf\u2248U\u01310$, and $p\u0131\xaf=ni0Ti0$.

To compute the polarization terms analytically, we first compute the gyro-averaged potential as $\varphi R=\u2211k\varphi keik\u22c5RJ0k\u22a52\rho i2$, where $\rho i=1\Omega iTi0mi=cmiTi0ZiB$ and $J0$ is the Bessel function.^{41} Employing the long-wavelength approximation ( $k\u22a5\rho i\u226a1$), we can express $\varphi \u0303$ using the Padé approximation, $\varphi \u0303x=11\u2212\rho i2\u2207\u22a52\varphi x$. Subsequently, the polarization density can be approximated as $ni,pol\u2248Zini0Ti0\rho i2\u2207\u22a52\varphi $. Intriguingly, in the long-wavelength limit, the particle parallel flow velocity is the same as the guiding center flow velocity, $Ui\u2225=Ui0$, as long as we correctly account for the polarization density. Following a similar approach, we can obtain $pi\u2225,pol\u2248ni,polTi0$ and $pi\u22a5,pol\u224832ni,polTi0$. The particle density and pressure, incorporating polarization corrections, now become $Ni=ni01+ZiTi0\rho i2\u2207\u22a52\varphi $, $Pi\u2225=ni0Ti01+ZiTi0\rho i2\u2207\u22a52\varphi $, and $Pi\u22a5=ni0Ti01+32ZiTi0\rho i2\u2207\u22a52\varphi $.^{42} Although a heuristic shifted Maxwellian distribution is used in our derivation, these expressions serve as a useful lowest-order approximation for polarization effects without loss of generality.

### B. Fluid model for electron species in equilibrium simulation

For classical transport induced by electron–ion collisions, the distortion of the electron distribution function results in different parallel and perpendicular resistivity with $\eta \u2225=0.51\eta \u22a5$ for an ion charge $Zi=e$.^{43} However, this ratio might differ when considering anomalous transport, such as turbulence-induced transport and other resistive effects. When scanning a wide range of resistivity values in the later sections of this paper, we use the same parallel and perpendicular resistivity ( $\eta =\eta \u2225=\eta \u22a5$) for simplicity.

In the current work, only isotropic pressure, $Pe=Pe\u2225\u2248Pe\u22a5$, is used in the fluid equations. When accounting for temperature variations in different direction, it becomes necessary to substitute the pressure gradient term $\u2207Pe$ with the divergence $\u2207\u22c5Pe$, where $Pe=Pe\u2225\u2212Pe\u22a5bb+Pe\u22a5I$ is the pressure tensor with distinct parallel and perpendicular pressure components. This modification will introduce an additional mirror force term related to the magnetic field changes, which will be discussed in Sec. IV. For derivations in this section and Appendix A, we will assume an isotropic pressure in all the equations.

To calculate the potential and address the other unknown variables in Eq. (7), we introduce simplifications in our model. We treat the magnetic field as static throughout the simulation, assuming that the magnetic field's evolution and dissipation occur over a longer timescale. We also neglect the change of the magnetic field due to change of the plasmas current, i.e., we assume that the divertor biasing does not significantly change the magnetic field. In this way, ions and electrons reach equilibrium with the electrostatic potential before any significant magnetic field alterations take place. Thus, our simulation results can be interpreted as an equilibrium state given a specific magnetic field structure.

To calculate the electron pressure, we consider the electron temperature $Te$, with the equation of state $Pe=NeTe$. Because electron parallel transit time is much shorter than the collisional time, we assume that electron is isothermal, i.e., electron temperature is a function of the flux surface only, $Te=Te\psi $. In the current simulations, we use a predefined radial profile for electron temperature to avoid the need for the electron energy equation. Neglecting the parallel electron temperature variation in our model results in the exclusion of the thermal force in Eq. (7) from the electron–ion collision, which potentially alter the electrostatic potential.

### C. Implementation of the fluid equations in field-line coordinates for GTC-X code

In GTC-X, we use a field-aligned mesh for the field solver.^{22} The coordinate system $\psi ,\zeta ,S$ includes the poloidal magnetic flux $\psi $, the angular coordinate $\zeta $, which is consistent with the cylindrical coordinates, and the normalized field-line distance coordinate $S$, which relates to the distance along each magnetic field-line. $S\u22080,1$ can be treated as a substitute for the parallel coordinate $Z\u2208\u2212Z0,+Z0$ when we constructed the GTC-X code.

^{22}Here, we have already simplified the equations by assuming symmetry in the toroidal direction, i.e., $\u2202\u2202\zeta =0$.

To numerically solve Eqs. (12)–(17), we convert the covariant components to the contravariant components, using the geometric tensor $g\alpha \beta =e\alpha \u22c5e\beta $, such as $J\alpha =g\alpha \beta J\beta $. This ensures that all the equations involve only six unknown contravariant variables. For simplicity, we have not included the detailed forms of the convection term, $Ui\u22c5\u2207Ui$, and the diffusion term, $\nu \u22072Ui$, in the field-line coordinates. The detailed computation of the geometric tensor and the convection/diffusion terms will be presented in Appendix A. Throughout the simulation in this paper, we have not included the diffusion term $\nu \u22072Ui$ due to its complex form in the field-line coordinates, which may be implemented in the future work. Nevertheless, the effects of the convection term $Ui\u22c5\u2207Ui$ are discussed in detail in Sec. IV.

To improve clarity, we have organized our equations in such a way that a specific variable is listed on the left-hand side of each equation. This indicates which variable that each equation addresses. Equation (12) reveals that the radial current in our model arises from flow convection, drag force, and diffusion. Due to the presence of radial current, the quasi-neutral condition, Eq. (13), necessitates a parallel current. The parallel Ohm's law then determines the 2D potential structure based on the resistive parallel current and pressure gradient. Eqs. (15)–(17) give the value of $J\zeta $, $Ui\psi $, and $Ui\zeta $ that these flow velocities again serve as important sources to generate the radial current. These six equations form a closed equation set that can be simultaneously solved using the HYPRE solver implemented in the GTC-X code.^{44}

## III. VERIFICATION OF EQUILIBRIUM POTENTIAL MODEL WITH NO RESISTIVITY AND DRAG FORCE

^{45}

The potential boundary $\varphi 0(\psi )$ here is right at the edge of the quasi-neutral plasmas, implying that we need to add the additional potential drop across the Debye sheath layer when considering the external potential set by the biasing system on the divertor walls. The true potential on the wall $\varphi w$ should be the sum of both the boundary potential $\varphi 0$ in this paper and the Debye sheath potential drop $\varphi D$, $\varphi w\psi =\varphi 0\psi +\varphi D\psi $. Incorporating this Debye sheath potential with a proper sheath boundary model is needed for calculating the structures of the parallel current and flow near the divertor.

In our simulation, the potential drop across the sheath, $\varphi D$, can be either calculated using the logical sheath boundary^{31} or estimated through simple Debye sheath models. For example, $\varphi D$ can be modeled as a function of electron temperature, $\varphi D\psi =\u2212Te\psi elnmi4\pi me$. This formulation assumes zero current through the Debye sheath, cold ions moving at a speed of $2Temi$, and electrons with a half-Maxwellian distribution flowing out of the region.^{46} In our present model, with both $\varphi w$ and $Te$ as preassigned values, $\varphi D\psi $ will be a fixed offset to the wall potential $\varphi w$ on each flux surface. Consequently, we can obtain a fixed boundary condition, $\varphi 0=\varphi w+Teelnmi4\pi me$ for the quasi-neutral plasma.

In future simulations, if the electron temperature is modeled self-consistently or a more comprehensive Debye sheath model is employed, $\varphi D$ may depends on variables such as the current density $J$. In such cases, the boundary condition for the quasi-neutral plasma $\varphi 0$ must align with these simulated quantities rather than being a fixed value as currently used. In Sec. VI, we will discuss the selection of appropriate Debye sheath models for our simulations and will examine the implications of assuming zero current through the Debye sheath. Before that, $\varphi 0$ will remain a fixed boundary condition, based on the simplified Debye sheath model, to facilitate discussions on the effects of other parameters.

To establish an equilibrium, it is necessary to include a particle source in the simulation, compensating for particles leaving the FRC SOL region in the parallel direction. A straightforward approach is to refill the central region with the same number of particles lost at the boundaries, which conserves the total number of particles in simulations. A practical implementation of this concept could involve a local Maxwellian velocity distribution $fM$, with similar profiles as the initial density $ni0\psi $ and temperature $Ti0\psi $. Such model maintains similar thermal properties as the initial state, which can be regarded as an external particle injection with zero flow velocity and fixed density and temperature. In this paper, we use such a simple source model for verifying the 2D equilibrium model in the FRC SOL. In the future study, we will use more realistic particle and energy sources from experimental measurements or modeling that include the particle transport from the core to the SOL and external injections such as the neutral beam injection and the particle fueling.

Following these considerations, the density source rate is chosen as $n\u0307src\psi ,S\Delta t=AsrcSni0\psi $. $AsrcS$ is a normalization factor that ensures the total refilled particle number $Nsrc$ equals to the lost particle number $Nloss$ within a time interval $\Delta t$, i.e., $Nsrc=\u222bdRAsrcSni0\psi =Nloss$. In the simulation, $AsrcS$ is held constant within the region $Z/R0\u2208\u22124.0,\u20094.0$ and set to zero outside the region. This will result in a uniform density source rate within the region, as the plateau shown in Fig. 3(b). It is worth noting that this is a simplified model, and more realistic source density profiles may be explored in future studies, particularly as more experimental data becomes available.

When refilling the lost particles back to the center with a local Maxwellian, part of them will be trapped by the mirror throat of FRC. When the system reaches equilibrium, the particle source rate is balanced with the collisional loss rate toward the divertors. To simulate the portion of ions constantly escaping the magnetic potential well, we use an ion–ion pitch-angle scattering operator to randomize the ion velocities, corresponding to $Cfi$ in Eq. (1). In this way, some of ion particles can be scattered into the loss cone and escape the region in the parallel direction. This scattering operator is crucial in this 1D simulation to establish such dynamic balance with the source $n\u0307src$.

Finally, the ion flow boundary we used was consistent with the so-called logical sheath boundary^{31} to satisfy the quasi-neutral condition. A simplified estimation of the Debye sheath potential drop is to assume a certain form of the electron distribution function and compute the flux as $\Gamma evce=\u222bvce\u221edv\u2225\u222b\u2212\u221e\u221edv\u22a5\u2009Fev$, where $vce$ is the cutoff velocity for the slowest electron. Since the electron flux is the same as the ion flux, $\Gamma e=\Gamma i$, we can obtain the critical velocity of reflected electron $vce$ from the ion flux and then estimate the sheath potential drop as $\varphi D=\u2212m2evce2$.^{31} This sheath potential should be added to the potential boundary $\varphi 0$ in our model to calculate the potential drop between the divertor wall and the plasma. Future research will focus on simulating the actual electron distribution to refine the sheath potential drop calculation in order to compare with experiments.

To verify Eqs. (20) and (21), the GTC-X simulation uses the 2D axisymmetric FRC magnetic geometry provided by the MHD equilibrium code LR_eqMI,^{47} with the FRC magnetic field structure described in the original GTC-X formulation paper.^{22} The poloidal flux surfaces and magnetic field geometry within the simulation domain in this paper is plotted in Fig. 2. The radial simulation domain spans from $\psi R=1.43,Z=0$ to $\psi R=2.40,Z=0$, with $R$, $Z$ value normalized to the radius of the equilibrium magnetic axis $R0=26.8\u2009cm$.

The parallel domain encompasses $Z=\xb120.0$. Notably, the inner boundary approaches the separatrix $\psi sepR=1.423,Z=0$, and the parallel boundary includes the inner mirror throat location around $Z=\xb19.37$. For the C2-W device, a second mirror throat is located near $Z=\xb130.0$. The 1D KSOL simulations by TAE team focus specifically on the second mirror location, where the magnetic field-line expands rapidly after the mirror throat, and a kinetic electron model is used to explore the deviation from the simple Boltzmann electron response.^{26} Due to the drastic density drop after the second mirror throat, including this region in GTC-X simulations up to the divertor would require substantially more computational resources. For this reason, the area near the second magnetic mirror location was not included in our current simulations. However, the formulation presented in this paper remains applicable for possible extensions into the magnetic expansion region toward the divertors.

The ion temperature is a function of $\psi $, exhibiting a sharp decrease from $861$ to 8 eV toward the outer radius, and the electron temperature is a constant $Te=80\u2009eV$ in this equilibrium. Although $Ti$ near the separatrix is approximately ten times that of $Te$, $Ti$ quickly drops to a comparatively low value, resulting in a mean ion temperature of $206\u2009eV$. The density profile also undergoes a significant drop from $0.857n0$ to 0.081 $n0$, where $n0=2.44\xd71013\u2009cm\u22123$ represents the peak density at the magnetic axis.

In our simulations, we use an explicit time integration with a 2nd-order Runge–Kutta method to evolve the dynamical system including the gyrokinetic and fluid equations. The time step, $4.3\xd710\u22125\u2009ms$, is chosen based on the convergence study in our previous paper^{8} to accurately resolve ion guiding center motion, which sets the shortest timescale of the simulated system. The simulation time for the system to evolve toward an equilibrium solution is on the order of 1 ms, which is much longer than that ion pitch-angle scattering time.

Figures 3(b) and 3(c) present the parallel variation of each term in Eqs. (20) and (21), evaluated at a mid-radius flux surface where the initial $Ti$ is comparable to $Te$. The amplitude of magnetic field, in Fig. 3(a), is normalized to $B0=530.65\u2009G$ at the magnetic axis. A pronounced peak in magnetic field strength is observable at the mirror throat locations $Z=\xb19.37$. The large gradient leads to significant variation of the $\u2207\u2225B$ terms in these equations.

As we can see in Fig. 3, the correct calculation of these mirror force effects can be important when achieving the force balance and the particle conservation. The source term $n\u0307src$ manifests in Fig. 3(b) as a plateau within $Z=\xb14$, marking the region where particles lost in the parallel direction are uniformly replenished.

Figure 4(a) shows the 2D potential structure, and Fig. 5(a) shows the density distribution in the simulation. Due to the particle loss at $Z=\xb120$, a low-density region can be observed outside the mirror throat. Despite a density peak near the inner boundary around $Z=0$, the highest potential values are found near the outer boundary. This is because we choose a boundary condition of $\varphi W=\varphi D$ ( $\varphi 0=0$) at $Z=\xb120$. In this case, the potential change only depends on the relative density drop $Ni\psi ,SNi\psi ,1$ along each field-line, instead of the absolute value $Ni\psi ,S$. Larger potential value at the outer boundary only indicates a larger relative ratio of density on those flux surfaces.

By verifying Eqs. (20) and (21) through our simulation, GTC-X's ability to accurately represent the parallel physics of ion species is confirmed. The next step is to extend the electron model, specifically Eq. (18), to couple the radial and parallel physics for a more comprehensive understanding of electrostatic potential structure in SOL. This will be addressed through an exploration of Eqs. (12)–(17) in Secs. IV–VI.

## IV. EFFECTS OF RESISTIVITY AND DRAG FORCE IN THE 2D FRC SOL MODEL

The drag force can originate from various collisional processes or turbulent fluctuations. We use the ion–ion collisional frequency, $\nu ii=4.8\xd710\u22128niln\Lambda Ti3/2mi/mp$, as a reference for this drag force in our paper. With a typical ion temperature of $Ti=200\u2009eV$ and $mi/mp=2$, we have $\nu ii=4.45\xd7103\u2009s\u22121$. The corresponding value for the drag force coefficient is $\sigma 0=miNi\nu ii$.

As we mentioned in Sec. II B, the fluid model Eqs. (12)–(17) were derived based on the assumption of isotropic pressure. However, when taking into account the anisotropy in pressure, a mirror force term will appear in the fluid equation. This term is important when we explain the force balance shown in Fig. 3(c). Following the algebra shown in Appendix B, the divergence of pressure tensor will give us an additional term associated with $\u2207\u22c5bb$. When the parallel ( $P\u2225$) and perpendicular ( $P\u22a5$) pressures are different, we need to modify the pressure gradients, $\u2202P\u2202S$ by $\u2202P\u2225\u2202S+P\u22a5\u2212P\u22251B\u2202B\u2202S$, and $\u2202P\u2202\psi $ by $\u2202P\u22a5\u2202\psi +P\u2225\u2212P\u22a5\u2207\xd7b\zeta RB$ for both electrons and ions. In this paper, as $Te$ is assumed to be a constant, the electron pressure is always isotropic, $Pe\u2225=Pe\u22a5$. The only modification under this assumption is in Eq. (17), the ion radial force balance. Furthermore, modifications would be necessary to accommodate anisotropic electron pressure in future studies.

In the 1D simulation, the parallel equation is decoupled from other equations and only the potential $\varphi 0$ at the plasma edges on each flux surface is needed to solve the electrostatic potential. However, more boundary conditions are needed in the 2D simulation with resistivity and drag forces. In the parallel direction, we still provide a fixed $\varphi 0\psi $ at the plasma edges. The parallel current at the midpoint is assumed to be zero, i.e., $J\u2225=J\u22c5b=0$ at $Z=0$ due to the symmetry between positive and negative Z regions. In the radial direction, we apply additional radial boundary conditions for potential and current by assuming, $\u2202\u2202\psi 2=0$, i.e., $\varphi $ and $J\psi $ changes linearly at the radial boundaries when calculating the radial derivatives in Eqs. (13), (15), and (17).

After determining all the necessary parameters in the fluid electron model, we conduct a simulation with $\eta =\eta 0$ and $\sigma =\sigma 0$. Similar to the 1D case, several ion collisional times are simulated to obtain the steady-state equilibrium. Figures 4(b) and 5(b) display the resulting electrostatic potential and the density structure. Comparatively, the density profile from this simulation appears less peaked in the center region than what was observed in the 1D parallel simulation, which can be attributed to the influences of resistive and collisional effects. Especially, the drift velocity $vF$ due to the drag force diffuse ion particles away from the core region. This phenomenon is more distinctly visible in the potential plot of Fig. 4(b), where more contour lines extend toward the outer radial boundary and the parallel divertor boundary, beyond $Z>9.37$. These results underscore the significance of incorporating 2D effects into the simulation to capture the complex behavior of the plasma, especially in the presence of resistive and collisional interactions.

Figure 6(a) illustrates the current density flow within the FRC SOL, revealing a pronounced current density near the axis. This enhancement at smaller radii is partly attributed to the FRC's cylindrical geometry, where the volume increases with radius, leading to a lower current density at a larger radius. The current density is presented in GTC units, $JGTC=en0R0\Omega cp=5.33\xd7106\u2009A/m2$, where $\Omega cp=5.08\xd7106\u2009s\u22121$ is the gyrofrequency of a proton at the magnetic axis $B0$. Maximum amplitude for poloidal current is $2.4\xd710\u22123$, while the toroidal current ranges from $\u22123.1\xd710\u22122$ to $1.3\xd710\u22122$.

The inclusion of the convection term, represented by $Ui\u22c5\u2207Ui\u22c5e\zeta $ in Eq. (12), can play a crucial role in predicting the current distribution in the SOL. This term can be a main source for radial current $J\psi $, when ion flow velocity is significant. Including this term ensures the model accounts for the influence of ion flow inertia, particularly in regions beyond the mirror throat where the parallel flow can reach values comparable to the sound speed. Similar level of parallel ion flow was also reported in previous Q2D simulations.^{25}

When the convection term is omitted in the simulation, the results show a significant change in the current density distribution, even though the potential and ion density remain similar to the previous case with the convection term included. In Fig. 6(b) (without the convection term for $\eta =\eta 0$ and $\sigma =\sigma 0$), the current direction is reversed near the axis when the convection term is not included. Maximum amplitude for poloidal current in this case is $1.3\xd710\u22124$, while the toroidal current ranges from $\u22123.0\xd710\u22122$ to $8.9\xd710\u22123$.

The convection term contributes an additional radial current $J\psi =miNiBS\omega i\psi UiS$, as derived in Appendix A. Due to a relatively large ion toroidal flow (shown by the contours in Fig. 6) and parallel flow, a stronger radial current is generated and modifies the current structure in FRC SOL. Comparing Figs. 6(a) and 6(b), the poloidal current structure can be correlated with the toroidal flow when the convective effects are considered in the simulation. Therefore, it can be essential to include the convection term when predicting the current distribution.

In GTC-X, the convection term is approximated by using the flow velocity from previous time step. Although such treatment can be used to achieve a steady-state solution, a more rigorous nonlinear solver for these terms should be developed in the future works. The inertial forces effects associated with strong perpendicular flows are neglected in the current gyrokinetic model in the GTC-X. Systematic derivation of such effects exists for tokamaks with strong toroidal rotation.^{48–50} These effects will be incorporated in our future work for a more consistent simulation of the convection effects in the FRC.

Next, to investigate how the resistivity and drag force can modify the 2D potential structure, we vary the resistivity $\eta $ and the drag force $\sigma $ by applying a multiplicative factor. In Fig. 7, the radial profile of the potential $\varphi $ at $Z=0$ is presented for changing drag force coefficient $\sigma $ from $1\sigma 0$ to $50\sigma 0$, while fixing the resistivity at $\eta 0$. The dashed line serves as a reference line from the 1D simulation in Sec. III.

With an increase in $\sigma $, ions are subjected to stronger drag force, leading to their redistribution toward the outer radius and the parallel boundaries. This redistribution moderates the density gradient observed in the 1D simulation, resulting in a less pronounced potential drop in the parallel direction. When we fix the potential as zero (end-shorting boundary) at $Z=\xb120$, lower potential is observed at the center of the simulation domain, as shown in Fig. 7(a).

Since the $E\xd7B$ shearing rate is more relevant in the turbulence suppression, we also plot the shearing rate, $\omega s=R2Bd2\varphi d\psi 2$, in Fig. 7(b). When the resistivity and drag force coefficient are set to their baseline values, $\eta =\eta 0$ and $\sigma =\sigma 0$, the inclusion of the drag force and the resistivity smooth out the potential profile which slightly reduces the shearing rate. However, as $\sigma $ continues to increase, the impact of these collisional factors on the shearing rate becomes more and more obvious, due to particle density redistribution. This means, under strong collisional effects, an accurate estimation of particle density in SOL can be important for predicting the $E\xd7B$ shearing rate.

In Fig. 8, we carry out a similar parameter scan for the resistivity $\eta $ from $1\eta 0$ to $50\eta 0$, while fixing the drag force coefficient at $\sigma 0$. The overall potential drop between the 1D simulation and the $\eta =\eta 0$, $\sigma =\sigma 0$ case is mainly due to the inclusion of a drag force $\sigma =\sigma 0$ as we have discussed in Fig. 7. For the case “^{*} $\eta =50\eta 0,\u2009\sigma =\sigma 0$,” we do not include the convection term in Eqs. (12) and (17) due to the numerical stability when including convection term with such a large resistivity.

The inclusion of $\eta $ near the classical resistivity value has relatively minor effects on the $E\xd7B$ shearing rate, as compared to the inclusion of $\sigma $. As we will see in the following discussion, this is attributed to the fact that the resistive current term $\eta J$ plays a minimal role in balancing the parallel forces. The predominant factor influencing the potential change is the alternation of the density profile, which is mainly affected by the drag force in the current parameter range.

Based on the parameter scan results, we choose a resistivity value that is five times the classical value, $\eta =5\eta 0$, for subsequent discussions. Past research has indicated that such a resistivity value is more pertinent for accurately forecasting the C2 FRC plasma's experimental evolution.^{25} The drag force coefficient is equivalent to the ion–ion collisional timescale, $\sigma =\sigma 0$, allowing us to observe the impact of this added drag force within reasonable simulation timeframe. In future studies, we anticipate adopting more realistic resistivity and drag force coefficients once we obtain additional experimental measurements.

To better understand the electron fluid model, we plot out the electron force balance in Fig. 9 with $\eta =5\eta 0$, $\sigma =\sigma 0$. Along the parallel direction at the middle flux surface $R=1.89$ in Fig. 9(c), the electrostatic potential, labeled as “ $\u2212\u2202S\varphi $,” is mainly balanced by the electron pressure, labeled as “ $\u2202SPe$.” This observation is similar to the 1D case, where the Boltzmann-like response, Eq. (19), can be effectively employed to estimate $\varphi $ based on the ion density. The resistive current effects “ $\u2212\eta JS$,” on the other hand, can be negligible when $\eta $ is around classical value.

The electron radial force balance at $Z=0$ in Fig. 9(a) reveals significant contributions from $Ui\xd7B$ and $J\xd7B$, which together give the electron flow $Ue\xd7B$. This toroidal electron flow is then comparable to the $\u2207Pe$ and the $\u2207\varphi $ term in Eq. (7). In the toroidal force balance at $R=1.89$ in Fig. 9(b), the absence of toroidal variations due to the symmetry condition $\u2202\u2202\zeta =0$ results in a balance between the radial current $J\psi $, the ion radial flow $U\psi $ and the resistive current $J\zeta $. There is an interesting separation of two regions for the toroidal force balance: within the mirror throat ( $Z=9.37$), the radial ion flow and the toroidal resistivity are the main contributions; outside the mirror throat, the radial ion flow and the radial resistive current dominates the toroidal force balance. This could be attributed to the increase in ion flow toward the divertors, which in turn generates a higher plasma current as the particles exit the SOL.

## V. PENETRATION OF PREASSIGNED BOUNDARY POTENTIAL PROFILE ALONG MAGNETIC FIELD-LINES

After the detailed analysis of the 2D potential model, the focus now shifts to understanding how the edge biasing voltage impacts the potential structure in the confining vessel of the FRC SOL. In experiments, the divertor biasing is achieved by an array of concentric annular electrodes to achieve discrete voltages at different radial locations, and these electrodes are insulated from each other.^{4,5,7} To achieve smooth potential structures in simulations, we manually select a continuous radial profile $\varphi w$ instead of discrete potential steps to qualitatively resemble the potential drops created by electrodes.

Assuming the Debye sheath drop $\varphi D$ as a function of $Te$ and the current $J$, we have the boundary condition $\varphi 0=\varphi w\u2212\varphi DTe,J$. In parctice, the wall potential can be controlled by the electrodes, allowing $\varphi 0$ to be determined self-consistently based on temperature and current data. For the simple Debye sheath we discussed, $\varphi 0=\varphi w+Teelnmi4\pi me$ can be a fixed potential boundary in the simulations. However, this model assumes zero current through the Debye sheath. In Sec. VI, we will explore the validity of this assumption and discuss under what conditions it may hold.

In this section, we explore how changes to the simplified boundary potential affect the 2D equilibrium potential. Notice that the Debye sheath potential $\varphi D$ in this form is no different than providing an additional shift of $\varphi 0$ profile. We can produce the desired potential boundary profiles by selecting appropriate wall potentials. However, one should always keep in mind that this relationship may not hold if $\varphi D$ also depends on the current, in which case its effects on $\varphi 0$ would be more complex.

To implement the fixed potential boundary, we first use an analytic profile $\varphi 0\psi =\varphi 0,max21\u2212tanh0.34\u2212\psi \xaf\u2212\psi \xaf00.29$ as the potential boundary condition at the quasi-neutral plasma edge $Z=\xb120$, where $\psi \xaf=\psi \psi 1\u2212\psi 0$ is the poloidal flux normalized by the simulation domain size of $\psi 0,\psi 1$, and $\varphi 0,max$ is the maximum potential of the radial profile. For simplification, we describe our simulations by $\varphi 0$ in most cases since it is directly applied at the boundaries of our simulation domain. A discussion on the difference between $\varphi 0$ and $\varphi w$ is provided in Sec. VI.

When there is no resistivity and drag force ( $\eta =0$, $\sigma =0$), the potential at the boundary is simply a constant shift of potential value along each flux surface. Same amount of $E\xd7B$ shear will be added in the parallel direction, saying that the boundary $\varphi 0$ can be accurately transfer to the center region of FRC. When the resistivity and drag force are included, there can be 2D variations due to the radial coupling, so the introduced $E\xd7B$ shear from the potential boundary may change when penetrating toward the middle region of FRC.

Figure 11 illustrates how the radial potential profile at the plasma center of $Z=0$ can change when we apply different potentials at the plasma edge of $Z=\xb120$. In general, the potential profile we applied at the boundary introduces a potential shift throughout FRC SOL. With $\eta =5\eta 0$, $\sigma =\sigma 0$, the resistive $\eta J$ term has little effects in the parallel electron force balance, as we discussed in Sec. IV. Thus, a solution similar to Eq. (19) can be obtained, meaning that the total potential can be separated into two parts: the boundary condition $\varphi 0\psi $ and the Boltzmann-like response due to $\u2207Pe$.

In Fig. 11(a), the potential difference between the quasi-neutral plasma edge ( $Z=20$) and the center profile ( $Z=0$) already exists when a zero potential boundary $e\varphi 0,max=0$ is used. A density structure with $e\varphi 0,max=200\u2009eV$ is also provided in Fig. 12(b), which can be qualitatively compared to Fig. 5(b). The similarity in density structures, whether or not a non-zero potential profile is applied, leads to a similar potential drop along the parallel direction in both scenarios.

However, if we take a closer look at the $E\xd7B$ shearing rate in Fig. 12(b), the situation is more complex than a straightforward superposition of these effects. Specifically, the $E\xd7B$ shearing at $Z=0$ with $e\varphi 0,max=200\u2009eV$, denoted as $\omega SZ=0,\u2009200\u2009eV$, cannot be simply expressed as the sum of the shearing due to the potential $\varphi 0\psi $, i.e., $\omega SZ=20,\u2009200\u2009eV$, and the shearing caused by original density variation with zero edge profile, $\omega SZ=0,\u20090\u2009eV$. The density structure is also modified when an edge potential profile is introduced into the FRC SOL, and more subtle changes exist in the radial potential profiles.

To see a stronger effects of the resistive current, we also run a non-zero $\varphi 0$ case with $\eta =100\eta 0$, $\sigma =10\sigma 0$ in Fig. 13. Due to stronger effects of the resistivity and the drag force, the radial profile of the electrostatic potential at $Z=0$ can differ more from the applied $\varphi 0$ boundary. The toroidal electron force balance and the drift velocity in ion dynamics, $vF$, can both modify the density distribution in FRC SOL and finally changes the radial potential profile.

## VI. DEBYE SHEATH MODEL FOR A BOUNDARY CONDITION

To conclude the discussion, we want to show that Debye sheath potential drop can possibly affect the 2D equilibrium potential structure. As discussed in Secs. III–V, the boundary condition $\varphi 0\psi $ is the potential profile at the edge of the quasi-neutral plasma, so an additional Debye sheath potential drop should be included if we want to compare the simulation with the actual potential value $\varphi w$ at the divertor walls. As defined in Sec. III, $\varphi 0=\varphi w\u2212\varphi D$, the zero boundary conditions $\varphi 0=0$ in Secs. III–V correspond to $\varphi w=\varphi D$. That is, an externally applied potential on the divertor wall is set to be identical to the Debye sheath potential drop.

In Secs. III–V, a simple Debye sheath model was assumed, $\varphi D\psi =\u2212Te\psi elnmi4\pi me$, based on a textbook derivation.^{46} Given that the electron temperature is predetermined in our model, the required additional potential drop $\varphi D$ is also a fixed profile throughout the simulation. Essentially, this approach acts as an offset correction to the profile $\varphi W$ and leads to a fixed boundary $\varphi 0$ at the quasi-neutral plasma.

For comparison between the simulations with $\varphi w=\varphi D$ ( $\varphi 0=0$) and with $\varphi w=0$ ( $\varphi 0=\u2212\varphi D$), we use an electron temperature profile $Te\psi =Te01.0+0.25tanh0.34\u2212\psi \xaf\u2212\psi \xaf00.29\u22121$, with $Te0=80\u2009eV$. Figure 14(b) illustrates the potential structure under the $\varphi w=0$ conditions. The potential $\varphi 0$ at $Z/R0=\xb120$ is larger at the inner radius, due to the higher electron temperature from the specified profile. This leads to a higher Debye sheath potential drop $\varphi D$.

When analyzing the differences in the 2D potential structure with $\varphi w=\varphi D$ and $\varphi w=0$, the variations can be elucidated in a manner akin to what is depicted in Fig. 12(a), but now with a $\varphi 0$ profile decreasing in the radial direction. Figure 15(b) illustrates the radial force balance incorporating the sheath potential drop as a boundary condition. Here, an extra $E\xd7B$ flow can be observed if we look at the radial profile labeled as “ $\u2212\u2202\psi \varphi $” in Figs. 15(a) and 15(b). This flow emerges predominantly due to the Debye sheath potential boundary, where $\u2202\varphi 0\u2202\psi =\u2212\u2202\varphi D\u2202\psi $.

The Debye sheath potential drop in our simulation, as a simple function of electron temperature, essentially acts as an additional adjustment to the $\varphi 0$ profile at the boundary. This model assumes zero net current through the Debye sheath, $J\u2225=0$, by balancing the ion and electron contributions. Specifically, at the boundary of the plasma, the parallel current $J\u2225$ can be expressed as the sum of the contributions from ions and electrons, $J\u2225=Zi\Gamma i\u2225\u2212e\Gamma e\u2225$, where $\Gamma i\u2225$ and $\Gamma e\u2225$ are parallel particle fluxes of ions and electrons, respectively.

Given our maintenance of the quasi-neutral condition, the total ion flow from all boundaries must equal the total electron flow. This means that any local current entering the region at one boundary must exit at another. Illustratively, if we examine the radial profile of the parallel current at $Z=20$ shown in Fig. 6(a), we observe current entering the region from the larger radius and exiting at the smaller radius, as further detailed in Fig. 16(a).

Figure 16(b) presents the current flow when a boundary potential profile with $e\varphi 0,max=200\u2009eV$ was applied. Although a lower potential at the smaller radius tends to direct currents toward that area, this trend was not pronounced in our simulations. The current pattern remained largely unchanged after altering the potential boundary profile, consistent with our earlier findings that most potential variations at the boundary primarily influence $E\xd7B$ flow. Additionally, the effects of resistive current terms are relatively minor in the overall force balance.

This observation suggests that assuming zero current through the Debye sheath and applying a simplified Deba sheath model at the larger radius may be a reasonable approximation. However, this should not be considered a definitive conclusion, as our fluid model includes several simplifications which limits our ability to compare with actual experiments. In the future, more comprehensive Debye sheath models will be necessary to accurately match the current and potential drop near the boundary.

For example, considering a Debye sheath model that includes non-zero current through the Debye sheath boundary. This necessitates accounting for the current flow balance $J\u2225=Zi\Gamma i\u2225\u2212e\Gamma e\u2225Te,\varphi D$ at the boundary. Assuming electrons flow out of the boundary with a half-Maxwellian distribution, we obtain an electron flux, $\Gamma e\u2225Te,\varphi D=ne0Te2\pi meexpe\varphi DTe$, from which the Debye sheath potential can be estimated with the simulated $J\u2225$ and $\Gamma i\u2225$. However, this still requires several assumptions, such as a definite flow direction of electrons, i.e., $J\u2225<Zi\Gamma i\u2225$. Also, special attention is meeded when $\Gamma e\u2225$ approaches zero, since $\varphi D\u2192\u2212\u221e$ given finite $ne0$ and $Te$. These potential issues underline the need for future tests and implementations of a self-consistent Debye sheath boundary in the model.

Future investigations will utilize current flow data at the boundary to explore more sophisticated Debye sheath models and validate them against experimental data, improving the model's fidelity. Potential useful references for developing such sheath boundaries could be the Chapter 3 in a recent textbook by Rozhansky, 2023^{51} and the related treatment in the SOLPS-ITER code.^{28,29} The simple Debye sheath model used in this study mainly acts as a preliminary tool, highlighting the importance of incorporating a precise sheath model for more accurate simulations.

## VII. CONCLUSION

In this paper, an axisymmetric 2D equilibrium potential in the FRC SOL is simulated by assuming the timescale separation between the turbulence evolution and the steady-state equilibrium with a preassigned potential boundaries. The formulation to calculate the presheath equilibrium includes a full-f gyrokinetic ion model and a massless electron model for a quasi-neutral plasma before entering the microscopic Debye sheath layer or thin magnetic presheath layer. Due to the particle flow through the SOL toward the divertors, the equilibrium presheath potential is intrinsically 2D, which involves the balance between radial and parallel transport.

The model was first verified in a simplified simulation when the resistivity and the drag force on ions are neglected, which reduces the equilibrium to 1D. The simplified model uses a Boltzmann-like electron response and can restore the ion parallel force balance and continuity equation on each flux surface. The potential boundary in this 1D equilibrium can perfectly transfer to the center region of the FRC SOL from the divertors.

After verifying the 1D physics in the parallel direction, parameter scans were performed to analyze the influences of resistivity and drag force on the 2D equilibrium. When considering classical resistivity corresponding to electron–ion collisions, the resistive current exhibited minimal effects in electron parallel force balance. The drag force, initially approximated using the ion–ion collision frequency, induced an outward transport of particles directed toward the divertors. This results in a slightly lower potential profile in the center region of the SOL. A more pronounced impact on electron toroidal force balance and subsequent radial transport behavior was observed when the resistivity was increased to two orders of magnitude above the classical value, combined with a drag force equivalent to ten times the ion–ion collision frequency.

The parallel electron force balance is always Boltzmann-like within the parameter ranges in this paper, which provides a simple solution to estimate the potential structure in the FRC SOL. The resistive current can be important in explaining the electron toroidal force balance, but not in the electron parallel and radial force balance. The collisional effects mainly appear through the density profile change, and then modify the potential structure through electron pressure gradient.

The effects of a changing radial profile of electrostatic potential at the quasi-neutral plasma boundary were investigated. However, an additional potential drop, associated with the Debye sheath layer, should be added for accurate comparisons with the actual potential at the divertor walls. In this paper, we only explored a simple Debye sheath model by including a Debye sheath potential drop proportional to electron temperature. This simplified Debye sheath model makes the Debye sheath potential drop an additional correction to the potential boundary applied at the edge of quasi-neutral plasma. This can be a rough estimation since the present simulation settings predicted a relatively small current flow. However, more accurate Debye sheath model should be considered to self-consistently match the sheath drop and the current flow, in order to compare with experimental measurements in the future.

Eventually, the interplay between the potential edge profiles and the parallel variation of potential due to the density structure within the SOL will determine the radial electric field at the center of the FRC SOL. The $E\xd7B$ shearing rates can also be drastically changed under different resistivity and drag force. An accurate estimation of density profiles can be a key point in order to correctly predict the electrostatic equilibrium potential.

Although we have scanned a certain range of parameter space of our 2D model, simulations with a more realistic setting are still required. In this paper, a uniform density sourcing rate within a fixed simulation region is used, but the real particle source should be particle transported from the core region into the SOL region, together with other particle injections, such as NBI systems and fueling system. More accurate resistivity and drag force coefficients are also needed, which requires more experimental data such as the neutrals density. Additionally, current information at the boundary should also be considered for a self-consistent calculation of Debye sheath potential drop as we just mentioned. A well-defined Debye sheath model could significantly enhance the capability to accurately estimate the potential changes near the divertor walls.

Though not discussed in this paper, our full-f scheme has successfully replicated the ion temperature gradient (ITG) instability observed in previous delta-f simulations. This comprehensive full-f scheme holds promise for future investigations into turbulent transport in FRC, incorporating the dynamic 2D equilibrium obtained in this paper.

## ACKNOWLEDGMENTS

The authors would like to thank Dr. Laura Galeotti for providing the FRC equilibrium data, Dr. Robert D. Falgout for setting up the HYPRE linear solver in the GTC-X code, and Dr. Peter Yushmanov for helpful discussions. This work was supported by TAE Grant No. TAE-200441, DOE SciDAC, and INCITE programs and (T.T.) U.S. DoE ECP (Exascale Computing Project). Simulations used resources on the Oak Ridge Leadership Computing Facility at Oak Ridge National Laboratory (DOE Contract No. DEAC05-00OR22725) and the National Energy Research Scientific Computing Center (DOE Contract No. DE-AC02-05CH11231).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**W. H. Wang:** Formal analysis (equal); Writing – original draft (equal). **X. S. Wei:** Validation (equal); Writing – review & editing (supporting). **Z. Lin:** Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). **C. Lau:** Visualization (supporting); Writing – review & editing (supporting). **S. Dettrick:** Visualization (supporting); Writing – review & editing (supporting). **T. Tajima:** Resources (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

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

### APPENDIX A: CALCULATION OF GEOMETRIC TENSOR AND EQUATIONS IN GENERAL COORDINATES

^{22}Using transformation functions in cylindrical coordinates, we can calculate the values of geometric tensor in the physical space

*et al.*

^{22}

Note that the parallel components in the $b$ direction differ from the projection onto the field-line basis $eS$. In the FRC geometry, the magnetic field without toroidal component is represented as $B=\u2207\psi \xd7\u2207\zeta =BRR\u0302+BZZ\u0302$, with $eS=g\u2207\psi \xd7\u2207\zeta $. This leads to the relationship, $b=BB=1BgeS$. Therefore, the parallel flow can be related to the covariant component as $Ui,\u2225=Ui\u22c5b=1BgUS$ in the FRC geometry. By expressing $UiS=gS\psi Ui\psi +gSSUiS$, we can substitute $UiS=gSS\u22121UiS\u2212gS\psi Ui\psi $ into the equation. The equation set can be further simplified by recognizing $B=Bb=1geS$. Thus, $BS=1g$ and $BS=B2g$.

Up to now, we have successfully derived a set of linear expressions for our equations, except for the nonlinear convection term, $Ui\u22c5\u2207Ui$. To address this challenge, one approach is to utilize the flow velocity from the previous time step and calculate the complicated convection term as a known source term in the current time step. Because we are interested in the steady-state solution, this approach can be reasonable, provided that the simulation can numerically converge. In the future, we may explore the development of a nonlinear algorithm to rigorously solve these nonlinear terms.

### APPENDIX B: CALCULATION OF PRESSURE TENSOR IN FRC GEOMETRY

^{22}so we can write

We simply need to replace $\u2202P\u2202S$ by $\u2202P\u2225\u2202S+P\u22a5\u2212P\u22251B\u2202B\u2202S$, and $\u2202P\u2202\psi $ by $\u2202P\u22a5\u2202\psi +P\u2225\u2212P\u22a5\u2207\xd7b\zeta RB$ for any species in Eqs. (12)–(17), to account for the additional effects from $\u2207\u22c5bb$ term in the pressure tensor.

## REFERENCES

*Electrostatic Turbulence and Transport in the Field-Reversed Configuration*

*Introduction to Plasma Physics and Controlled Fusion*

*Computational Plasma Physics: With Applications to Fusion and Astrophysics*

*f*particle-simulation method

*Plasma Theory: An Advanced Guide for Graduate Students*