Structure of pressure-gradient-driven current singularity in ideal magnetohydrodynamic equilibrium

Singular currents typically appear on rational surfaces in non-axisymmetric ideal magnetohydrodynamic (MHD) equilibria with a continuum of nested flux surfaces and a continuous rotational transition. These currents have two components: a surface current (Dirac δ-function in flux surface labeling) that prevents the formation of magnetic islands, and an algebraically divergent Pfirsch–Schlüter current density when a pressure gradient is present across the rational surface. On flux surfaces adjacent to the rational surface, the traditional treatment gives the Pfirsch–Schlüter current density scaling as J∼1/Δι , where Δι is the difference of the rotational transform relative to the rational surface. If the distance s between flux surfaces is proportional to Δι , the scaling relation J∼1/Δι∼1/s will lead to a paradox that the Pfirsch–Schlüter current is not integrable. In this work, we investigate this issue by considering the pressure-gradient-driven singular current in the Hahm–Kulsrud–Taylor problem, which is a prototype for singular currents arising from resonant magnetic perturbations. We show that not only the Pfirsch–Schlüter current density but also the diamagnetic current density are divergent as ∼1/Δι . However, due to the formation of a Dirac δ-function current sheet at the rational surface, the neighboring flux surfaces are strongly packed with s∼(Δι)2 . Consequently, the singular current density J∼1/s , making the total current finite, thus resolving the paradox. Furthermore, the strong packing of flux surfaces causes a steepening of the pressure gradient near the rational surface, with ∇p∼dp/ds∼1/s . In general non-axisymmetric MHD equilibrium, contrary to Grad’s conjecture that the pressure profile is flat around densely distributed rational surfaces, our result suggests a pressure profile that densely steepens around them.


