8pt plus 2pt minus 2pt

*,

  *

# A new tool for two-dimensional field-reversed configuration equilibrium
study

H. J. Ma1,2,3, H. S. Xie2,3∗, B. H. Deng2,3, Y. K. Bai2,3, S. K. Cheng2,3, Y.
Li2,3, B. Chen2,3, M. Tuszewski4, H. Y. Zhao2,3, and J.Y. Liu1∗ 1 Key
Laboratory of Materials Modification by Laser, Ion and Electron Beams
(Ministry of Education), School of Physics, Dalian University of Technology,
Dalian 116024, People’s Republic of China 2 Hebei Key Laboratory of Compact
Fusion, Langfang 065001, People’s Republic of China 3 ENN Science and
Technology Development Co., Ltd., Langfang 065001, People’s Republic of China
4 ENN Consultant, Riverside, CA 92506, United States of America
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A new tool (GSEQ-FRC) for solving two-dimensional (2D) equilibrium of field-
reversed configuration (FRC) based on fixed boundary and free boundary
conditions with external coils included is developed. Benefiting from the two-
parameter modified rigid rotor (MRR) radial equilibrium model and the
numerical approaches presented by [Ma et al, Nucl. Fusion, 2021], GSEQ-FRC are
used to study the equilibrium properties of FRC quantitatively and will be
used for fast FRC equilibrium reconstruction. In GSEQ-FRC, the FRC equilibrium
can be conveniently determined by two parameters, i.e., the ratio between
thermal pressure and magnetic pressure at the seperatrix $\beta_{s}$, and the
normalized scrape of layer (SOL) width $\delta_{s}$. Examples with fixed and
free boundary conditions are given to demonstrate the capability of GSEQ-FRC
in the equilibrium calculations. This new tool is used to quantitatively study
the factors affecting the shape of the FRC separatrix, revealing how the FRC
changes from racetrack-like to ellipse-like.

††:

##  1 Introduction

Field-Reversed Configurations (FRC) is considered as one of the possible
approaches to achieve fusion energy, either as a magnetic confinement approach
or as a target plasma for magnetized target fusion [1, 2, 3, 4, 5]. The plasma
equilibrium is an essential element to understand the FRC properties [6], and
is a foundation for studying various plasma phenomena such as
magnetohydrodynamic (MHD) instabilities and plasma transport.

Understanding the interior properties of hot, dense FRC plasma is still a
challenge. In early FRC experiments, profiles information were deduced from
measurements of the excluded flux array [7] and single-chord interferometry
[8], and the shape of the separatrix can be estimated. In recent FRC
experiments, multi-point Thomson scattering [9] and multi-chord interferometry
[10] provide more detailed FRC profile data. However, non-perturbative
internal magnetic field measurement is still difficult. Therefore,
theoretical/model equilibrium calculations are particularly important in the
study of FRCs, and many equilibrium models such as Rigid-Rotor (RR) [1], two-
point equilibrium (2PE) [11], three-point equilibrium (3PE) [12], symmetric
(SYM) [13] and Modified Rigid-Rotor (MRR ) [6] are proposed.

The Grad-Shafranov (G-S) equation [14, 15] can be utilized to describe the
traditional FRC equilibrium [6]. And in cylindrical coordinate system $(r,z)$
with axial symmetry, it can be written as:

$\Delta^{*}\psi\equiv r\frac{\partial}{\partial
r}\left(\frac{1}{r}\frac{\partial\psi}{\partial
r}\right)+\frac{\partial^{2}\psi}{\partial
z^{2}}=-\mu_{0}r^{2}P^{\prime}(\psi),$ (1)

with

$\displaystyle J_{\theta}=rP^{\prime}(\psi),$ (2) $\displaystyle
B_{z}=\frac{1}{r}\frac{\partial{\psi}}{\partial{r}};\quad
B_{r}=-\frac{1}{r}\frac{\partial{\psi}}{\partial{z}},$ (3)

where, $\psi=\Psi_{p}/2\pi=\int_{0}^{r}B_{z}rdr$ is the normalized poloidal
magnetic flux, $J_{\theta}$ is toroidal current density of plasma,
$P^{\prime}$ is the derivative of pressure with respect to $\psi$ and
$\mu_{0}$ is the vacuum permeability. In this work, zero toroidal magnetic
field $B_{\theta}=0$ is assumed.

