Reconstruction of Coronal Magnetic Fields Using a Poloidal-Toroidal Representation

A new method for reconstruction of coronal magnetic fields as force-free fields (FFFs) is presented. Our method employs poloidal and toroidal functions to describe divergence-free magnetic fields. This magnetic field representation naturally enables us to implement the boundary conditions at the photospheric boundary, i.e., the normal magnetic field and the normal current density there, in a straightforward manner. At the upper boundary of the corona, a source-surface condition can be employed, which accommodates magnetic flux imbalance at the bottom boundary. Although our iteration algorithm is inspired by extant variational methods, it is non-variational and requires far fewer iteration steps than most of them. The computational code based on our new method is tested against the analytical FFF solutions by Titov&D\'{e}moulin (1999). It is found to excel in reproducing a tightly wound flux rope, a bald patch and quasi-separatrix layers with a hyperbolic flux tube.


INTRODUCTION
To understand physical processes in the solar corona and their connection with dynamics of the heliosphere, information about the coronal magnetic field is indispensable. Measurements of coronal magnetic fields have been made by use of radio observations (Alissandrakis & Chiuderi Drago 1995;White et al. 1991;Brosius et al. 1997;Lee et al. 1998;White 2005), by spectropolarimetry of near infrared lines (Lin et al. 2004;Tomczyk et al. 2008;Judge et al. 2013;Plowman 2014;Dima & Schad 2020) or by utilizing a magnetic-induced transition line in Fe X and other Fe X and Fe XI lines in EUV (Landi et al. 2020). These techniques as of today, however, can only provide a coarse 2D map of the lineof-sight magnetic field, the magnitude of the magnetic field or a pointwise vector field. Coronal magnetic field strengths have also been indirectly estimated from observations of coronal loop oscillations (Nakariakov & Ofman 2001;Van Doorsselaere et al. 2008) and of coronal Corresponding author: G. S. Choe gchoe@khu.ac.kr mass ejections (CMEs) (Jang et al. 2009;Gopalswamy et al. 2012). These indirect measurements are too locally confined to give us a geometrical picture of the coronal magnetic field. On the other hand, the magnetic field vectors in the photosphere have been measured by spectropolarimetry with a sufficiently high resolution and accuracy to produce a vector magnetogram, which provides a 2D cut near the solar surface through the 3D vector field (Beckers 1971;Harvey 1985;Solanki 1993;Lites 2000).
With photospheric vector magnetograms available, there have been efforts to reconstruct coronal magnetic fields. Although the real corona is fully dynamic, the magnetic field that can be generated by one vector magnetogram without knowing its prehistory is a static field. A magnetohydrostatic field under gravity requires the information of pressure (and temperature if not isothermal) at the photospheric level in addition to the vector magnetogram (Grad & Rubin 1958;Wiegelmann et al. 2007;Zhu & Wiegelmann 2018), but such information is not readily available yet. Since the plasma β, the ratio of plasma pressure and magnetic pressure, in the corona, especially in active regions, is much less than unity (Iwai et al. 2014, cf. Rodríguez Gómez et al. 2019, the approximation of coronal magnetic field to a force-free field (FFF) has been prevalent. An FFF in a domain V is a vector field B such that As in equation (1), we will omit constant coefficients in Maxwell's equations by a proper normalization throughout this paper. From equations (1) and (2), we have i.e., α is constant along each field line. If the scalar field α is constant everywhere in V , equations (1) and (2) form a linear vector Helmholtz equation and its solution, an FFF with a constant α, is a linear force-free field (Aly 1992). If α at the boundary ∂V is non-constant, the scalar field α(r) for r ∈ V is an unknown, and such an FFF is a nonlinear force-free field (NLFFF). Since the scalar α in the photosphere is far from constant, reconstruction of a coronal magnetic field amounts to seeking an NLFFF. The earliest attempts of constructing an NLFFF used an algorithm, in which the three components of the vector field are computed from the photospheric boundary successively upward plane by plane (Nakagawa 1974;Wu et al. 1985Wu et al. , 1990) as if a time-dependent hyperbolic partial differential equation (e.g., an advection equation) were solved marching along the time coordinate, whose role is taken by the vertical coordinate z in the FFF solver. As already pointed out by Grad & Rubin (1958), a magnetohydrostatic (MHS) equation has imaginary characteristics like an elliptic equation (e.g., Poisson equation) as well as real characteristics like a hyperbolic equation. Thus, solving a force-free equation, which is semi-elliptic, as a Cauchy problem is ill-posed so that the successive integration algorithm cannot avoid exponentially growing errors as going up to a higher altitude (Aly 1989;Amari et al. 1997;McClymont 1997). For this reason, the algorithm is now seldom used, but it has left a legacy of the much-used term "extrapolation." In this paper, we will use the term "reconstruction" instead of the somewhat misleading "extrapolation" to refer to solving for a coronal magnetic field with certain boundary conditions. Among a variety of coronal NLFFF solvers that have been proposed so far (for a review see Aly & Amari 2007;Wiegelmann & Sakurai 2021), there are two major groups of methods that are in practical these days. They are Grad-Rubin methods (current-field iteration methods) and variational methods. The former were originally proposed by Grad & Rubin (1958) and have been applied in diverse formulations and algorithms (Sakurai 1981;Amari et al. 1999Amari et al. , 2006Wheatland 2006). The Grad-Rubin methods in common employ an iteration procedure, in which the domain is first loaded with a field-aligned electric current J m+1 = α m B m satisfying equation (3) for B m at the iteration step m and α = J ·n/B ·n specified at the boundary points (footpoints) having one sign of the normal field component B ·n in ∂V , and then B m+1 is updated by solving the equation ∇×B m+1 = J m+1 . Although the Grad-Rubin problem, which is to solve equations (1)-(2) with B n given at every point of ∂V and α only in the part of ∂V with one sign of B n , is known to be well-posed for |α| < α m < ∞ (Bineau 1972;Boulmezaoud & Amari 2000), it has not yet been rigorously proved whether the Grad-Rubin iteration procedure always converges to a solution or not. However, numerical codes based on the Grad-Rubin method have successfully been applied to real solar problems demonstrating its usefulness (Bleybel et al. 2002;Régnier & Amari 2004;Régnier & Priest 2007;Petrie et al. 2011;Mastrano et al. 2020).
In the variational methods, we solve for the magnetic field that extremizes (actually minimizes) a certain functional, which is usually a volume integral involving magnetic field, with certain boundary conditions and some additional constraints if any. For example, if a certain field line connectivity is imposed in V and the conjugate footpoints of each field line are fixed in ∂V , and if the magnetic field in V is varied under the ideal magnetohydrodynamic (MHD) condition without footpoint motions in ∂V so that the first two constraints may be maintained, the magnetic field that minimizes the functional W = 1 2 V B · B dV , the total magnetic energy in the domain V , is a force-free field (Grad & Rubin 1958;Chodura & Schlüter 1981). A sufficient condition for δW ≤ 0 is that the fictitious plasma velocity is proportional to the Lorentz force, i.e., v ∝ J × B. This physically implies that the force-free state can be approached by removing the kinetic energy, into which the excessive potential energy is converted, from the system possibly by a hypothetical friction and/or viscous diffusion. However, the problem of reconstructing coronal magnetic field is quite different from that treated by Chodura & Schlüter (1981). Since we do not know the field connectivity beforehand, we need to lift the ideal MHD condition and the constraint of field connectivity intentionally. Instead we have to impose the normal current density J n or three components of B at ∂V . Then, an energy-decreasing evolution tends to deplete magnetic helicity through its dissipation within the system and transport through the boundary (Berger 1984). Thus, maintaining J n or B at ∂V requires winding up field line footpoints there. Magnetofrictional methods for coronal FFF reconstruction either alternate stressing and relaxing steps explicitly (Mikić & McClymont 1994;Roumeliotis 1996;Jiao et al. 1997) or are inherently equipped with a rather automatic re-stressing mechanism (Valori et al. 2005(Valori et al. , 2010Inoue et al. 2014;Guo et al. 2016;Jiang & Feng 2016). Since the magnetofrictional codes are more or less modified forms of MHD solvers, they can employ free boundary conditions, also called open boundary conditions, at the outer boundary (lateral and top boundaries for a box-shaped domain) (e.g., Valori et al. 2010). This flexibility in the outer boundary conditions helps the system to evolve toward a force-free state in most cases, but it is uncertain what mathematical problem the resulting force-free state is the solution of.
Another group of most widely used variational methods are the so-called optimization methods, in which the functional is to be minimized (Wheatland et al. 2000). The simplest choice of the boundary condition well-posing this variational problem is ∂B/∂t = 0, which can straightforwardly be implemented. There has been a concern that a boundary condition giving all three components of magnetic field results in an overspecified problem in contrast to the Grad-Rubin formulation (Grad & Rubin 1958;Boulmezaoud & Amari 2000). However, available observational data of photospheric magnetic field and the boundary conditions employed at the outer boundary (lateral and top boundaries for a box-shaped domain) can hardly satisfy the compatibility relations of force-free fields (Aly 1989;Wiegelmann & Sakurai 2021). The mathematically clear formulation of the optimization methods tells that specifying all three components of magnetic field at the entire boundary works for minimizing the functional L in equation (5) to a value, which is not necessarily zero, even if the compatibility conditions are not met. This robustness is a great merit of the optimization method as is its well-posedness as a variational problem. It is, however, a shortcoming of the method that the outcome strongly depends on what B is given at the outer boundary (Wiegelmann et al. 2006a). It should be also noted that in the optimization method, the nonzero ∇ · B term arises not merely from numerical discretization errors, but rather from the nondivergence-free form of ∂B/∂t required for reducing the functional L given above 1 . As a remedy for this concern, it has been suggested to impose ∇·B = 0 as a constraint using a Lagrange multiplier in the variational procedure (Nasiri & Wiegelmann 2019). It has also been tried to purge the final solution of nonzero ∇ · B by a postprocessing (Rudenko & Dmitrienko 2020), which, however, entails an unignorable change in either the normal or tangential component of magnetic field at the bottom boundary.
Besides the two classes of methods above, the boundary integral method, also called boundary element method, has also been applied in the solar context (He & Wang 2008;He et al. 2020;Guo et al. 2019). This method is quite similar to the Green function method for partial differential equations, and a surface integral involving a reference function over the photospheric boundary needs to be evaluated to obtain B at each coronal point at every iteration step (Yan & Sakurai 2000;Yan & Li 2006). It is a merit of the boundary integral method that a finite computational domain does not need to be set up, nor any artificial boundary conditions at the outer boundary, because it assumes a semiinfinite domain and a finite total magnetic energy in it.
All the aforementioned methods solve an FFF problem with B n and J n or B as the bottom boundary condition. The problems so posed aim to reconstruct coronal magnetic fields with vector magnetograms at hand. Our new method presented here also tackles such problems. However, FFF problems may be posed in different ways depending on one's interest. If it is necessary or desirable to impose a certain field connectivity, one can use a magnetofrictional method using Euler potentials (e.g., Choe & Cheng 2002) or the fluxon method using thin, piecewise linear flux tubes called fluxons for magnetic field description (DeForest & Kankelborg 2007). In another problem setting, an FFF solution was sought with a flux rope initially placed at a desired location (van Ballegooijen 2004;van Ballegooijen et al. 2007). Here, only B n is imposed at the bottom boundary and J n comes out of the solution.
With such a variety of methods for coronal FFF reconstruction available today, we still want to present a new method, which has the following properties.
(1) The magnetic field is described in a way ensuring divergence-freeness.
(2) The boundary conditions at the bottom boundary, B z (z = 0), and J z (z = 0), are straightforwardly implemented once and for all.
(3) The lateral and top boundary conditions can accommodate magnetic flux imbalance at the photospheric boundary.
(4) Fewer iteration steps are required for convergence than in most other methods.
(5) The numerical code is robust and equally operative for simple and complex field geometries.
Let us make some remarks on the above items. The most fundamental way of describing a divergence-free (solenoidal) vector field is using a vector potential such that B = ∇ × A. Imposing B z at the boundary z = 0 would be readily done by fixing appropriate A x and A y at the boundary once and for all. However, J z or B x and B y cannot be determined by the values of A at the boundary only. Therefore, we have to adjust some components of A at z = 0 at every iteration step in keeping with the evolution inside the domain (e.g., Roumeliotis 1996). Thus, items (1) and (2) cannot be realized together when vector potentials being used. In this paper, we describe magnetic field with two scalar functions Φ and A z , named poloidal and toroidal functions, which will be explained in detail in Section 2. Such a field description has never been tried in numerical computation of solar magnetic fields whether static (FFF included) or dynamic. In our formulation, the divergence-freeness of B is guaranteed, and B z is solely represented by Φ, and J z by A z . Thus, setting boundary conditions at z = 0 is done once and for all. With regard to item (3), we note that the total positive and negative fluxes are generally not equal in magnitude in the magnetogram of any active region. The imbalanced fluxes are often coerced into achieving a balance by preprocessing (Wiegelmann et al. 2006b), or are accommodated in FFF computations by assuming certain symmetries across the lateral boundaries (Seehafer 1978;Otto et al. 2007). In our model, we set up the lateral boundaries as rigid conducting walls so that no magnetic flux may escape the domain through the boundaries and field lines tangential to them may slide freely during iterations. As for the top boundary, we employ a source surface condition, in which the magnetic field should have only the normal component (B z ), but no tangential components (B x = B y = 0). Thus, the unpaired extra magnetic flux at the bottom boundary is connected to the top boundary so that the condition S B ·n dS = 0 should be met. The source surface boundary condition at the outer corona was first suggested by Altschuler & Newkirk (1969) and has since been widely applied to potential field models. A theory of NLFFFs with a source surface boundary condition was put forward by Aly & Seehafer (1993), but FFF reconstruction with B n and J n together or B as the bottom boundary condition and with the source surface condition at the top boundary has not been attempted before the present paper 2 .
New NLFFF solvers must be tested against known analytical FFFs before being applied to solar vector magnetograms. Currently two analytical NLFFF solutions are widely used as reference fields. The FFF models by Low & Lou (1990) (hereafter LL) are exact analytical solutions, which involve modestly sheared magnetic fields without flux ropes. The models by Titov & Démoulin (1999) (hereafter TD) are approximate analytical solutions involving a flux rope and a background magnetic field in equilibrium. Most NLFFF solvers presented so far have well reproduced the LL fields, particularly when the analytic solutions are used as boundary conditions at all six boundaries (Schrijver et al. 2006). The TD models are more difficult to reconstruct, especially in generating a single flux rope structure, but a few codes have done the job well (Valori et al. 2010;Jiang & Feng 2016). We have also tested our new code against those analytical models, focusing more on the TD models, which have much more complex field topology than the LL models.
In this paper, we present a new method of coronal FFF reconstruction and its test against analytical models. Its application to a solar active region will be given in a sequel. In Section 2, the poloidal and toroidal representation of magnetic field is expounded. Then, we define the problems to be solved and explain our numerical algorithm in Section 3. In Section 4, we present the tests of our new method for TD models in comparison with other methods. Lastly, a discussion and summary are provided in Section 5.

POLOIDAL-TOROIDAL REPRESENTATION OF MAGNETIC FIELD
Let us now consider a domain, which encloses a star, but excludes the star, and set up a coordinate system such that the stellar surface is a coordinate surface of one coordinate named ξ. For example, if the domain is the space exterior to a spherical star of radius R, the stellar boundary is the surface r = R, i.e., ξ = r, and the natural choice of the coordinates would be spherical coordinates (r, θ, ϕ). If the domain is a semi-infinite space above a plane, the planar stellar boundary is the surface z = 0, i.e., ξ = z, and the natural choice of the coordinates would be Cartesian coordinates (x, y, z) or cylindrical coordinates (ρ, ϕ, z). A magnetic field (a solenoidal vector field) in such a domain can in general be decomposed into a poloidal field B P and a toroidal field B T (Elsasser 1946;Lüst & Schlüter 1954;Chandrasekhar & Kendall 1957;Chandrasekhar 1961;Backus 1958;Rädler 1974;Stern 1976;Backus 1986;Low 2006Low , 2015Berger & Hornig 2018;Yi & Choe 2022), i.e., in which Here ∇ξ =r in spherical coordinates and ∇ξ =ẑ in Cartesian or cylindrical coordinates. The scalar fields Φ and Ψ, respectively, are called poloidal and toroidal scalar functions (Backus 1986) or Chandrasekhar-Kendall functions (Montgomery et al. 1978;Low 2006). We will simply name them the poloidal and toroidal functions in this paper, and we will call the magnetic field description given by equations (6)-(8) the poloidaltoroidal representation (hereafter PT representation) 3 . As can be seen in equations (7) and (8), B P and B T are individually divergence-free (solenoidal), and so is B naturally. This is one of the merits of the PT representation when used in numerical computation. Even if a numerical expression of B is not exactly divergencefree, the discretization error does not increase with time (iterations) nor accumulate in some places. Thus, we are relieved of the concern about ∇ · B. We can decompose a magnetic field B(r) into two parts as 3 The PT representation has also been called the Mie representation (Backus 1986) or the Chandrasekhar-Kendall representation (Montgomery et al. 1978;Low 2006) crediting Mie (1908) and Chandrasekhar & Kendall (1957), respectively. However, the credited works only addressed solutions of linear vector Helmholtz equations, to which linear FFFs also belong. For such a field, Φ and Ψ are not independent of each other while for a general magnetic field, Φ and Ψ are independent. The existence and uniqueness of Φ and Ψ for an arbitrary B was first treated by Backus (1958), and Chandrasekhar's description of general magnetic fields by two independent scalar functions was first given in his single-authored book (Chandrasekhar 1961). in which B n (r) =n(n · B), wheren = ∇ξ/|∇ξ|, is the component of B(r) normal to a constant-ξ surface containing the point r, and B t (r) = B −n(n · B) is the tangential component. Equation (8) tells that B T has only tangential components to the constant-ξ surfaces, i.e.,n · B T = 0 , and any field line of B T entirely lies in a coordinate surface of ξ. On the other hand, B P has both tangential and normal components in general, and This property is valid for an arbitrary scalar field ξ. However, not all formulations in the form of equations (6)-(8) with an arbitrary ξ are qualified to be called a standard PT representation, which additionally demands that the curl of a poloidal field be a toroidal field as the curl of a toroidal field is a poloidal field, i.e., where Θ is another scalar field. This requirement is met when the constant-ξ surfaces are either parallel planes or concentric spheres (Rädler 1974;Yi & Choe 2022). In other words, it is required that ξ = f (r), a function of r only in spherical coordinates, or ξ = f (z), a function of z only in Cartesian and cylindrical coordinates. Fortunately, stars, the sun included, are almost perfectly of a spherical shape, and the base of an individual active region can well be approximated as a plane. In a standard PT representation, equation (12) tells that and it follows that In the coronal FFF reconstruction, what is good about equations (11) and (14) is that B n and J n can be fully described by the functions Φ and Ψ in a constant-ξ surface only. In a Cartesian coordinate system, where is the 2D Laplacian operator in a z = const. plane. In a spherical coordinate system, where is the 2D Laplacian operator in an r = const. surface (sphere). Since the forms of ∇ 2 xy and ∇ 2 θϕ do not include any normal derivatives, B n and J n are fully described by the boundary values of Φ and Ψ only. Thus, implementing the boundary conditions B n and J n in a numerical grid is quite straightforward.
In this paper, we will confine ourselves to rectangular domains with Cartesian coordinates. In the vector potential description of magnetic field, where ∇ xy =x ∂ ∂x +ŷ ∂ ∂y is the 2D ∇-operator and The above equation shows that at least first order normal derivatives of two components of A are required for describing J z (z = 0). Thus, the boundary values of some components of A must be changed according to the variation of A inside the domain. This laborious adjustment of A(z = 0) can be avoided when the gauge is employed, not only at the boundary z = 0, but also inside the domain. The gauge (22) is satisfied when where Φ is an arbitrary scalar field defined in the domain. Thus, the vector potential under this gauge is in the form and the resultant magnetic field is in the form Comparing the above equation with equations (6)-(8), we can see that the Φ in equation (25) is nothing but the poloidal function Φ in equation (7) and the A z is −Ψ in equation (8). In this paper, the Φ above is our poloidal function, but for our toroidal function, A z is taken instead of −Ψ. Therefore, we have the following expressions for the respective three components of B and J , to be used for our numerical computation.
The first terms in the righthand side of equations (26) and (27) are the 2D curl-free (irrotational) part of the tangential magnetic field B t = B xy = B xx + B yŷ , and the second terms form its 2D divergence-free (solenoidal) part. The same is true for J t = J xy = J xx + J yŷ given by equations (29) and (30). A remark should be made that the PT representation in a Cartesian coordinate system demands a special treatment when the magnetic field is periodic in x and y (Schmitt & von Wahl 1992). In our FFF calculations of this paper, we do not employ a periodic boundary condition and thus are relieved of such a special care. See Appendix A for some details.
Before the present paper, the PT representation has been used in construction of force-free fields or magnetohydrostatic (MHS) equilibria. However, these works are related to LFFFs or a special class of MHS equilibria. For LFFFs, the toroidal function and poloidal functions are related by Ψ = αΦ, where α = const. and the poloidal function is a solution of a linear Helmholtz equation ∇ 2 Φ + α 2 Φ = 0. This linear problem was treated by an eigenfunction expansion (Chandrasekhar & Kendall 1957;Nakagawa & Raadu 1972) and by constructing Green's function based on eigenfunctions (Chiu & Hilton 1977). An MHS equilibrium is described by where p is the plasma pressure and ψ the gravitational potential. One can then express the current density in the form of J =αB + ∇µ × ∇ψ , in whichα = (J · ∇ψ)/(B · ∇ψ) and µ is a function of B · ∇ψ and ψ (see Low 1991, for details). In general MHS equilibria, it holds that B · ∇α = 0, i.e.,α is a function of each field line. In a special class of MHS equilibria, in whichα = const. everywhere, it holds that Ψ =αΦ, and the poloidal function Φ satisfies an equation ∇ 2 Φ +α 2 Φ − F = 0, where F is a certain function of µ (Neukirch & Rastätter 1999). This problem was solved with a Green's function method with an eigenfunction decomposition by Petrie & Neukirch (2000). Thus, all the previous works employing the PT representation have solved linear problems for the poloidal function only. Construction of NLFFFs employing independent poloidal and toroidal functions was proposed by Neukirch (1999), but has not been worked out in detail.

Numerical Algorithm
As pointed out by Grad & Rubin (1958), an FFF problem given by equations (1) and (2) is semi-elliptic and semi-hyperbolic. The hyperbolic nature of the problem lies in equation (3), which is called a magnetic differential equation (Kruskal & Kulsrud 1958). The Grad-Rubin type procedure includes a solver of the magnetic differential equation, practically putting more weight on the hyperbolic nature of the problem. In variational methods, the elliptic nature is more emphasized and the hyperbolic nature is rather implicitly considered. Our new algorithm is more inclined to variational methods in that we solve elliptic equations at every iteration step with the hyperbolic nature of the problem considered in the source terms. However, it is not variational because we do not try to extremize any functional. The rationale and details of our new method are given below.
Let us simply consider a magnetofrictional evolution of magnetic field under ideal MHD condition. Then, the electric field is given by E = −v × B after a proper normalization removing constants involved in unit systems.
where ν(r, t) is an arbitrary scalar field, which may be set to unity in nondimensionalized equations. Then in whichb = B/B is a unit vector in the direction of B, and J =b(b · J ) and J ⊥ = J − J are respectively component vectors of J parallel and perpendicular to B. Therefore, the induction equation (Faraday's law) reads Since the righthand side of this equation involves second order spatial derivatives of B, it is a parabolic equation.
As t → ∞, both sides of the parabolic equation go to zero. The asymptotic state ∇ × J ⊥ = 0 means that J ⊥ = ∇φ, where φ is a certain scalar field. If the field line footpoints at the boundary are fixed, E = 0 at the boundary and ∇φ = 0 everywhere in the domain. In the coronal FFF problems, the stressing (winding) and the relaxing (unwinding) are alternating or they indistinguishably coexist. In any case, as t → ∞, the stressing is reduced to zero, and the asymptotic state is such that J ⊥ = 0, i.e., J = J , a force-free state. Now the PT representation has two variables describing magnetic field, Φ and A z (= −Ψ), which are directly related to B z and J z by equations (28) and (31) respectively. From equation (34), we can write down the evolutionary equation of B z and J z as where t is a sort of "pseudo-time." Setting the righthand sides of these equations to zero (t = ∞) gives the following force-free conditions: It would be more than desirable to solve for B z and J z at once, but the righthand sides of the above equations are so highly nonlinear as to make such an attempt appear hopelessly difficult. The lefthand sides of the equations are linear and purely elliptic while the hyperbolic nature of the FFF equations is contained in the righthand sides. With this understanding, we propose the following iteration procedure: in which the superscript denotes the iteration step. Here we are to solve simple 3D Poisson equations for B z and J z with source terms evaluated with known values. Then, the poloidal function Φ and the toroidal function A z are updated by solving 2D Poisson equations at each z = const. plane, and all the other variables at the (m+1)-th step are calculated from Φ m+1 and A z m+1 by equations (26), (27), (29) and (30). The iteration procedure given by equations (39)-(42) is not variational because it does not try to extremize any functional. The iteration step m is far from representing any pseudotime because we are already at t = ∞. Although we currently cannot present a mathematical proof of the convergence of our algorithm, our tests with many different real and artificial magnetograms have never failed in convergence. A similar iterative method was used for solving 2D steady-state Navier-Stokes equations described in stream function and vorticity (Roache 1975).

Boundary and Initial Conditions
In this paper, we consider a rectangular domain Its boundary consists of six planes. The z = 0 plane corresponds to the coronal base, and the z = L z plane is an artificial top boundary of the corona. The lateral boundary of an active region corona consists of four planes: x = 0, x = L x , y = 0 and y = L y . A numerical grid is set up in the computational domain such that G = {(i, j, k)| i = 0, 1, 2, . . . , N x , j = 0, 1, 2, . . . , N y , k = 0, 1, 2, . . . , N z }.
At the bottom boundary (z = 0), B z and J z are given, which are readily translated into Φ and A z by solving 2D Poisson equations ∇ 2 xy Φ = B z and ∇ 2 xy A z = −J z once for all. Although imposing B z and J z at z = 0 is straightforward owing to the PT representation of magnetic field, evaluating the righthand side of equation (41) at the grid plane k = 1 is somewhat tricky because the finite difference form of the z-derivative term ∂ ∂z ∇ xy ·J at k = 1 (the first grid plane above the bottom boundary) requires the value of J at k = 0 (z = 0). As can be seen in equations (29) and (30), J x and J y involve a second order z-derivative ∂ 2 Φ/∂z 2 , which cannot unambiguously be defined at z = 0. However, the term J (z = 0) does not need to be evaluated at all. Since the equations we are solving (eqs. [37] and [38]) already represent a stationary state (t = ∞), we should rather set J ⊥ (z = 0) = 0. Using the relation ∇ · J = −∇ · J ⊥ , we can rewrite equation (38) as The stationary condition J ⊥ (z = 0) = 0 allows us to evaluate the second term in the rightmost hand side of the above equation straightforwardly. The implementation of the boundary condition at z = 0 is thus completed. We have tried another method, which uses a plausible relation in which α = J z /B z and B * xy is the observed horizontal magnetic field. Using this relation corresponds to employing the observed three components of magnetic field at z = 0 and enforcing the force-free condition at the bottom boundary. The results of both methods are compared for many cases and they turn out to be surprisingly similar to each other.
The magnetogram of an active region generally does not have the same positive and negative magnetic fluxes because some field lines are connected to the outside of the active region as closed fields or to the solar wind as open fields. We want to accommodate such a flux imbalance in the new FFF construction method. Previous attempts to consider a flux imbalance have placed image polarities of opposite signs outside the real domain (Seehafer 1978;Otto et al. 2007). These methods, however, may make a considerable flux go out of or come into the real domain through the lateral boundary and thus the constructed field configuration may undesirably depend on the field connection between the real domain and the image domains rather than on the field connection within the real domain alone. In our study, we rather confine all the flux within the computational domain except the flux passing through the top boundary. Magnetic field lines thus cannot penetrate the lateral boundary, but can be tangential to it. The effects of other boundary conditions allowing nonzero magnetic flux through the lateral boundary are discussed in Section 5.
In regard to the lateral boundary (x = 0, L x , y = 0, L y ), the unit outward normal vector to the boundary is denoted byn, the normal component vector of a vector field F by F n =n(n · F ), its tangential component vector by F =n × (F ×n) = F − F n , and the normal derivative by ∂/∂n =n · ∇. Since we have to solve four Poisson equations, i.e., equations (39)-(42), at every iteration step, we have set up four boundary conditions respectively on B z , Φ, J z and A z at the lateral boundary. To make this boundary impenetrable to magnetic field, i.e., B n = 0, we choose the following boundary conditions.
As can be seen in equations (26) and (27), the conditions (46) and (48) result in B n = 0. The 2D Poisson equation (40) gives the following constraint: in which dl is a unsigned line element surrounding the 2D domain at z = const. and dS is the area element in the plane. Then, p 1 in equation (46) is given by If an infinite plane is a perfectly conducting rigid wall impenetrable to magnetic field, the magnetic field in one side of the wall can be regenerated by replacing the wall by an image field beyond the wall in such a way that B t is symmetric and B n is anti-symmetric across the boundary. In our case, the lateral boundary consists of four finite planes and the magnetic flux in each z = const. plane is unmatched, which does not realize an exact symmetry or anti-symmetry across each boundary plane. Nevertheless, the local symmetry condition (45) turns out to be in practice a good choice for our situation too. The image field beyond an infinite wall also gives J t = 0 and ∂J n /∂n = 0. The boundary condition (47) is also inspired by this picture. At the top boundary, we employ a source surface condition (Altschuler & Newkirk 1969;Aly & Seehafer 1993). The purpose of this boundary condition is not only to mimic the upper corona, where the solar wind carries out the open magnetic field, but also to provide an exit for the imbalanced magnetic flux. The source surface boundary condition is simply given by The condition (51) results in implying that no current can escape or come in through the top boundary. Since J n = 0 at the lateral boundary in the final force-free state, it follows that z=0 J z dx dy = 0 , i.e., the positive and negative currents through the bottom boundary should be balanced. This constraint must be fulfilled by preprocessing of the magnetogram data. The source surface condition at the top boundary is known to create a magnetic field null-line at the boundary and a current sheet connected to it in the domain (Zwingmann et al. 1985;Platt & Neukirch 1994) except for potential fields. This point will be discussed in relation to our numerical solutions in Section 5.
To start the iteration procedure given by equations (39)-(42), we need zeroth step (n = 0) values of B z , Φ, J z and A z , which may be called initial conditions. So far many numerical calculations for FFFs have taken potential fields for their initial conditions (Schrijver et al. 2006). In the PT representation, a potential field is prescribed by in which the last line results from equation (55c) and our lateral boundary condition (48). Since equations (55c) and (55d) do not match our boundary conditions at z = 0, the potential field as an initial condition will have awkward bending of field lines just above the bottom boundary. Instead of starting with a current-free state, we rather want to load the entire domain with currents. The simplest way for this would be to set This field is prescribed by with ∇ 2 xy Φ = B z and ∇ 2 xy A z = −J z . In our experience, this initial field reduces the number of iteration steps required for convergence compared with the potential field as an initial condition.  a The field line twist at the surface of the entire closed torus, which is larger than the twist of the coronal part of the torus. b Used in Wiegelmann et al. (2006a).
As with other numerical methods, our method is to be tested with a computational code built upon it against known analytical solutions. So far most existing codes have been tested against the analytical model by Low & Lou (1990), which represents a class of moderately sheared magnetic fields, and have shown good to excellent performances (Schrijver et al. 2006). Our own variational code VFVP (Variational FFF Code in Vector Potential Formulation) we developed earlier (see Appendix B) and the code based on our new method NFPT (Non-variational FFF Code in Poloidal-Toroidal Formulation) also are found to be working as well for the Low & Lou model as other codes are (Schrijver et al. 2006;Inoue et al. 2014), and we do not feel the necessity of presenting the results here. In this paper, we present the performance of the codes in reproducing the NLFFF model by Titov & Démoulin (1999). In this analytical model, a flux rope carrying a helical force-free field and ambient fields generated by two magnetic charges and a line current underneath the surface are in equilibrium. Each Titov-Démoulin (TD) model is prescribed by six free parameters: R the major radius of the torus, a the minor radius, L the half-distance between two magnetic charges, q the magnitude of the magnetic charge, d the distance between the photosphere and the subsurface line current lying along the symmetry axis of the torus, and I 0 the subsurface line current. The toroidal current I in the torus is determined by the force balance condition. The field line twist number N t at the surface of the entire (closed) torus also comes out from these parameters. The coronal part of the torus, which lies above the photospheric plane, has a twist N cor at its surface Table 1 lists the parameters prescribing the three TD models to be presented in this paper. It is to be noted that since the TD model employs a thin flux tube approximation assuming a high aspect ratio (A = R/a 1), the minor radius of the flux tube may not be uniform along its axis in numerical models with a modest A, particularly when lateral and top boundary conditions are differently prescribed from those of the analytical models.
Comparison of numerical results and reference models can first be made by their appearance. In each TD model, there is only one flux rope involved. It is thus a touchstone of different numerical models whether one flux rope is well reproduced with the characteristic features of the analytical model. Figure 1 shows the field lines of the flux rope in model TD1 (refer to Table 1) obtained from (a) the analytical model, (b) our new code NFPT, (c) our earlier code VFVP, and (d) the optimization code in the SolarSoftWare (SSW), which is available in the public domain 5 . All field lines are traced from the same footpoints, one group in the positive polarity area and the other in the negative polarity area of the flux rope in the analytical model. The field lines traced from the positive and negative footpoints of the magnetic axis of the flux rope in the analytical model are rendered thicker than other field lines. The field lines of NFPT and VFVP apparently express a single flux rope while those of SSW show two separated flux tubes. The two magnetic axis field lines obtained from our new code NFPT overlap with each other as in the analytical model, but those from our earlier code VFVP are slightly off if not as much as those from SSW. Since we do not have observational data at the lateral and top boundaries in practical situations, all numerical computations here use the data of the analytical model at the bottom boundary only, and proper artificial boundary conditions are used for the lateral and top boundaries. For the same TD model, Wiegelmann et al. (2006a) also reported two separated flux tubes crossing each other when the analytical solution is used only at the bottom boundary, whereas their constructed field bears a remarkable resemblance to the analytical model when the analytical solution is imposed at the lateral and top boundaries too.
Beyond mere appearance, we can quantitatively compare numerical models with analytical solutions using so-called "figures of merit" devised by Schrijver et al. (2006), which will also be called "performance metrics" in this paper. Although we will use the same notations for them as in Schrijver et al. (2006), the definitions of those and other metrics are reiterated here for readers' convenience. are not known to those who construct B. Also, the analytical exactness is different from numerical exactness in a grid of finite resolution. Furthermore, the analytical model of Titov & Démoulin (1999) is an approximate solution. For these reasons, some authors have (e.g., Valori et al. 2010;Guo et al. 2016) relaxed the analytical solution in a numerical grid to obtain a reference field. In this paper, however, we use the unprocessed analytical solutions as reference fields in spite of the possible loss of performance scores. The metric CW sin was originally devised by Wheatland et al. (2000), in which a different notation σ J was given to it, but the current naming has become more popular afterwards. While the metrics C vec , C CS , E n , E m , and are based on the comparison of the numerical solution with the reference field, the metric CW sin only judges the force-freeness of the numerical solution itself. If the reference field and the numerical solution are exactly identical, the former five metrics should be all 1. For an exact force-free field, CW sin should be 0. In some literature and also in this paper, another quantity p is presented, which is the magnetic energy of a solution in units of the potential field energy.
where B p the potential magnetic field. This metric uses the potential magnetic field as the reference field and indicates the nonpotentiality of the solution, which may be due to the deviation of the solution from the exact FFF as well as to the difference of the exact FFF and the potential field. Table 2 lists the seven performance metrics for TD1 derived from the analytical solution and the three numerical solutions, NFPT, VFVP and SSW. The presented value of CW sin for the analytical solution is obtained by assigning the analytical magnetic field to a numerical grid and evaluating the current density by a finite-difference method, and hence it deviates from 0 inevitably. As for the first five metrics, the scores are ranked in order of VFVP, NFPT and SSW. However, we have already seen that the field obtained by VFVP does not show as much resemblance to the analytical solution as that by NFPT. It is only in TD1 that VFVP's scores are better than NFPT's. For other TD models, NFPT turns out to exceed VFVP in all metrics, as we shall see further below (Tables 3 and 4).. In contrast to the metrics by Schrijver et al. (2006), the scores of the metric CW sin are ranked in order of NFPT, VFVP and SSW, reflecting best the resemblance of the numerical solutions to the TD field featured by one flux rope. The metric CW sin is independent of the choice of the lateral and top boundary conditions and evaluates the force-freeness of a solution itself. It is interesting that the apparent resemblance depends more on the exactness (force-freeness) of the solution than its numerical proximity to a reference solution. In other TD models too, the NFPT code conspicuously outperforms others in CW sin.  Figure 1 shows that field lines are more and more loosely wound as we go from (b) to (d). To compare the performance of different numerical methods, we have also compared the twist of the numerical solutions with that of the analytical one. The twist of a field line denoted by C B is only meaningfully defined about another field line denoted by C A , and we thus denote the twist of C B around C A by T w(C A , C B ). In the case of a flux rope, its magnetic axis would be the most meaningful choice for C A . In numerical solutions, the magnetic axes traced from the positive polarity and from the negative polarity can be different if more than one flux tube appear in a numerical solution as in the cases of VFVP and SSW. Thus we have to calculate two sets of twist  with two magnetic axes traced from different polarity areas. The twist is solely defined by serial local connections between two curves, and the following formula is valid for open curves as well as for closed curves (Tyson & Strogatz 1991;Berger & Prior 2006).
in whicht A = dr A ds , where s is the arclength of the curve C A , is a unit tangent vector to C A , andû = u/|u|, where u is a vector connecting a point in C A to a close point in C B in such a way that u ·t A = 0. With this method, we have measured the twist of 24 field lines, whose footpoints are in a circle of radius r centered at a magnetic axis footpoint at z = 0 , and have taken an average over them.
in which (r, θ i ) is the polar coordinates of the i-th field line footpoint, and θ i = iπ/12. Circles of five different radii, r/r 0 = 1, 2, 3, 4, 5 where r 0 = 7 Mm, are taken in the positive and negative polarity areas, respectively, so that ten mean values of twist in units of turns (2π radians) are obtained for each solutions and are plotted in Figure 2. As shown in the figure, the mean twists of field lines traced from the positive polarity side and from the negative side for the same r are indistinguishable. It is remarkable that the twists of the NFPT solution are closest to those of the analytical solution. The twists obtained by VFVP fall a little short of the twists of the former two solutions, and the SSW yields quite small twists as seen in Figure 1. Previously, the twist of a field line was sometimes measured by the following formula (Inoue et al. 2011(Inoue et al. , 2014 in which α = B · J /B 2 and and l is the arclength of the field line C B . This method is possibly subject to criticism in that a twist cannot be defined for one curve alone. With an implicit assumption that a magnetic axis of a flux rope takes the role of C A , the equation can be made meaningful with a correction that the line integral should be taken over the magnetic axis C A , dl being the arclength of C A , and necessarily under the condition that the toroidal current density should be uniform over each cross-sectional area of the flux tube where the field line C B lies. We have compared the twists by equation (66) and those by equation (68) without modification for all the solutions of TD1, and have found that the latter method tends to overestimate the twists, which is attributed to the larger length of C B than C A . The dis- crepancy is found to grow with r, because the length ratio of C B to C A increases with r. The model TD2 (see Table 1) characteristically has a bald patch, which is a segment in a polarity inversion line, where locally concave upward field lines touch the photosphere (Titov et al. 1993). A bald patch is a possible location of a solar prominence of inverse polarity (Lee et al. 1995;Mackay et al. 2010), and it can develop into a current sheet (Low 1992;Cheng & Choe 1998), where magnetic reconnection may take place (Delannée & Aulanier 1999). In TD models, a bald patch appears only with a sufficiently large twist of the flux rope (Titov & Démoulin 1999). TD models with bald patches have been reproduced in several numerical solutions (Valori et al. 2010;Jiang & Feng 2016;Demcsak et al. 2020).
The parameters of our TD2 are the same as those of Case 5 in Demcsak et al. (2020), in which a Grad-Rubin type code (Wheatland 2007) was used. Figure 3 shows the top and side views of four TD2 solutions, (a) analytical, (b) NFPT, (c) VFVP and (d) SSW. As in the case of TD1, the NFPT solution reproduces one flux rope with a single magnetic axis of the analytical solution. The solutions by VFVP and SSW produce two flux tubes. The two flux tubes by SSW are widely separated. All four models show bald patches, but the endpoints of the bald patch field lines are quite differently located. The footpoint positions of the bald patch field lines apparently indicate their writhes (Berger & Prior 2006). The footpoints and the S-shape appearances of the bald patch field lines in the NFPT solution and the analytical solution are quite close although the field lines in the latter are slightly longer and their apex reaches a higher altitude. While the separation of the magnetic axes of two flux tubes is smaller in VFVP than in SSW, the writhe of the bald patch field lines is larger in SSW than in VFVP. Compared with the solution by Demcsak et al. (2020), the analytical solution and and the NFPT solution show larger writhes of bald patch field lines than their result while the VFVP and SSW solutions show smaller writhes. Limited to TD models, nonvariational methods (NVPT and a Grad-Rubin type code) seem to be better at revealing topological features of the field than variational methods (VFVP and SSW).
The performance metrics of the four solutions for TD2 are listed in Table 3. The NFPT solution exceeds other numerical solutions in all metric scores, and it particularly stands out in CW sin. It is again suggested that the force-freeness represented by CW sin is more important in reproducing topological features of the fields, e.g., flux ropes and bald patches, than other metrics.

Hyperbolic Flux Tubes in the Titov-Démoulin Models
If a magnetic field is purely two-dimensional, i.e., if it depends on two coordinates and its field lines lie in parallel planes, an X-shaped field configuration has two surfaces called separatrices intersecting each other at a line called a separator. The field lines in each quadrant may locally be approximated by hyperbolas near the Xpoint and they have discontinuous connectivities across separatrices. The separator may be deformed into a current sheet, where magnetic reconnection can take place  (Syrovatskii 1981). If we add a smooth third component of magnetic field dependent on two coordinates to a purely 2D magnetic field, two field lines, whose endpoints are very close with a former separatrix between them, have a large separation at the other ends (Longcope & Strauss 1994). Thus, the field connectivity is everywhere continuous. In a 3D situation, separatrices and separators appear if nullpoints or bald patches are present in the domain. Otherwise, the field connectivity is everywhere continuous, but there may exist quasi-separatrix layers (QSLs), across which the separation between the ends of field lines, which are very close somewhere in the domain, is very large (Priest & Forbes 1992;Démoulin et al. 1996). Quite similarly to the 2.5D situation, the field lines in the vicinity of the intersection of two QSLs look like hyperbolas when they are projected onto the plane normal to the intersection. This structure is called a hyperbolic flux tubes (HFTs) and its magnetic axis is a quasi-separator (Titov et al. 2002(Titov et al. , 2003. A current sheet is likely to be formed in an HFT by magnetic pinching induced by suitable footpoint motions and consequently magnetic reconnection ensues (Titov et al. 2003). To find an HFT, Titov et al. (2002) proposed a squashing factor (or degree) Q, which represents how much a cross-section of a flux tube is deformed at the other side, as follows.
in which [X ∓ (x ± , y ± ), Y ∓ (x ± , y ± )] is a vector function connecting a footpoint (x ± , y ± ) to its conjugate footpoint (x ∓ , y ∓ ). The flux surface of an HFT is an isosurface of Q with Q 2. To see how well numerical FFF solvers reproduce a TD field containing an HFT, we have constructed the model TD3 in Table 1 and calculated the squashing factor Q by tracing field lines to the boundary from every grid-point in the domain. The upper row of Figure 4 shows the Q value in the y = 0 plane, which is a vertical plane cutting the TD flux rope at its apex. Apparently the Q color maps show an X-shaped configuration, which is a typical feature of an HFT. The X-shaped structure manifested by the same color scheme of Q is most conspicuous in the analytical solution, then in the NFPT solution and rather faint in the SSW solution. Figure 4 also shows field lines passing through three boxed areas in the y = 0 plane, a black box above the axis of HFT, a light green box in the vicinity of the HFT axis and an orange box just below it. The upper field lines above the HFT form one flux rope in the analytical and NFPT solutions and two flux ropes in the SSW solution. The light green field lines are supposed to cover a part of the QSL. The structure of the QSL is best revealed in the analytical solution and then in the NFPT solution. Figure 5 shows the Q value in the z = 0 plane in black and white brightness, superposed with the B z level contours, positive in red and negative in blue. In all three solutions, the Q distributions below the central part of the flux ropes (−1 < y < 1) are quite similar, but in the vicinity of the flux rope legs, the Q distributions are noticeably different. The Q map of the analytical solution has hook-like structures, which almost surround the flux rope footpoint areas. The result of our NFPT also shows hook-like structures, which, however, are a little smaller and thinner than the analytical ones. In the Q map by SSW, the hooks are not fully wound, but are broken in the middle. Apparently the difference in the Q maps of TD3 is somewhat similar to the difference in the bald patch field lines of TD2 shown in Figure 3. To assess the proximity of the overall Q distribution of a numerical model to the analytical one, we devise the following metric where i is the grid-point index, M the total number of grid-points, q the squashing factor of the analytic model, and Q the squashing factor of a numerical solution. If the numerical solution is identical to the analytical one, the metric E q should be zero. We have evaluated the values of E q in the z = 0 plane (E q,bot ) as well as in the whole computational domain (E q,whole ) and have listed them in Table 4 along with the performance metrics. The error metric E q,bot of NFPT is about one fifth that of SSW. However, the error metric E q,whole does not show so much difference as E q,bot , which may be at-tributed to the larger portion of grid-points with small q (Q also) in the whole domain than in the bottom plane only.
In an HFT, the Q value takes a maximum in a line called a quasi-separator (Titov et al. 2002). It corre- b The error metric for the Q values in the whole computational domain. sponds to an X-line in a 2.5D magnetic field and may be called the magnetic axis of an HFT. If one purposes to find a quasi-separator (QS) and the HFT in its neighborhood rather than to find the entire structure of a QSL, one may not want to take the trouble of evaluating Q, which requires finding the conjugate footpoint pairs by field line tracing. Here we propose a rather simple "local" method of locating a QS and HFT. In a plane normal to a QS, field lines projected onto this plane are like hyperbolas and the magnetic field in this plane has the following property, in which x and y are arbitrary Cartesian coordinates in this plane. Since we do not know the location of the QS from the beginning, we just calculate J ⊥ in a plane normal to the local magnetic field. Then, we set so that we may focus on candidate points for a QS. Although the isocontours of H Q would not follow the shape of a QSL, the maximum of H Q is expected to be located on the QS. Figure 6 shows a color map of H Q in the y = 0 plane made from the analytical solution of TD3 superposed with white isocontours of Q. While Q contours are useful for identifying a QSL structure, an H Q map (or contours) is advantageous for locating a QS. As both Q and H Q are useful for finding an HFT, their distributions in the bottom plane (z = 0) are expected to be similar. This is confirmed by Figure 7, in which the maximum value of H Q in a field line is expressed by brightness of brown color at its footpoint in z = 0 and the isocontours of Q = 100 are given in blue for the analytical solution, the NFPT solution and the SSW solution of TD3. The figure shows that the distributions of the two quantities are quite consistent. For practical purposes, we do not need to find the maximum of H Q is each field line. One may roughly draw field lines globally and construct an H Q distribution in a suspected region. Then, H Q can be a quick tool for locating a QS and an HFT.

Efficiency and Effectiveness of the Computational Code
Our NFPT code requires solving two 3D Poisson equations (39) and (41) and two 2D Poisson equations (40) and (42) for k = 1, 2, ..., N z planes in each iteration step. In the code, the Poisson equations are solved by a direct solver package FISHPACK (Swarztrauber & Sweet 1975) quite efficiently and accurately. The computational time per iteration step required by NFPT is about ten times that of VFVP and about twice that of SSW. However, the number of iteration steps required for convergence by NFPT is about one tenth of that by VFVP or SSW. In our TD2 calculations, the NFPT code re-quired about 200 iteration steps while the other codes demanded at least 2500 steps. In terms of the computational resources required for producing a solution, the NFPT is comparable to the VFVP and uses far less resources than the SSW.
A problem described by a set of partial differential equations is well-posed if the imposed boundary conditions and constraints yield a unique solution. Although a problem is mathematically well-posed, its numerical implementation may not be so for diverse reasons. Among others is the way of information transfer between the boundary and the inner computational domain, inherent to the numerical method. In order to see whether the boundary condition is well reflected in the solution consistently throughout the computational domain, we have devised the following test. Suppose that a numerical solution has been constructed in a domain Then, one can construct a numerical solution in a domain reduced in z, i.e., in with a numerical grid G k0 = {(i, j, k) | i = 0, 1, 2, . . . , N x ; j = 0, 1, 2, . . . , N y ; k = k 0 , k 0 + 1, k 0 + 2, . . . , N z }, where k 0 ∆z = z 0 with ∆z = L z /N z , using the former numerical solution at z = z 0 (k = k 0 ) as a boundary condition at the bottom boundary of the new domain. At this time, we should not use the former solution as the initial condition, nor impose any known solutions as the boundary conditions at the lateral and top boundaries of G k0 . If the numerical solutions with grids G 0 and G k0 are identical in the domain G k0 , it can be said that the bottom boundary and the solution in the domain are consistently connected. For TD2, we have employed a numerical grid with (N x , N y , N z ) = (150, 250, 100). To perform the above test, we have chosen k 0 = 10, 20, 30, 40, 50 and constructed numerical solutions with the NFPT code in five different grids G k0 , in all of which the numerical solutions obtained in the original grid G 0 are used as bottom boundary conditions at k = k 0 of G k0 . To compare the numerical solutions in the original grid and in each reduced grid, we have evaluated five performance metrics by Schrijver et al. (2006) in the domain G k0 , using the solution obtained in each reduced grid G k0 for B and the solution in the original grid G 0 for b in equations (59)-(63). The results are given in Table 5. The metrics C vec and C CS are 0.999 for all k 0 while other metrics show a slight tendency of degradation with increasing k 0 . These scores of NFPT are unrivaled with those of variational codes (VFVP and SSW), all the more so with increasing k 0 , unless fixed values of the field are prescribed at all boundaries. Thus, it can be said that the PT representation and its natural way of imposing boundary conditions are self-consistent to generate an unambiguous solution in the whole domain. It is also suggested that this test be performed on other numerical methods to be developed.

Discussion on Lateral and Top Boundary Conditions
The lateral boundary condition we have used in this paper does not allow any magnetic flux to cross the lateral boundary. Any imbalanced flux at the photospheric boundary is thus ducted through the top boundary. Previously, other boundary conditions have been proposed, which allow nonzero normal magnetic field at the lateral boundary (e.g., Seehafer 1978;Otto et al. 2007). These boundary conditions assume image polarities across the lateral boundary by a plane symmetry or a line symmetry of a certain parity so that there may be field lines Note-For all the metrics, the numerical solution in the grid G0 is taken for the reference solution b and the numerical solution in the grid G k 0 is taken for B.
connecting the polarity patches in the real domain and those in the image domains. Thus, the field configuration in the real domain is more or less affected by the field connection across the lateral boundary. Since there are several methods of placing image polarities, the resulting field configuration depends on the method chosen (Otto et al. 2007). By mathematical reasoning alone, one can hardly tell which boundary condition is superior to others . Our non-penetrating wall condition, which is also an artificial boundary condition, is not an exception. As far as a more or less isolated active region is concerned, however, the connections between the  Note-The total unsigned magnetic flux through the bottom boundary is denoted by F b , the total unsigned flux through the lateral boundary by F l , and the total unsigned flux through the top boundary by Ft. The total open flux is F l + Ft. The three flux ratios obtained in a domain with a double size in z are designated by the subscript d.
real (observed) polarities and the image polarities may rather be considered undesirable. For TD2, we have tested four different lateral boundary conditions, one of which is our non-penetrating condition, in order to see how much magnetic flux passes through the lateral and top boundaries and how close the numerical solution is to the reference field. The four lateral boundary conditions are as follows.
In BC1, which is taken in the present paper, B n = 0, and J n is not necessarily zero, but must be zero in a force-free state. In all the other cases, B n = 0. In BC2, B t = 0, J n = 0, and J t is not necessarily zero, but must be zero in a force-free state. Thus, the magnetic field passing through the lateral boundary in BC2 is a potential field. For BC1 and BC2 both, the force-free coefficient α is anti-symmetric across the lateral boundary, and neither condition can be accommodated in LFFF models. Since there is no current flowing through the lateral and top boundaries with BC1 and BC2, the current at the bottom boundary must be preprocessed so that the net current may be zero. In BC3 and BC4, there are normal magnetic currents as well as normal magnetic fields crossing the lateral boundary. Since the force-free coefficient α is symmetric across the lateral boundary, these boundary conditions can be used for LFFF model as BC3 was adopted by Seehafer (1978). The line symmetry boundary condition by Otto et al. (2007) has not been tried here. In Table 6, the total unsigned magnetic flux through the lateral boundary F l , the total unsigned flux through the top boundary F t , the total unsigned open flux F open = F l + F t are given in units of the total unsigned flux through the bottom boundary F b . Also, performance metrics are given for each case. The total unsigned open flux is the smallest with BC1, and the lateral flux of other cases is almost comparable to that. The performance metrics are best with BC1, which implies that the numerical solution is closest to the reference field among other cases. Since the lateral magnetic flux is also expected to depend on the domain's aspect ratio, we have also tried to run the code in a domain, whose height is twice the original L z while the x-and y-sizes are unchanged. The ratios of boundary fluxes to the bottom flux in this case are given in the last three columns of Table 6. The flux through the top boundary is significantly reduced for all cases compared with the value in the original domain. However, the total open flux through the lateral and top boundaries is either slightly decreased (BC3 and BC4) or even increased a little (BC2). With BC1, the total open flux is significantly decreased by doubling the domain size. Although not given in Table 6, the performance metrics are best with BC1 in the double-sized domain too. Based on this investigation about the boundary conditions and the code performance, it is carefully suggested that our boundary condition (BC1) would be the best choice for reconstruction of a rather isolated active region and that it would relieve one from worrying about where and how to put image polarities and current sources outside the domain.
It should be noted that the TD models have a perfect flux balance at the bottom boundary and the open flux is caused by the boundary conditions, not by a flux imbalance. When we use a source surface condition given by equations (51) and (52) at the top boundary, there must be a polarity inversion line where B z = 0 at this boundary unless the bottom magnetic field is totally unidirectional. From the source surface condition and the force-free condition, J (z = L z ) = 0. Thus, the open field is a potential field. A surface emanating from the null-line to the bottom boundary is a separatrix sur- face, across which field topology discontinuously changes from "closed" to "open" as well as "non-potential" to "potential." Since there is no reason for the field direction to be continuous across this separatrix, this surface must be a current sheet (Zwingmann et al. 1985;Platt & Neukirch 1994). As we do not know in advance where the footpoints of the separatrix are located at the bottom boundary, the footpoint line may lie where J z = 0 unless we have a wide enough buffer area with J z = 0 and B z = 0 at the bottom boundary. This is an intrinsic problem of our model, which causes those performance metrics evaluating the proximity of the numerical solution to the reference field to deviate from 1 even though the solution in the closed flux region is quite similar to the reference field. The metric CW sin measuring the force-freeness, however, is minimally affected except near the separatrix. To see whether the null-line with a singular transverse current exists at the top boundary and whether a separatrix current sheet exists in the domain, we have plotted the distributions of B z and J = |J | at the top boundary and that of J in the vertical plane x = 0 in Figure 8. Because the value of numerically evaluated current density highly depends on magnetic field strength and the null-line and the separatrix are in a weak field region, the features we are seeking do not stand out conspicuously, but they are still identifiable in the figure. As we have expected, the major field structure with strong currents lies in the closed flux region, and we have thus been able to reproduce the significant geometrical features of the reference TD fields as given in the previous section.

Summary
In summary, we have presented a new method of constructing a coronal force-free field based on a poloidaltoroidal representation of magnetic field. The PT representation allows us to impose the boundary conditions B n and J n at the photospheric boundary once and for all with only the boundary values of the poloidal and toroidal functions. With a rigid, conducting, slip wall boundary conditions at the lateral boundaries and a source surface condition at the top boundary, a magnetic flux imbalance at the bottom boundary can be accommodated. Since no current can escape the computational domain, the current at the bottom boundary must be preprocessed so that the net current through it may be zero. At the top boundary, a rigid, conducting, slip wall condition can also be used instead of the source surface condition. With this condition, however, a flux imbalance at the bottom boundary is not allowed and thus the B n data there must be preprocessed. Our new method is nonvariational in the sense that the converging sequence toward the solution does not extremize any conceivable functional. It rather directly targets a solution, but iterations are needed because of the high degree of nonlinearity. Thus, it requires far fewer iteration steps than variational methods although one iteration step requires more computational resources than the latter.
We have tested the NFPT code based on the new method against the analytical FFF models by Titov & Démoulin (1999) with other available variational codes, our own VFVP code and the optimization code in So-larSoft (SSW). Our new NFPT code excels relative to others in reproducing characteristic features of TD models, for example, one flux rope with a proper twist, a bald patch with a proper writhe and quasi-separatrix layers with a hyperbolic flux tube. The NFPT code also produces the best scores in most performance metrics (Schrijver et al. 2006), especially in CW sin measuring the solution's own force-freeness. The application of the NFPT code has also been made to the vector magne-tograms of a real active region, which will be reported in a sequel paper shortly.

ACKNOWLEDGMENTS
The authors thank the reviewer for his/her constructive comments, which have helped improve the paper a lot.
This work has been supported by the National Research Foundation of Korea (NRF) grants 2016R1D1A1B03936050 and 2019R1F1A1060887 funded by the Korean government. It was also partially supported by the Korea Astronomy and Space Science Institute under the R&D program (Project No. 2017-1-850-07). Part of this work was performed when S. Y. and G. S. C. visited the Max-Planck-Institut für Sonnensystemforschung, Göttingen, Germany, and they appreciate the fruitful interactions with researchers there.

A. ON THE PT REPRESENTATION OF A PERIODIC FIELD IN CARTESIAN COORDINATES
Here we will plainly expound why a special treatment is needed in a PT representation in a Cartesian coordinate system when the magnetic field is periodic in two directions spanning a toroidal field surface. Our account here is a readable revision of previous studies (e.g., Schmitt & von Wahl 1992, and references therein).
Consider a magnetic field periodic in x and y with wavelengths (periods) λ x and λ y , respectively. Taking ξ = z, we have We may define two vector fields F and G such that Note that F and G are also periodic in x and y. We can then write B y =ẑ · (∇ × G) .
We consider a 2D rectangular domain S z = {(x, y, z)|0 ≤ x ≤ λ x , 0 ≤ y ≤ λ y } and apply Stokes' theorem to have In the line segments x = 0 and x = λ x , F (x, y, z) or G(x, y, z) is the same, but dr is in the opposite direction. The same is true for the line segments y = 0 and y = λ x . Thus, the contour integrals should be zero. However, the surface integrals may not be zero even if B is periodic. Such a field cannot simply be expressed in the form of equations (A1) and (A2). In such cases, we define the following functions of z only B 0y (z) = 1 λ x λ y Sz B y dS z .
Then, the surface integrals of B x − B 0x and B y − B 0y are zero, and only B − B 0xx − B 0xx can be expressed by the standard PT representation in the form of equations (A1) and (A2). A similar argument applies to B z given by equation (28). Applying the divergence theorem to the 2D domain S z , we have in which dl = |dr|. For Φ periodic in x and y, the contour integral is zero, but the surface integral of B z may not be zero even if B z is periodic. We then define Note that B 0z is not a function of z, but a constant owing to the constraint ∇ · B = 0 combined with the periodicity of B in x and y.
To sum up, a magnetic field B periodic in Cartesian coordinates x and y is generally expressed as in which B 0 = B 0x (z)x + B 0y (z)ŷ + B 0zẑ (A14) is so called a "mean flow" (Schmitt & von Wahl 1992). When periodic boundary conditions are explicitly imposed in two Cartesian coordinates, one should consider the mean flow. Since we do not employ a periodic boundary condition, a mean flow is not a consideration in our paper.

B. ON OUR EARLIER VARIATIONAL NLFFF CODE IN VECTOR POTENTIAL FORMULATION
Before the NLFFF code based on a PT representation, we had developed and used a variational NLFFF code using a vector potential formulation of magnetic field. Since this code is based on a magnetofrictional method (Chodura & Schlüter 1981), the algorithm is as simple as in which t is a pseudo-time, ν(r, t) a proper coefficient maxizing the convergence rate, and J ⊥ = B × (J × B)/B 2 . To expedite the convergence, we equip the code with a gradient descent algorithm (Chodura & Schlüter 1981). When a vector potential A is used to describe the magnetic field, we cannot set all three components of A at z = 0 fixed in order to impose B z and J z there. In a magnetofrictional code by Roumeliotis (1996), A z (x, y, 0) was set fixed and A x (x, y, 0) and A y (x, y, 0) were varied at every time-step. In our variational code, we set A x (x, y, 0) and A y (x, y, 0) once for all to fix B z (x, y, 0) and the solution of the following 2D Poisson equation is given as A z (x, y, 0) at every time-step, in which J z,obs is the boundary condition of J z derived from an observation and the z-derivative is evaluated using a one-side finite differencing. The computed J z (x, y, 0) is thus equated with J z,obs at every time-step. Our variational code is working very well for moderately sheared fields (e.g., Low & Lou 1990), but shows a little weakness for magnetic fields with flux ropes as with other variational codes. This motivated us to devise the new formulation presented in this paper, in which the imposition of the bottom boundary condition is tidy and effective.