Introduction
Magnetic field configurations with nested flux surfaces are desirable for plasma confinement fusion devices, due to the much faster motions of charged particles along the magnetic field lines than in transverse directions.Fusion devices with nested flux surfaces effectively insulate hot plasmas from device walls, facilitating good confinement.Axisymmetric fusion devices such as tokamaks can have a continuum of nested flux surfaces.However, for devices that are inherently threedimensional (3D), such as stellarators [1][2][3][4], the existence of a continuum of nested flux surfaces is not guaranteed.
Nevertheless, many existing magnetohydrodynamic (MHD) equilibrium solvers for stellarators, e.g., VMEC [5], NSTAB [6,7], and DESC [8], assume a continuum of nested flux surfaces as a point of departure.These solvers seek equilibria satisfying the MHD force balance equation either by minimizing the MHD potential energy (i.e., the sum of the magnetic and the thermal energies) or by directly imposing the force balance as a constraint.Here, p, B, and J = ∇ × B are standard notations for the plasma pressure, the magnetic field, and the current density.Toroidal 3D MHD equilibria with a continuum of nested flux surfaces can potentially give rise to current singularities at rational surfaces, where magnetic field lines form closed loops rather than fill the flux surfaces ergodically.To understand the issue, note that the force balance equation (1) implies that the component of the current density perpendicular to the magnetic field is given by By writing the condition ∇ • J = 0 implies and in general ∇ • J ⊥ ̸ = 0 if ∇p ̸ = 0.In a straight-field-line coordinate system (Ψ, θ, φ), the magnetic field is given by where Ψ is the toroidal flux function, θ is a poloidal angle, φ is a toroidal angle, and ι = ι (Ψ) is the rotational transform.In these coordinates, the magnetic differential operator B • ∇ is given by where J is the Jacobian of the coordinates.By using Fourier representations and the magnetic differential equation ( 4) can be written as Therefore, J ∥ /B can be expressed as where u 00 is a flux function and ∆ mn are coefficients to be determined [9,10].Equation ( 9) exhibits two types of singularities at rational surfaces, where the rotational transform ι are rational numbers.Let x = ι (Ψ) m − n.The first term in Eq. ( 9), known as the Pfirsch-Schlüter current, has a 1/x-type singularity at the rational surface with ι = n/m.The second term is a Dirac δ-function singularity, which is permitted because xδ(x) = 0.
The physical quantity associated with the current density is the total current through a surface.The Dirac δ-function type singularities are sheet currents at rational surfaces.Even though the current density becomes infinite, the total current is finite because the Dirac δ function is integrable.The 1/x-type singularities of the Pfirsch-Schlüter current pose a more serious problem.If the distances between the rational surface and neighboring flux surfaces are proportional to x, then the total current through a constant-φ surface element enclosed by θ, θ + dθ, x = ϵ 1 , and x = ϵ 2 is proportional to ´ϵ2 ϵ1 dx ′ /x ′ , which diverges logarithmically when ϵ 1 → 0. For this reason, the 1/x-type singularities are considered unphysical and should be avoided.
The 1/x-type singularities can be eliminated if the plasma pressure is globally flat, but that is not of interest to magnetic confinement fusion.Grad noted that magnetically confined MHD equilibrium solutions may be constructed by considering pressure profiles that are flat in a small neighborhood of each rational surface.However, if the rotational transform ι (ψ) is continuous, the magnetic differential equation ( 4) is densely singular and there will be infinitely many flatpressure regions throughout the entire plasma volume.Grad described such a pressure distribution as pathological [11].Along this line of thought, Hudson and Kraus constructed fractal-like pressure profiles with gradients localized to where the rotational transforms are sufficiently irrational, i.e., those that satisfy a Diophantine condition [12,13].Alternatively, one may allow the rotational transform ι to have discontinuities such that rational surfaces do not exist.
The central difficulty associated with the 1/x-type singularities is that they appear to give rise to divergent currents.However, this conclusion is obtained through heuristic arguments as we have discussed above rather than through a rigorous mathematical proof.The primary objective of this work is to reassess this conclusion through a simple model problem -the Hahm-Kulsrud-Taylor (HKT) problem [14,15].The HKT problem was originally posed as a forced reconnection problem, but in the ideal MHD limit, it also serves as a prototype for current singularity formation driven by resonant magnetic perturbations [16].The HKT problem is amenable to analytic solutions [17,18] and has been studied with various numerical codes including a Grad-Shafranov (GS) solver [19], a fully Lagrangian code [20], and the Stepped Pressure Equilibrium Code (SPEC) [21] for the case with p = 0 [16][17][18]22].In this study, we extend the established analytic and numerical solutions to include a non-vanishing pressure gradient.
This paper is organized as follows.Section 2 lays out the Grad-Shafranov formulation of the HKT problem and provides an asymptotic analytic solution.Section 3 compares the analytic solution with numerical solutions.We conclude and discuss future outlook in Section 4. y x ∇p Figure 1.The Hahm-Kulsrud-Taylor problem with a pressure gradient.The in-plane components of the magnetic field reverse directions at the mid-plane (the dashed line).The upper and lower boundaries are shaped by mirror-symmetric sinusoidal perturbations.In response to the perturbation, a δ-function singular current sheet develops at the mid-plane, which resonates with the boundary perturbations.The pressure gradient of the plasma drives the Pfirsch-Schlüter current density, which diverges algebraically near the resonant surface.The Hahm-Kulsrud-Taylor (HKT) problem, illustrated in Fig. 1, has a magnetized plasma enclosed by two conducting walls in slab geometry.Before the conducting walls are perturbed, the initial magnetic field is a smooth function of space in the domain with −a ≤ x ≤ a and 0 ≤ y ≤ L. The in-plane component points in the y direction with B y0 = x; together with a non-uniform B z0 component and the pressure p 0 , the system is in force balance; i.e., the equation is satisfied.Assuming a periodic boundary condition along the y direction, we then impose a sinusoidal perturbation with an up-down symmetry that displaces the conducting walls at x = ±a to x = ± (a + δ cos (ky)), where k = 2π/L is the wavenumber.As the system evolves to a new equilibrium subject to the ideal MHD constraint, current singularities, including a Dirac δ-function current sheet and a pressure-gradient-driven current singularity, will develop at and around the resonant surface x = 0 (the dashed line in Figure 1).In response to the boundary perturbation, the flux surfaces adjust their geometries while preserving the magnetic fluxes between them in accordance with the ideal MHD constraints.The new MHD equilibrium satisfies the Grad-Shafranov equation in Cartesian geometry where Here, the magnetic field in Cartesian geometry, with z being the direction of translational symmetry, is expressed in terms of the poloidal (in-plane) flux function ψ as In an equilibrium state, both the out-of-plane component B z and the plasma pressure p are functions of ψ.
Because the magnetic fluxes are conserved, the magnetic field is determined by the geometry of flux surfaces.To describe the geometry of flux surfaces, we first need to choose a variable to label them.Let x 0 be the positions of flux surfaces along the x direction before the boundary perturbation is imposed.The initial condition B y0 = x 0 yields the poloidal (in-plane) magnetic flux function ψ (x 0 ) = x 2 0 /2.While it is customary in toroidal confinement devices to label flux surfaces by the magnetic flux function, for the HKT problem it is convenient to label the flux surfaces with their initial positions x 0 because a single value of ψ can correspond to more than one flux surface.
With this labeling, the geometry of flux surfaces is described by a mapping from (x 0 , y) to (x, y) via a function x (x 0 , y).Using the chain rule, we can express the partial derivatives with respect to the Cartesian coordinates in terms of the partial derivatives to the coordinates (x 0 , y): Here, the subscripts of the partial derivatives on the left-hand side indicate the coordinates that are held fixed; the partial derivatives on the right-hand side are with respect to the (x 0 , y) coordinates.
Hereafter, partial derivatives are taken to be with respect to the (x 0 , y) coordinates by default, unless otherwise indicated by the subscripts.Using these relations for partial derivatives, the Cartesian components of the in-plane magnetic field are given by and The out-of-plane component B z is determined by the conservation of the out-of-plane magnetic flux as Here, B z0 is the initial z-component of the magnetic field and the flux surface average ⟨f ⟩ is defined as for an arbitrary function f (x 0 , y), where L is the period of the system along the y direction.
Likewise, the pressure is determined by the adiabatic energy equation as where p 0 (x 0 ) is the initial pressure profile and γ is the adiabatic index.
From the relation J = ∇ × B, the out-of-plane component of the current density is given by and the poloidal (in-plane) component J p is where B p = ẑ × ∇ψ is the poloidal (in-plane) component of the magnetic field.Substituting Eq. ( 21) in Eq. ( 11), the GS equation can be written as