Two-dimensional (2D) numerical equilibrium calculation of FRC maybe simpler
than that for tokamak with omitted $B_{\theta}$. However, the requirement of
including the plasma in the scrape of layer (SOL) region imposes some extra
difficulties. The existence of bifurcation solutions for the FRC equilibrium
calculations was a major roadblock encountered by earlier FRC equilibrium
research, and a 2D solution can be difficult to reach numerically [16, 17].
Latter, the problem was overcome by using the total toroidal current
($I_{\theta}$) [18] in the computational region or the area of plasma
surrounded by the separatrix ($S$) [19, 20] as the global constraints in
equilibrium calculations. Based on $I_{\theta}$ and $S$, it is possible to
solve FRC equilibrium with arbitrary elongation ratios. In these traditional
equilibrium calculations [19, 20, 12, 21, 22, 23, 24], iterative algorithms
are used, the free model parameters usually are not physical parameters,
therefore, post processing is required during and after the iteration to
derive physical parameters for comparing with experimental data, or for merit
assessment. The indirect algorithms are somehow inconvenient and inefficient
in equilibrium studies, for example, in FRC equilibrium design. Sometimes more
than two parameters are iterated, while as shown in reference [6], there are
only two independent parameters needed to specify an FRC equilibrium in most
models. The coupling of the model parameters could be the reason leading to
bifurcated solutions.

In the recently proposed MRR model [6], FRC equilibrium can be determined by
two physical quantities $\beta_{s}$ and $\delta_{s}$, where $\beta_{s}$ is the
ratio of the plasma pressure to the magnetic pressure at the separatrix, and
$\delta_{s}$ is the normalized SOL width. The General Solver for EQuilibrium
(GSEQ) code [25] also utilizes the algorithm proposed in reference [6] to
solve 2D FRC equilibrium, with three steps: specifying the free parameters of
the MRR model $\beta_{s}$ and $\delta_{s}$; solving other model parameters
from 1D equilibrium; and solving 2D equilibrium. Benefiting from this
algorithm, the equilibrium calculation process almost becomes an algebraic
equation solving process with high efficiency. In addition, the FRC
characteristics are readily obtained when the two model parameters $\beta_{s}$
and $\delta_{s}$ are determined. This greatly improve the efficiency to arrive
at the desired equilibrium. The convenience and efficiency of this GSEQ code
can be very beneficial in the future for equilibrium fitting of experimental
data. One unique feature of the GSEQ-FRC tool is that two key quantities of an
FRC equilibrium, the maximum vacuum magnetic field $B_{e}$ (the detailed
meaning of $B_{e}$ is discussed in A) and the separatrix radius $R_{s}$ are
scaling factors, which are not involved in the equilibrium solving process,
and can be specified after the solution. This can significantly improve the
computational efficiency in equilibrium design of FRCs.

This paper is organized as follows. Section 2 introduces the algorithms and
procedures used in GSEQ-FRC. In section 3, the 2D equilibrium with fixed
boundary condition is presented. Examples of free boundary equilibrium is
shown in section 4. In section 5, the factors affecting the shape of the FRC
separatrix are discussed. Summary and discussion are presented in section 6.

##  2  The algorithm in GSEQ-FRC

The method for solving the free parameters of equilibrium models, and the
initial result of 2D FRC equilibrium using the GSEQ code was presented in
reference [6]. In the following, the 2D equilibrium solver is described in
section 2.1. In Section 2.2, a step-by-step method for the GSEQ-FRC tool is
presented.

###  2.1 Introduction to the G-S solver

The GSEQ code, a plasma equilibrium solver with external coils, is developed
for advanced study of plasma equilibrium. The accuracy of the result from GSEQ
code has been demonstrated in B. This code supports two distinctive ‘upper’
($r\to r_{\rm wall}$) boundary conditions: fixed boundary condition for metal
wall as flux-conserver, and free boundary condition for quartz chamber. Most
of the algorithm for solving the G-S equation are well known, therefore, only
a brief outline will be given. In the flux-conserving metal wall, $\psi$ of
the ‘upper’ boundary condition is calculated using Green’s function [24]:

$\psi_{b}(R,Z)=\iint
G\left(R,Z;R^{\prime},Z^{\prime}\right)J_{\theta}\left(R^{\prime},Z^{\prime}\right)dR^{\prime}dZ^{\prime},$
(4)

where

$G\left(R,Z;R^{\prime},Z^{\prime}\right)=\frac{\mu_{0}}{2\pi}\frac{\sqrt{RR^{\prime}}}{k}\left[\left(2-k^{2}\right)K(k)-2E(k)\right]$
(5)

is the free space Green’s function, which gives the poloidal flux at ($R$,$Z$)
from a unit coil current source at ($R^{\prime}$,$Z^{\prime}$), $K(k)$ and
$E(k)$ are the elliptic integrals of the first and the second kind
respectively, and

$k^{2}\equiv\frac{4RR^{\prime}}{\left(R+R^{\prime}\right)^{2}+\left(Z-Z^{\prime}\right)^{2}}.$
(6)

In the free boundary condition case, a detailed procedure can be found in
references [26, 25]. The ‘upper’ boundary varies due to the changes of the
plasma during the numerical iterations, and the G-S equation is solved inside
the plasma-vacuum boundary. The $\psi_{b}$ on the edge of a computational
domain is calculated as a combination of two parts:

$\psi_{b}=\psi_{b}^{{coil}s}+\psi_{b}^{m},$ (7)

where $\psi_{b}^{coils}$ represents the contribution from the currents in the
external toroidal coils, and $\psi_{b}^{m}$ is the contributions from the
plasma currents using the latest approximation to $J_{\theta}$ in the
calculation proceedure. Both terms in equation 7 are computed using the
Green’s function for a toroidal current source.

Next, we consider the computational procedures of the plasma subregion. The
plasma equilibrium in the GSEQ code is found by solving finite difference
approximations to the G-S equation using the successive over relaxation (SOR)
method[18, 27]. For the FRC, the MHD equilibrium maybe more difficult with
$r\to 0$ as the ‘lower’ boundary condition. The numerical procedure is
demonstrated in figure1 and comprises of three steps: In the first step, the
information of the external coils including the current, the position and
shapes, the computational domain and $I_{\theta}$ as the global constraints on
the plasma equilibrium is given as the input data, and the guessed $\psi(r,z)$
is used as the initial $\psi$ distribution. In the second step, the G-S
equation solver is run. The new $\psi$ is updated using $P(\psi_{old})$ and
the SOR method [$\psi=w\psi_{old}+(1-w)\psi$, $w\in(0,1)$ is the acceleration
parameter] with fixed boundary or free boundary conditions respectively. The
loop is based on the convergence of $\psi$ and $\psi_{old}$. In the third
step, the updating $\psi$ is output and post-processing when equation 8 is
satisfied

$d=max\lvert\psi(r,z)-\psi_{old}(r,z)\lvert<\varepsilon,$ (8)

where $\varepsilon$ is a specified small number, typically $10^{-9}$.

Figure 1: Flow chart of the GSEQ code.

###  2.2 A step-by-step method for GSEQ-FRC

In this subsection, GSEQ-FRC tool is presented, which combines the method for
solving the free parameters of MRR model [6] and the GSEQ code. The flowchart
of the GSEQ-FRC tool is illustrated in figure 2.

Figure 2: The algorithm of GSEQ-FRC tool.

As described in reference [6], the MRR-1 model is:

$P(\psi)=\frac{B_{e}^{2}}{2\mu_{0}}\beta_{s}\cdot\exp\left(-\alpha\frac{\psi}{B_{e}R_{s}^{2}}\right)\cdot\left[\sigma\Big{(}\frac{\psi}{B_{e}R_{s}^{2}}\Big{)}^{q}+1\right]^{n},$
(9)

with $q=1$, $n=2$ as default, $\alpha$ and $\sigma$ are the free parameters.
The free parameters in equation 9 are solved from the following normalized
constraint equations [6]:

$\displaystyle\delta_{s}=\left.-\frac{p}{\frac{dp}{dr}}\right|_{r=R_{s}}=\left.-\frac{p(\psi)}{\frac{dp(\psi)}{d\psi}\cdot\frac{d\psi}{dr}}\right|_{r=R_{s}},$
(10) $\displaystyle p\left(\psi_{m}\right)-1=0,$ (11)
$\displaystyle\int_{0}^{\psi_{m}}\frac{d\psi}{\sqrt{\left[1-p\left(\psi_{m}\right)\right]}}+\frac{1}{4}=0,$
(12)

where $\delta_{s}$ is used to characterize the SOL width, and $\psi_{m}$ is
the trapped flux. The dimensionless forms are used for the radius $r/R_{s}\to
r$, the flux $\psi/(B_{e}R_{s}^{2})\to\psi$, and the pressure $P/P_{m}\to p$,
where $P_{m}=B_{e}^{2}/2\mu_{0}$. The detailed of equations 10-12 has been
shown in reference [6]. Using equations 10-12, the FRC equilibrium is uniquely
determined from $\beta_{s}$ and $\delta_{s}$.

In GSEQ-FRC tool, the important quantities $B_{e}$, $R_{s}$ and $R_{w}$ (the
chamber radius) are not involved in the equilibrium solving process, which is
because the $I_{1D}=2\pi\int_{0}^{R_{w}}j_{\theta}rdr$ as a global constraint
instead of $I_{\theta}$ in the 2D equilibrium solving of GSEQ-FRC. The
following steps are used: Specifying two parameters $\beta_{s}$ and
$\delta_{s}$ to obtain the normalized 1D equilibrium profiles. With the target
values of $B_{e}$, $R_{s}$ and $R_{w}$, the poloidal flux at the chamber wall
$\psi_{bm}$ is determined, which is used to calculate the external coil
current of the fixed boundary condition. Meanwhile, $I_{1D}$ is also
calculated to replace $I_{\theta}$ in GSEQ module. The GSEQ module is run and
the 2D equilibrium is uniquely determined when the $\psi$ value converges.
Subsequently, the 2D equilibrium properties are obtained.