Outer-Region Solution
For a small boundary perturbation, we may linearize the GS equation in terms of the displacements of the flux surfaces along the x direction ξ (x 0 , y) ≡ x (x 0 , y) − x 0 .To the leading order of ξ, the magnetic field components are and The linearized pressure is Using Eqs.(24-27) and the relation ψ = x 2 0 /2 in Eq. ( 23), the linearized GS equation now reads If we adopt the ansatz ξ = ξ(x 0 ) cos(ky), then ⟨∂ξ/∂x 0 ⟩ = 0 and the linearized GS equation reduces to The general solution of Eq. ( 29) is a linear superposition of two independent solutions and the boundary condition ξ(±a) = ±δ requires The independent solution cosh (kx 0 ) /x 0 in Eq. (30) diverges at x 0 = 0.If we insist that the linear solution is valid everywhere, then we must set c 2 = 0; the boundary condition (31) then determines the coefficient c 1 = aδ/ sinh (ka).However, this solution is not physically tenable.Because lim x0→0 sinh (k |x 0 |) /x 0 = k sgn(x 0 ), the geometry of the flux surfaces is approximately given by x ≃ x 0 + sgn(x 0 ) (kaδ/ sinh (ka)) cos (ky) in the vicinity of x 0 = 0. Hence, the flux surfaces with |x 0 | ≲ kaδ/ sinh (ka) overlap, causing a condition that is physically not permitted.We conclude that the linear solution is not valid within an inner region with |x 0 | ≲ O (kaδ/ sinh (ka)) and a nonlinear solution must be sought.

Inner-Layer Solution
A general nonlinear solution for the inner layer near a resonant surface was first developed by Rosenbluth, Dagazian, and Rutherford (hereafter RDR) for the bifurcated equilibrium after an ideal internal kink instability [23].The RDR solution was later adapted to the HKT problem with p = 0 [17,18].Here, we generalize previous approaches to incorporate the pressure-gradient effects.
Because the inner region is a thin layer, we assume the conditions |∂x/∂y| ≪ 1 and |∂/∂y|≪ |∂/∂x 0 | are satisfied.Under these assumptions, the dominant balance of the GS equation ( 23) is given by Integrating Eq. (32) yields where and g(y) is an arbitrary function that will be determined later by asymptotic matching of the innerlayer solution to the outer-region solution.The sgn (dψ/dx 0 ) factor comes from the requirement that ∂x/∂x 0 > 0 must be satisfied to avoid overlapping flux surfaces.Without loss of generality, we are free to set f (0) = 0, and the tangential discontinuity of B y at x 0 = 0 is then given by This tangential discontinuity corresponds to a Dirac δ-function singularity in the current density.
Using the flux function ψ = x 2 0 /2 for the HKT problem and integrating Eq. (33) again yields the inner-layer solution of RDR where the function h(y), which describes the geometry of the resonant surface, is yet to be determined.

Incompressibility Constraint
The functions f (x 0 ) and g (y) are not independent, but are implicitly related through Eq. ( 34) as we will see below.First, we determine B z and p in the inner layer by substituting x RDR (x 0 , y) for x (x 0 , y) in Eqs. ( 18) and (20).Then, applying the obtained B z and p in Eq. ( 34), the resulting relation gives a relation between f (x 0 ) and g (y).
We now show that Eq. ( 37) is approximately equivalent to the incompressible constraint in the strong-guide-field limit with B 2 p ∼ p ≪ B 2 z .From Eq. (37), the differentials df and dx 0 are related by where we have used the relation In the strong-guidefield limit, the first term in the right-hand side of Eq. ( 38) is the dominant term.As the leading order approximation, we may ignore other terms in the equation, resulting in the relation Integrating Eq. (39) yields where C is a constant.Finally, we apply Eq. (40) in Eq. ( 36) and obtain Equation (41) has a simple physical meaning: The plasma in the inner layer is uniformly compressed or expanded, and the constant C is the compression ratio.Because the inner-layer solution is to be matched to the outer-region solution and the latter is not compressive (up to the linear approximation), the compression ratio must assume the value C = 1; that is, the plasma approximately satisfies the incompressibility constraint.It should be noted, however, that the incompressibility constraint will not be valid if the boundary perturbation has an m = 0 Fourier component.In this case, the constant C needs to be determined by the asymptotic matching as well.