##  3  2D equilibrium with fixed boundary condition

In this section, we apply the GSEQ-FRC to study some typical FRC equilibria.

###  3.1 2D equilibrium with fixed boundary condition

A simple FRC machine with three coils is shown in figure 3. The central
confinement region has 0.28 m inner diameter and 1.5 m length. The settings of
the external coils are presented in Table 1. In this subsection, the wall is
assumed to be ideal conductor (here we mean that the flux $\psi$ at wall is
taken as the vacuum flux $\psi^{coil}$ from the outside coils), so the
equilibrium can be solved with the procedure for fixed boundary case as
discussed above.

The target values of $B_{e}$, $R_{s}$ and $R_{w}$ are not involved in the
equilibrium solving when using GSEQ-FRC for the design of the 2D equilibrium.
This new tool eliminates the requirement for numerous calculations to match
with the design values. In the process of GSEQ-FRC solving 2D equilibrium, the
$I_{1D}$ is used as the global constraint instead of $I_{\theta}$, which is
the first time to link 1D equilibrium characteristics with 2D. In the
following, $B_{e}=1T$ and $R_{s}=7.5cm$ are taken as the target values.

Figure 3: The schematic of the FRC machine, where the gray boxes are the
external coils, and the solid red line is the wall. The upper half shows the
magnetic fluxes with the FRC included, and the lower half shows the magnetic
fluxes of vacuum field. Table 1: The external coils parameters.

  * Coil No. | Z (m) | R (m) | Width (m) | Hight (m)
---|---|---|---|---
1 | -0.7 | 0.15 | 0.1 | 0.006
2 | 0.0 | 0.15 | 1.2 | 0.006
3 | 0.7 | 0.15 | 0.1 | 0.006

The 2D equilibrium of FRC using fixed boundary condition is demonstrated in
figure 4. The contour of the 2D equilibrium is shown in figure 4 (a). The
axial distribution of $\psi$ and magnetic field $B_{z}$ at the wall are shown
in figures 4 (b) and (c), respectively, with external currents of $6.45\times
10^{5}$ (A) for each coils. The 1D equilibrium is the input data in the
process of solving 2D equilibrium using GSEQ-FRC tool, so the comparison of
pressure and current density profiles from solve_MRR module (1D) and the
midplane profiles (2D) of 2D equilibrium are shown in figures 4 (d) and (e)
respectively as a confidence check, which is satisfactorily.

Figure 4: The 2D fixed boundary condition equilibrium with the $B_{e}=1T$,
$R_{s}=0.075m$, $\beta_{s}=0.6$ and $\delta_{s}=0.4\delta_{s}^{RR}$, where
$\delta_{s}^{RR}$ is the $\delta_{s}$ of RR model. (a) is the contour of
magnetic fluxes of 2D equilibrium from GSEQ-FRC. (b), (c) are the profiles of
$\psi(z)$ and $B_{z}(z)$ respectively, where the red solid line is the profile
in vacuum and the green dash line is the profile of final convergent results.
(d), (e) are the equilibrium pressure and current density profiles of 1D and
the 2D in midplane respectively. The external coil currents are all
$6.45\times 10^{5}$ (A) for each coils.

###  3.2 Equilibrium with complex fixed boundary condition

In this subsection, the complex equilibrium with irregular boundary and
asymmetric coil currents is reconstructed based on the FIX device [28], which
is shown in figure 5. The equilibrium of the two mergeing FRCs is demonstrated
in figure 5 (a). It should be noticed that the currents of the external coils
are asymmetric, which can be seen from $\psi$ and $B_{z}$ profile at the walls
of figures 5 (b) and (c), respectively. The complex equilibrium is solved
using the $I_{1D}$ from MRR_module as a constraint with the total current at
the O-point of left-side FRC. The pressure and current density profiles of 1D
equilibrium are consistent with the poloidal profiles of 2D equilibrium at the
O-point, which is presented in figures 5 (d) and (e).

Figure 5: The complexed 2D equilibrium with non-regular boundary and
asymmetric coil currents. (a) is the contour of magnetic fluxes of 2D
equilibrium, (b), (c) is the profiles of $\psi$ and $B_{z}$ with $z$, (d), (e)
is the pressure and current density profiles of 1D and O-point of left-side
FRC respectively.

##  4  2D equilibrium with free boundary condition

In this section, the 2D equilibrium with free boundary condition is calculated
using GSEQ-FRC. As shown in figure 6, the 2D free boundary condition
equilibrium is presented with the same settings as those in the fixed boundary
condition case. Comparing figure 6 (a) with figure 4 (a), it is found that the
2D equilibrium of free boundary and fixed boundary of FRC are significantly
different, which is reasonable because that the boundary fluxes are very
different as can be seen from figure 6 (b) and figure 4 (b). In figure 6 (c),
it is found that $B_{z}$ at the midplane is not the preset $B_{z}=1T$, which
is explained in A. Also, the pressure and current density profiles of 1D
equilibrium match well with the midplane profiles of 2D equilibrium, which is
illustrated in figures 6 (d) and (e).