Asymptotic Matching
The inner-layer solution must be asymptotically matched with the outer-region solution to have a complete solution.For simplicity, we will assume the strong-guide-field limit and replace the relation (37) by the incompressibility constraint In this limit, to the leading order approximation, the effect of pressure gradient on the geometry of flux surfaces is negligible.Moreover, because the pressure gradient also does not affect the outerregion solution, the asymptotic matching procedure, which we will outline below, is identical to that for the case with p = 0 [17].
The outer-region solution is expressed in terms of the displacements ξ of flux surfaces.Therefore, to facilitate the matching, we express the flux surface displacement of the inner-layer solution as and examine its asymptotic behavior as x 0 → ±∞.
From the incompressibility constraint (42), we can infer that f → x 2 0 as x 0 → ±∞.Hence, in the asymptotic limit of x 0 → ±∞ we can split the integral in Eq. (43) into two parts, yielding This asymptotic behavior of the inner-layer solution as x 0 → ±∞ should match the behavior of the outer-region solution in the limit of x 0 → 0 ± , which is Matching Eqs. ( 44) and (45) yields h(y) = 0 and the following two conditions: and 1 2 Here, we use the notation ↔ in the second condition because, as we will see later, Eq. (47) can not be matched exactly.
The matching condition (46) gives an integral equation to determine the function g(y) as follows: by eliminating x ′ in favor of f using the incompressibility constraint (42) in Eq. ( 46), we obtain However, analytically solving this equation to obtain g(y) is a daunting task.By studying the behavior of g(y) around y = 0 and y = π/k, RDR suggested a function form g(y) ∝ sin 8 (ky/2) as an approximate solution.This suggestion was confirmed by Loizu and Helander with a numerical solution of the integral equation [24].Moreover, by fitting the function form with the numerical solution, they also determined the coefficient in front, yielding Next, for the matching condition (47), it is clear that g(y) cannot be matched with cos(ky) exactly.Following RDR, we expand g(y) as a Fourier series and only match the m = 1 term; that gives Finally, we use the relation (51) to eliminate c 2 in the boundary condition (31), yielding an equation for c 1 : The solution is Here, the sign for the square root in Eq. ( 53) is chosen such that c 1 ≃ δ/2 sinh(k/2) in the limit of a small perturbation.Now, we have obtained c 1 , c 2 , h(y), and g(y).We can then solve Eq. ( 42) to obtain f (x 0 ), which completes the necessary information to calculate the RDR solution, Eq. (36).Solving f (x 0 ) from Eq. (42) in general requires a numerical treatment.However, the leading order behavior of f (x 0 ) in the limit of |x 0 | → 0 can be obtained analytically as [18] f where the coefficient c f can be expressed in terms of the gamma function Γ [25] as

Pressure-Gradient-Driven Current Singularity
Our analysis thus far has shown that the effect of pressure gradient on the geometry of flux surfaces is negligible.We are now in a position to discuss the effect of pressure gradient on the current density distribution.We are primarily interested in the current density in the vicinity of the resonant surface, i.e., the inner-layer solution.
In the strong-guide-field limit, the plasma is approximately incompressible; therefore, from Eqs. ( 18) and ( 20) we have p (x 0 ) ≃ p 0 (x 0 ) and B z (x 0 ) ≃ B z0 (x 0 ) as the leading order approximation.To calculate the current density requires B ′ z (x 0 ), but we cannot simply assume B ′ z (x 0 ) ≃ B ′ z0 (x 0 ) for the following reason: when the guide field is strong, both B z (x 0 ) and B z0 (x 0 ) are close to constant, with different small variations to account for the force balance under different conditions, and these differing small variations give rise to different derivatives for the two fields.To obtain an approximation for B ′ z (x 0 ), we take the derivative of Eq. ( 34) with respect to x 0 , yielding Hence, Using Eq. ( 57) in Eq. ( 22) we can obtain the poloidal current density In the vicinity of x 0 = 0, f ∼ |x 0 | 8/3 and df /dx 0 ∼ |x 0 | 5/3 .Assuming p ′ 0 ̸ = 0, the pressure gradient term dominates and the poloidal current density diverges as J p ∼ p ′ 0 /x 0 .Next, we calculate the the out-of-plane current density J z = ẑ • ∇ × B p .The dominant component of B p is B y ≃ sgn (x 0 ) f (x 0 ) + g(y).Therefore, the out-of-plane current density which is not affected by the pressure.Using the asymptotic behavior (54), the leading order behavior of J z near x 0 = 0 is Hence, the out-of-plane component J z → 0 as x 0 → 0. Now, we put our results in the context of the conventional treatment of the Pfirsch-Schlüter current density.For a solution of the GS equation, the perpendicular component of the current density is A crucial point is that because of the formation of the Dirac δ-function current singularity, the poloidal magnetic field B p does not go to zero as at the resonant surface; instead, B p → ± g(y)ŷ as x 0 → ±0.Consequently, the magnetic differential operator B • ∇ is invertible, and the magnetic differential equation ( 4) implies that the parallel current density J ∥ diverges in the same way as ∇ • J ⊥ .‡ How do we reconcile the fact that the operator B • ∇ is invertible at the resonant surface with its non-invertible appearance in a straight-field-line coordinate system?To answer this question, it is instructive to construct such a coordinate system.Because the initial magnetic field lines are straight, the Cartesian coordinates are a straight-field-line coordinate system.Now, if we take a Lagrangian perspective by describing the final state in terms of the mapping from the initial positions x 0 = (x 0 , y 0 , z 0 ) of fluid elements to their final positions x = (x, y, z), the initial coordinates (x 0 , y 0 , z 0 ) are also a straight-field-line coordinate system for the final state due to the ideal MHD frozen-in constraint.In the Lagrangian perspective, the magnetic field at x is determined by the initial magnetic field B 0 at x 0 and the mapping x (x 0 ) via the relation [20,26] where J = det (∂x/∂x 0 ) is the Jacobian of the mapping.‡ Strictly speaking, it may be more appropriate to describe the operator B • ∇ as conditionally invertible at x 0 = ±0.
Because the function g goes to zero at y = 0 and L, the equation ± g(y)dϱ/dy = ς is solvable only when ς approaches zero sufficiently fast as y approaches 0 and L such that the function ς/ √ g is integrable.This condition is satisfied for −∇ • J ⊥ on the right-hand side of the magnetic differential equation (4).
Although our GS formulation is not fully Lagrangian, we can reconstruct the Lagrangian mapping of fluid elements from the initial to the final state once the solution is obtained.In this 2D problem, the fluid motion is limited to an x-y plane.So, the z coordinates of the initial and final states are identical, z = z 0 .In an x-y plane, for each fluid element labeled by (x 0 , y) in the final state, we need to find its initial position (x 0 , y 0 ).This "inverse" Lagrangian mapping can be expressed as a function y 0 (x 0 , y).From the conservation of magnetic flux through an infinitesimal fluid element and using Eq. ( 18) to relate B z0 and B z , we can calculate and integrate it along each constant-x 0 contour to obtain y 0 (x 0 , y).Using the RDR solution (36) in Eq. ( 70) yields In the limit of x 0 → 0, ∂y 0 /∂y → 0 everywhere except at y = 0 and y = L, where ∂y 0 /∂y diverges as ∼ |x 0 | −1/3 .Therefore, the straight-field-line coordinate system (x 0 , y 0 , z 0 ) becomes ill-behaved at the resonant surface, and that explains why the invertible operator B • ∇ appears to be noninvertible in such a coordinate system.Figure 2 shows a constant-z 0 slice of the straight-field-line coordinate (x 0 , y 0 , z ) for the HKT problem with L = 2π, a = 0.5, and δ = 0.1, where the solid lines are constant-x 0 or constant-y 0 coordinate curves.The singular nature of the coordinate system at the resonant surface is evident, as all the constant-y 0 coordinate curves converge towards y = 0 or y = L. Now, the critical question is: Does the 1/x-type current singularity lead to an infinite current?Fortunately, the answer is no, and the reason can be attributed to the Dirac δ-function singularity that simultaneously appears at the resonant surface.The finite tangential discontinuity B y x0=0 = 2 g(y), where • denotes the jump across an interface, arises from a continuous initial magnetic field through squeezing the space between flux surfaces [16].As we can infer from the RDR solution (36), for flux surfaces sufficiently close to the resonant surface such that the condition is satisfied, the mapping from the flux surface label x 0 to the physical distance s is quadratic: Because f (x 0 ) ∼ |x 0 | 8/3 and g(y) ≃ 4c 2 1 k 2 /3 sin 8 (ky/2), where k = 2π/L, the condition (72) eventually will be satisfied with a sufficiently small x 0 for all y except at y = 0 and y = L.As such, the mapping eventually will become quadratic, s ∼ x 2 0 , although the transition to the quadratic mapping will occur at different x 0 for different y.Because of the quadratic mapping, the poloidal current density J p ∼ 1/x 0 becomes J p ∼ 1/ √ s.And since ´s 0 ds ′ / √ s ′ is integrable, the total current does not diverge.
What about at y = 0 and y = L where the mapping is not quadratic?Because g = 0 at those locations and f ∼ |x 0 | 8/3 , Eqs. (58) gives J p ∼ |x 0 | 1/3 , which does not diverge.Moreover, to compensate for the strong squeezing of the quadratic mapping, the flux surfaces in the downstream regions near y = 0 and y = L have to bulge out.Consequently, the scaling at y = 0 and y = L becomes s ∼ |x 0 | 2/3 ; see Ref. [18].Intuitively, the J ∼ 1/ √ s scaling can be understood as follows.The quadratic mapping s ∼ x 2 0 causes a steepening of the pressure gradient: Because the pressure gradient is balanced by the J × B force for an MHD equilibrium, the perpendicular current density J ⊥ ∼ |∇p| ∼ 1/ √ s, and ∇ • J ⊥ diverges in the same manner.Finally, the parallel current density satisfies the magnetic differential equation ( 4), and the operator B • ∇ is invertible; therefore, J ∥ ∼ 1/ √ s as well.