Figure 6: The 2D free boundary condition equilibrium with the setting of
$B_{e}=1T$, $R_{s}=0.075m$, $\beta_{s}=0.6$ and
$\delta_{s}=0.4\delta_{s}^{RR}$. (a) is the contour of magnetic fluxes of 2D
equilibrium, (b), (c) is the final convergent profiles of $\psi$ and $B_{z}$
with $z$, (d), (e) is the pressure and current density profiles of 1D and
midplane of 2D equilibrium respectively. Figure 7: The $B_{z}$ radial profiles
of fixed (magenta line) and free (black line) boundary condition equilibrium
near the left boundary ($z=-0.72m$).

The $B_{z}$ radial profiles near the left boundary ($z=-0.72m$) from fixed and
free boundary are shown in figure 7. It is found from figure 7 that $B_{z}$
from free boundary equilibrium is larger, i.e., the magnetic pressure is
greater, for which the plasma is more compressed axially, so that the
elongation ratio of FRC from free boundary condition equilibrium is smaller
than that from fixed boundary equilibrium.

We have found [6] that four independent parameters
$(\beta_{s},\delta_{s},B_{e},R_{s})$ can capture most of the FRC 1D
equilibrium profile features. Whereas, five independent parameters
$(\beta_{s},\delta_{s},B_{e},R_{s},l_{s})$ can probably roughly describe a 2D
FRC equilibrium. In the present work, we can control the first four parameters
$(\beta_{s},\delta_{s},B_{e},R_{s})$ well. However, how to control the FRC
length $l_{s}$ is still an open question, especially in the free boundary
condition case.

For the case of free boundary condition, the method used in the present
version of GSEQ-FRC is similar to the one that solve the free boundary
equilibrium of traditional tokamak [26]. One difference between FRC and
traditional tokamak is that the plasma exists in the region outside the
separatrix of the FRC, i.e., the coupling between FRC plasma and coil is
stronger. Due to this, the different computation positions of the left and
right boundaries will also affect the results. We also notice some previous
studies of 2D FRC equilibria [17, 27, 29, 23], which are relevant to free
boundary condition equilibria. However, they are not the same as our treatment
or the standard treatment in tokamak[26] of the free boundary equilibrium.
Therefore, the 2D free boundary condition equilibrium of FRC requires further
investigation.

##  5  The factors affecting the shape of FRC separatrix

In early 2D numerical equilibrium studies, the racetrack-like separatrix is
often obtained. However, both racetrack-like and ellipse-like separatrices
were shown in early experiments. Reference [30] investigated the factors
affecting the shape of the separatrix based on MHD equilibrium, and concluded
that steeper $P(\psi)$ can cause more ellipse-like separatrix. However,
references [31, 20] also studied this problem but obtained an opposite
conclusion, i.e., the shape of the separatrix becomes racetrack-like as
$P(\psi)$ becomes steeper at the separatrix. If we check how references [30,
31, 20] obtained their conclusions, we find those conclusions are only
qualitative. For example, the method that caused the steep $P(\psi)$ was not
mentioned clearly in reference [30]. Whereas, the parameter $\gamma$ increases
with the other parameters also be changed to achieve the ellipse-like shape in
references [31, 20] (cf. the model parameters of subplots (a) in figure 5 and
6 from reference [20]).

The factors affecting the shape of the separatrix are investigated
quantitatively in this subsection to resolve the above confusion. We use the
following equation to describe the shape of the separatrix, which is defined
in reference [32] with

$\frac{r^{2}}{a^{2}}+\frac{|z|^{m}}{b^{m}}=1,$ (13)

where $a$ is the radius of the separatrix at the midplane, $b$ is the half-
length of the separatrix, and $m$ is the shape index. The larger the parameter
$m$ is the more racetrack-like shape of separatrix. $m=2$ is the elliptical
shape in equation 13.

The fixed boundary condition equilibrium are shown in figure 8 with varying
$\delta_{s}$ but fixed $\beta_{s}$. As $\delta_{s}$ decreases, the shape of
the separatrix presented in figure 8 (a) changes from ellipse-like to
racetrack-like, and the parameter $m$ that best fits the shape of the
separatrix with equation13 increases as shown in figure 8 (b). In this
subsection, we use $-dp/dr$ to represent the pressure gradient. The pressure
gradient near the separatrix is lager with the shape of separatrix becoming
racetrack-like, which is demonstrated in figure 8 (c).