Numerical Verification of the Analytic Solutions
We now compare the analytic solutions in Section 2 with numerical solutions using a GS solver.
The GS solver has been extensively tested, showing good agreement with the solutions of a fully Lagrangian solver [16] and the SPEC code [18] for the HKT problem with p = 0. Here, we consider an initial equilibrium with a linear pressure profile and a magnetic field The parameters in the following numerical calculations are: the domain sizes a = 1/2 and L = 2π; the guide field B 0 = 10; the mean pressure p0 = 1; the pressure gradient r = 1; the perturbation amplitude δ = 0.1; the adiabatic index γ = 5/3.In our previous studies, we assumed a mirror symmetry across the mid-plane and only solved the GS equation in half of the domain with x 0 ≥ 0. In this study, because of the asymmetric pressure profile across the resonant surface, we cannot impose the mirror symmetry.We divide the domain into two regions separated by the flux surface labeled by x 0 = 0 and solve the GS equation in each region by a descent method.Simultaneously with the iteration in both regions, the x 0 = 0 flux surface is allowed to move until the force-balance condition B 2 /2 + p = 0 is satisfied.
Panel (a) of Figure 3 shows a numerical solution of the B y component in the entire domain.Panels (b-d) show a few selected one-dimensional (1D) cuts of the numerical solution, together with the corresponding analytical solution (33).The analytic and the numerical solutions are in good agreement.The magnetic field is discontinuous across the resonant surface, corresponding to a Dirac δ-function current singularity.Figure 4 compares the numerical solution of the J z component with the theoretical prediction using Eq.(59).Here, J z → 0 at the resonant surface as the theory predicts.Likewise, Figure 5 shows comparisons between the numerical solutions of J y along selected 1D cuts and the theoretical predictions using Eq.(58).In this figure, the horizontal axis |x − x r | is the distance between the flux (b) (a)  Finally, Figure 6 shows a numerical solution of the pressure p in the entire domain and along selected 1D cuts.Here, the steepening of the pressure gradient due to the squeezing of flux surfaces near the resonant surface is evident.Examining the solutions in the vicinity of the resonant surface indicates that the ∂p/∂x ∼ |x − x r | −1/2 , as expected from our analysis.As a remark, although the analytic theory is derived under the strong-guide-field limit, we find that the theoretical predictions remain close to numerical results even for a moderate guide field such as B 0 = 2.