Figure 8: The changes of 2D equilibrium for different $\delta_{s}$ with
$\beta_{s}=0.6$. (a) is the shapes of the separatrix, (b) is the changes of
$m$ with $\delta_{s}$, (c) is the normalized pressure profiles. The chamber
has a length of 1.5 m and a diameter of 0.3 m, with the external coils at 0.2
m. The other values are set to $B_{e}=1T$, $R_{s}=0.075m$.

Figure 9 illustrates the changes of fixed boundary condition equilibrium for
different $\beta_{s}$ with fixed $\delta_{s}$. When $\beta_{s}$ increases, the
shape of the separatrix changes from ellipse-like to racetrack-like, and the
parameter $m$ also increases, which are presented in figure 9 (a) and (b),
respectively. Once again, as shown in figure 9 (c), the pressure gradient
becomes larger with the shape of separatrix becomes racetrack-like.

Figure 9: The changes of 2D equilibrium for different $\beta_{s}$ with
$\delta_{s}=0.15$. (a) is the shapes of the separatrix, (b) is the changes of
$m$ with $\beta_{s}$, (c) is the normalized pressure profiles. The other
values are set to $B_{e}=1T$, $R_{s}=0.075m$.

From figures 8 and 9, we find that decreasing $\delta_{s}$ or increasing
$\beta_{s}$ make the pressure gradient near the separatrix lager, and the
shape of the separatrix becomes more racetrack-like.

Next, we do an analysis of the results in references [30, 31, 20], where the
derivative of the pressure with respect to $\psi$ is

$\frac{dP(\psi)}{d\psi}=\left\\{\begin{array}[]{ll}-c(1+\epsilon\psi)&\quad(\psi\leqslant
0)\\\ -ce^{-\gamma\psi}&\quad(\psi>0)\end{array}\right.,$ (14)

where $c$, $\epsilon$ and $\gamma$ are constants. $\delta_{s}$ from equation14
is expressed as

$\delta_{s}=\frac{1}{\gamma B_{e}R_{s}^{2}\sqrt{1-\beta_{s}}}.$ (15)

From equation 15, it is found that the increase in $\gamma$ causes a decrease
in $\delta_{s}$, thus a steeper pressure near the separatrix, and the shape of
separatrix changes from ellipse-like to racetrack-like. Equation 15 reasonably
explains the conclusions in references [31, 20]. The change of $\beta_{s}$ can
also cause a change of $\delta_{s}$, which is presented in equation 15.
However, when $\delta_{s}$ is fixed, the pressure gradient increases with
$\beta_{s}$, and the shape of the separatrix becomes racetrack-like, which
would be similar to figure 9. Thus, the results in references [31, 20] are
consistent with ours. However, our results are more clearly and conclusive due
to that we can control the pressure gradient parameter explicitly and
quantitatively by using $\delta_{s}$ and $\beta_{s}$. A simple physics picture
is that racetrack-like equilibrium has larger axial gradient than the ellipse-
like equilibrium, and therefore corresponds to higher radial gradient.

##  6  Summary and conclusion

The GSEQ-FRC tool, a 2D FRC G-S equilibrium simulation tool has been developed
and is applied to FRC equilibrium design. Unlike the conventional method, in
GSEQ-FRC, the equilibrium are solved with two parameters using the physical
properties of FRC. Several examples with both fixed boundary and free boundary
conditions have demonstrated its powerful capabilities. The properties of the
FRC equilibrium with fixed and free boundaries are investigated
systematically. Furthermore, the factors affecting the shape of the FRC
separatrix are discussed quantitatively, and it is found that steeper pressure
at the separatrix causes the shape of separatrix to change from ellipse-like
to racetrack-like. The GSEQ-FRC tool has advantages in equilibrium
calculations, which can also be extended to FRC equilibrium reconstruction.
The application of the GSEQ-FRC tools to specific experiments data will be
shown in the future publications.

This work is supported by the China central government guides the development
of local science and technology funding No.206Z4501G and the compact fusion
project of the ENN group. One of the author HSX would also like to thank T.
Takahashi in Gunma University for the help in benchmarking the GSEQ code.

##  Appendix A  The accuracy of the quasi-1D equilibrium

From the equation (B.4) in reference [6], the general 1D force balance
equation from ${\bm{J}}\times{\bm{B}}=\nabla P$ is

$P(r)+\frac{B_{z}^{2}(r)}{2\mu_{0}}=\frac{B_{e}^{2}}{2\mu_{0}}+\int\frac{B_{z}}{\mu_{0}}\left(\frac{\partial
B_{r}}{\partial z}\right)dr.$ (16)

The second term on the right side is the magnetic field curvature effect at
the midplane. In equation 16, we define $B_{e}\equiv\sqrt{2\mu_{0}P_{m}}$,
where $P_{m}$ is the maximal pressure, i.e., the pressure at the O-point. The
quasi-1D equilibrium
$P+\frac{B_{z}^{2}}{2\mu_{0}}=\frac{B_{e}^{2}}{2\mu_{0}}=P_{m}$ ignoring the
curvature effect can describe the elongated cylindrical FRC well, where
$B_{e}$ can be also defined as the vacuum field $B_{z}(r\to\infty)$ or
$B_{z}(r\to wall)$ where $P\to 0$. The quasi-1D equilibrium approximation is
commonly used in both theoretical and experimental studies.

Figure 10 shows the radial profiles of the free boundary equilibrium at the
midplane in figure 6. From the red solid line and blue dashed line in figure
10, it is found that the quasi-1D equilibrium matches well inside the
separatrix, but the deviation increases outside the separatrix. However, the
left and right sides of equation 16 balanced well, which is demonstrated by
the red solid line and the brown dashed line.

Figure 10: The radial profiles of the free boundary equilibrium at the
midplane in figure 6.

Although it is found from figure 10 that the magnetic field curvature effect
is significant outside the separatrix in the 2D free boundary condition
equilibrium, the 2D results of GSEQ for $P(r)$ and $J_{\theta}(r)$ still match
well with the 1D results, which is shown in figure 6 (d) and (e). This implies
that the solution approach used in GSEQ-FRC based on MRR model is robust for
even the midplane magnetic field curvature effect is non-ignorable.

##  Appendix B  The benchmark of the GSEQ code

The GSEQ code is benchmarked with grass_ft [33]. The grass_ft is a FRC G-S
equation solver with fixed boundary condition. In grass_ft code, the ‘upper’
boundary $\psi(r_{w},z)$ is set with $r_{w}=const$:

$\psi=\left\\{\begin{array}[]{l}\psi_{w};\quad 0\leqslant z\leqslant z_{c}\\\
\frac{r_{\psi}+1}{2}\psi_{w}+\frac{1-r_{\psi}}{2}\psi_{w}\cos\left(\frac{z-z_{c}}{z_{mir}-z_{c}}\pi\right);z>z_{c}\end{array}\right..$
(17)

Here, $z_{mir}$ is the axial length from the midplane to the mirror end,
$z_{c}$ is the axial position at which the mirror field critically influences,
and $r_{\psi}$ is a control parameter for the mirror ratio. The pressure
profile with $\psi$ is

$p(\psi)=\left\\{\begin{array}[]{l}p_{s}+\frac{p_{s}}{\psi_{w}}\ln\left(\frac{p_{w}}{p_{s}}\right)\psi\\\
+\frac{1}{2}\frac{p_{s}}{\psi_{w}^{2}}\left[\ln\left(\frac{p_{w}}{p_{s}}\right)^{2}\right]\psi^{2}\quad(\psi\geq
0)\\\
p_{s}\left(\frac{p_{w}}{p_{s}}\right)^{\psi/\psi_{w}}\quad(\psi<0)\end{array}\right.,$
(18)

where $p_{s}$, $\psi_{w}$, $p_{w}$ are coefficients. In equation 18, the
normalized parameter of $p(\psi)$ is $\psi_{w}^{2}/2\mu_{0}r_{w}^{4}$. Note
that in grass_ft code, $\psi>0$ is inside the separatrix as shown in equation
18, which is opposite to the default definition in GSEQ code. We thus change
the GSEQ code to do this benchmark.

The results calculated by GSEQ and grass_ft code are shown in figure 11, where
the free parameters in equation 18 are set with $p_{s}$=5.5, $\psi_{w}$=1 and
$p_{w}=1\times 10^{-5}$. Figure 11 (a) is the specific fixed boundary set by
equation 17. Figure 11 (b) and (c) shows the contour of $\psi(r,z)$ from GSEQ
and grass_ft equilibrium codes with the same parameters respectively. Figure
11(d) illustrates that the difference of $\psi(r,z)$ between GSEQ and grass_ft
code is less than $4\times 10^{-5}$, which means that the GSEQ code agrees
well with grass_ft code.

Figure 11: The comparison of the results of GSEQ and grass_ft code. (a) is the
“upper” boundary with $r_{\psi}$=-4, $\psi_{w}=$-1, $z_{c}$=0.25 and
$z_{mir}$=0.5, (b) is the contour of $\psi(r,z)$ from GSEQ code, (c) is the
contour of $\psi(r,z)$ from grass_ft code, and (d) is the difference of $\psi$
between GSEQ and grass_ft code. The other values are set to $R_{w}=0.17m$ and
$Z_{w}=1.5m$.

## References

  * [1] Armstrong W. T. et al. 1981, Field-reversed experiments (FRX) on compact toroids. Phys. Fluids, 24, 2068
  * [2] Finn J. M. et al. 1982, Field-reversed configurations with a component of energetic particles. Nucl. Fusion, 22, 1433
  * [3] Tuszewski M. 1988, Field reversed configurations. Nucl. Fusion, 28, 2033
  * [4] Steinhauer L. C. 2011, Review of field-reversed configurations. Phys. Plasmas, 18, 070501
  * [5] Guo H. Y. et al. 2015, Achieving a long-lived high-beta plasma state by energetic beam injection. Nat. Commun., 6, 6897
  * [6] Ma H. J. et al. 2021, Two-parameter modified rigid rotor radial equilibrium model for field-reversed configurations. Nucl. Fusion, 61, 036046
  * [7] Tuszewski M. 1981, Excluded flux analysis of a field-reversed plasma. Phys. Fluids, 24, 2126
  * [8] Okada S. et al. 1989, Reduction of the density profile of a field-reversed configuration plasma from detailed interferometric measurements. Jpn. J. Appl. Phys., 65, 4625
  * [9] Deng B. H. et al. 2012, Electron density and temperature profile diagnostics for C-2 field reversed configuration plasmas. Rev. Sci. Instrum., 83, 10E339
  * [10] Deng B. H. et al. 2016, High sensitivity far infrared laser diagnostics for the C-2U advanced beam-driven field-reversed configuration plasmas. Rev. Sci. Instrum., 87, 11E125
  * [11] Steinhauer L. C. et al. 2009, Equilibrium paradigm for field-reversed configurations and application to experiments. Phys. Plasmas, 16, 072501
  * [12] Steinhauer L. et al. 2014, Two-dimensional interpreter for field-reversed configurations. Phys. Plasmas, 21, 082516
  * [13] Lee K. Y. 2020, Generalized radial profile of field-reversed configurations based on symmetrical properties. Nucl. Fusion, 60, 046010
  * [14] Grad H. et al. 1958, Hydromagnetic equilibrium and Force-Free Fields. Journal of Nuclear Energy, 7, 284
  * [15] Shafranov V. D. et al. 1966, Plasma equilibrium in a magnetic field. Reviews of Plasma Physics, 2, 130
  * [16] Marder B. et al. 1970, A vifurcation problem in E-layer equilibrium. plasma physics, 12, 435
  * [17] Spencer R. L. et al. 1982, Free boundary field-reversed configuration (FRC) equilibrium in a conducting cylinder. Phys. Fluids, 25, 1365
  * [18] Hewett D. W. et al. 1983, Two-dimensional equilibrium of field-reversed configurations in a perfectly conducting cylindrical shell. Phys. Fluids, 26, 1299
  * [19] Suzuki K. 1991, Effect of the Mirror Field on the Averaged $\beta$ Value in Field Reversed Configuration. Jpn. J. Appl. Phys., 60, 3186
  * [20] Suzuki Y. et al. 2000, Two-dimensional numerical equilibrium of field-reversed configuration in the strong mirror field. Phys. Plasmas, 7, 4062
  * [21] Steinhauer L. C. et al. 2020, Anatomy of a field-reversed configuration. Phys. Plasmas, 27, 112508
  * [22] Kako M. et al. 1983, equilibrium of field-reversed configuration with subsidiary coils. J. Phys. Soc. Jpn., 52, 3056
  * [23] Kanki T. et al. 1999, Numerical simulation of magnetic compression on a field-reversed configuration plasma. Phys. Plasmas, 6, 4672
  * [24] Gerhardt S. P. et al. 2006, Equilibrium and stability studies of oblate field-reversed configurations in the Magnetic Reconnection Experiment. Phys. Plasmas, 13, 112508
  * [25] Xie H. S. 2019, MHD Equilibrium Solver GSEQ (in Chinese). ENN
  * [26] Johnson J. L. et al. 1979, Numerical determination of axisymmetric toroidal magnetohydrodynamic equilibria. J. Comput. Phys., 32, 2
  * [27] Kako M. et al. 1983, Equilibria of Field-Reversed Configuration with Subsidiary Coils. J. Phys. Soc. Jpn., 52, 3056
  * [28] Yambe K. et al. 2009, Equilibrium of Field-Reversed Configuration Plasma Sustained by Rotating Magnetic Field. J. Plasma Fusion Res. SERIES, 8, 971
  * [29] Fuentes N. et al. 1983, ECMC, a portable two-dimensional code for plasma equilibrium computation on coaxial-multiple-coil systems. Comput. Phys. Commun., 90, 169
  * [30] Spencer R. L. et al. 1985, experimental and computational equilibrium of field-reversed configurations. Phys. Fluids, 28, 1810
  * [31] Suzuki Y. et al. 1999, Analysis of Averaged B value in Two Dimensional equilibrium of a Field-Reversed configuration with End Mirror Fields. J. Plasma Fusion Res. Series., 218
  * [32] Ohkuma Y. et al. 2010, Separatrix shape of field-reversed configuration. Phys. Plasmas, 17, 042502
  * [33] Takahashi T. et al. 2004, Losses of neutral beam injected fast ions due to adiabaticity breaking processes in a field-reversed configuration. Phys. Plasmas, 11, 3131