Discussions and Conclusion
In summary, we have derived an analytic theory for the pressure-gradient-driven current singularity near the resonant surface of the ideal Hahm-Kulsrud-Taylor (HKT) problem and have validated the theory with numerical solutions.Our key finding is that although the current density diverges as J ∼ 1/∆ι, where ∆ι is the difference of the rotational transform relative to the resonant surface, this singularity does not lead to a divergent total current.The reason is that the distance s between the resonant surface and an adjacent flux surface is not proportional to ∆ι.Because of the formation of a Dirac δ-function current sheet at the resonant surface, the neighboring flux surfaces are strongly packed, and the distance s ∼ (∆ι) 2 .Consequently, the current density J ∼ 1/ √ s, which is integrable and the total current is finite.
Our analysis also finds that the Pfirsh-Schlüter current density J ∥ and the diamagnetic current density J ⊥ both diverge as 1/ √ s, where the diamagnetic current density is the dominant component.The diamagnetic current density diverges because the strong packing of flux surfaces near the resonant surface causes a steepening of the pressure gradient; therefore, J ⊥ ∼ dp/ds ∼ 1/ √ s.Furthermore, solving the magnetic differential equation ( 4) yields the parallel current density J ∥ diverging in the same manner as the perpendicular component J ⊥ .Notably, contrary to the conventional wisdom that the pressure needs to flatten at rational surfaces to have an integrable current density, here we show that the current singularity remains integrable despite a steepening of the pressure gradient.
Our analysis assumes that the initial pressure gradient dp/dx 0 does not diverge and is nonvanishing at x 0 = 0.More generally, we may consider an initial condition with J p ∼ dp/dx 0 ∼ |x 0 | α , where α > −1 such that the current density is integrable.After the boundary perturbation, the current density J ∼ dp/ds ∼ s (α−1)/2 .Because the power index (α − 1)/2 > −1 if α > −1, the current density is integrable.Furthermore, we may consider an initial magnetic field B y0 ∼ sgn (x 0 ) |x 0 | ν .After the boundary perturbation, the formation of a Dirac δ-function current sheet requires the distance between flux surfaces to scale as s ∼ |x 0 | ν+1 ; consequently, the current density J ∼ dp/ds ∼ s (α−ν)/(ν+1) , which is integrable.In all these cases, as long as the initial current density is integrable, the current density in the perturbed state is also integrable.
Although our conclusion is obtained with a simple prototype problem, a similar conclusion probably applies to more general magnetic fields and potentially resolves the paradox of nonintegrability of the Pfirsch-Schlüter current density.If that is to be the case, then a weak (i.e., non-smooth) ideal MHD equilibrium solution with a continuum of nested flux surfaces, a continuous pressure distribution with non-vanishing pressure gradients on rational surfaces, and a continuous rotational transform is not prohibited.Such an equilibrium may be pathological, to quote Grad, but nonetheless mathematically intriguing to contemplate.

Figure 2 .
Figure2.A straight-field-line coordinate system for the HKT problem with L = 2π, a = 0.5, and δ = 0.1.The coordinate system becomes ill-behaved at the resonant surface.

Figure 3 .
Figure 3. Magnetic field component By (a) in the entire domain; (b) along x 0 = 0 + ; (c) along y = π; (d) along y = π/2.Black solid lines in Panel (a) are flux surfaces, and vertical dashed lines indicate the locations of the one-dimensional cuts of Panels (c) and (d).Black solid lines in Panels (b-d) are numerical results, and the red dotted lines are the predictions of the analytic theory.

Figure 4 .
Figure 4. (a) A numerical solution of Jz in the entire domain; (b) 1D cuts along y = π in Panel (a).The black solid line in Panel (b) is the numerical solution, and the red dotted line is the prediction of the analytic theory.

Figure 5 .
Figure 5. Comparisons between the numerical solutions of Jy and the theoretical predictions along a few 1D cuts.Here, the solid dots represent the numerical solutions and the dashed lines are theoretical predictions.Panel (a) shows the region with x 0 < 0 and Panel (b) shows the region with x 0 > 0. Here, the horizontal axis |x − xr| is the distance between the flux surface at x and the resonant surface at xr.As predicted by the theory, Jy diverges as |x − xr| −1/2 near the resonant surface.

Figure 6 .
Figure 6.(a) A numerical solution of the pressure p in the entire domain; (b) 1D cuts along the dashed lines in Panel (a).