The Island Equilibrium and Transport (IslET) model is extended to allow for the effects of ion compressibility albeit for a constant electron temperature. The model describes fully relaxed, saturated, steadystate magnetic islands arising from the evolution of a tearing mode caused either by an instability or an externally imposed resonant perturbation. It assumes that the collision frequency lies in the semicollisional regime. A salient new feature of the extended model is an equation describing the transport of parallel momentum. The model is applied to evaluate the drag exerted on a thin island by the emission of driftacoustic waves.
I. INTRODUCTION
Magnetic islands are created by tearing modes through magnetic reconnection in the vicinity of field lines that close upon themselves.^{1,2} Due to their detrimental effect on confinement, they are of concern in a broad range of magnetic plasma confinement devices but particularly in tokamaks and stellarators. In confinement devices with auxiliary heating, islands are primarily observed as the result of either neoclassical tearing modes (NTM) or error fields. Avoidance, control, and suppression methods have been developed, in particular, using electron cyclotron current drive (ECCD). Due to the significant power required for suppression as well as concerns that control systems could be overwhelmed by multiple tearing modes, there is persistent interest in improving the understanding and prediction of magnetic island evolution.
Fitzpatrick and Waelbroeck^{7,8} solved the problem of determining the pair of functions $ V f ( W )$ and $ D ( W , V )$ in a twofluid plasma by separating the problem into two regimes based on the relative magnitude of the parallel phase velocity $ \omega * e / k  $, where $ \omega * e = k Y c s \rho s / L n$ is the electron diamagnetic frequency and $ k  $ and $ k Y$ are the wave number components in the parallel and Y directions, respectively. Here, $ c s = T e / m i$ is the ion sound speed, $ \rho s = c s / \omega c i$ is the ionsound Larmor radius, $ \omega c i$ is the ion cyclotron frequency, and $ L n = n /  \u2207 n $ is the density scale length. The two regimes where the problem was solved are as follows:

The hypersonic regime $ \omega * e / k   \u226b c s$ corresponding to thin islands, for which the width of the island is narrow compared to the wavelength of the driftacoustic wave that it excites;

The subsonic regime, corresponding to wide islands, for which the island is wide compared to the radial wavelength of the driftacoustic wave it excites.
The solutions provided in these papers resulted from the numerical integration of a reduced system of nonlinear ordinary differential equations. The papers described scaling of the solutions with the principal parameters of interest. The asymptotic reduction of the problem to ordinary differential equations, however, required the adoption of simplifying assumptions:

In the hypersonic regime, the island propagation frequency was assumed to be close to the electron diamagnetic frequency;

In the subsonic regime, ad hoc profiles were adopted to connect the solutions inside and outside the separatrix.
Whereas these assumptions and their consequences are in agreement with presently available numerical solutions of the complete initialvalue problem, it is desirable to expand their applicability by seeking solutions free of these assumptions. This is of interest, for example, in the presence of gradients in the electron temperature when the natural frequency of the island may differ from the diamagnetic frequency. The present paper is a first step in this direction based on the use of an alternative reduced model for island evolution, the IslET model,^{9–11} that is based entirely on an asymptotic reduction of the “four field” equations describing the evolution of the island.
The transport analysis embodied in the IslET model^{9–11} calculates the two functions $ D ( W , V )$ and $ F ( W , V )$ for steady state islands, that is, for islands such that W and V are constant in time, including electron temperature gradient effects. Equivalently, it calculates the external drive parameter $ \Delta \u2032$ that would cause the tearing mode to saturate when the island width reaches W and the propagation speed V if the island is not subjected to external forces. The model has contributed to the understanding of island evolution by clarifying the effects of island “magnetization” (describing the magnitude of the island width compared to the ionsound Larmor radius $ \rho s$) and density flattening, accounting for the changes in the island propagation velocity from the unmagnetized to the magnetized regime as well as the role of coupling to the drift wave^{10} and the effect of the electron temperature gradient on the island stability.^{11}
A limitation of previous versions of the IslET model, which inhibits comparison to Refs. 7 and 8, is its neglect of the ion parallel velocity. This results in an incompressible description of the ion fluid and suppression of the ionacoustic or sound wave. In the linear regime, coupling to the sound wave is known to have a stabilizing effect on the tearing mode.^{12–14} In the nonlinear regime, likewise, numerical solutions of the compressible equations^{15–17} as well as analytic solutions in limiting cases^{16,18,19} show that the sound wave can have important effects on island evolution. In particular, it is known to be responsible for the flattening of the density inside the island and the accompanying reduction in the island propagation velocity when its width W is such that $ k   c s \u226b \omega * e$^{15,17} (i.e., $ W \u226a \rho s$). The reduction in the relative velocity between the island and the surrounding ions has also a stabilizing effect through the polarization current.^{16–18} Furthermore, the analysis of island evolution in the large island limit shows that compressibility is necessary in order to account for the effect of curvature.
Compressible models for the island evolution have been presented and solved numerically^{16,17} and semianalytically^{16,19} by Fitzpatrick and collaborators. Additionally, previous studies have addressed the interaction of sheared flows with magnetic islands, particularly the deformations caused by viscosity of the magnetic flux contours under sheared flow conditions.^{20,21} Specifically, Ren et al. showed in Ref. 20 that a deformation of the $\psi $ contours occurs due to the viscous sheared flow and calculated its effect on the evolution of the magnetic islands. Smolyakov et al. studied in Ref. 21 the effect that the deformations of the $\psi $ contours have on the island stability. However, the present work is not focused on these effects. In the present paper, we follow the work in Ref. 11 by extending the IslET model to include ion compressibility in order to enable a higherfidelity description of island evolution that can shed light on unresolved questions, such as the cause of the bifurcation in the island dynamics observed in Ref. 22 and modeled in Ref. 7. In order to focus in the ion compressibility effects, we do not consider electron temperature gradients. Therefore, we are extending the IslET model by including compressible effects for constant electron temperature.
The remainder of the paper is organized as follows. Section II presents the dynamical system whose steadystate solutions the IslET model aims to describe. Section III describes the derivation of the compressible IslET model, solving the lowestorder limit of the dynamical model to determine the equilibrium solutions in terms of profile functions and the firstorder equations to obtain transport equations governing the profiles. In Sec. IV, we consider the solution of the IslET model in the thin (unmagnetized) island limit. In Sec. V, we discuss our results and describe the work that is needed to complete the analytic description of island evolution in an inhomogeneous plasma.
II. FORMULATION
A. The drift model
B. Boundary conditions
Solving the problem gives us the force as a function of $ v edge$, $ F y ( v edge ) = \mu lim x \u2192 \u221e U$. It is noteworthy that the calculation of this force holds significance independently, as realworld islands consistently encounter drag forces arising from their interaction with the resistive wall of the confinement device, error fields, and internal electromagnetic perturbations.
In the island frame of reference, the natural propagation velocity of the plasma is given by the root of the forcebalance condition $ F y ( v root ) = 0$. Conversely, in the frame where the background plasma is at rest, the natural propagation velocity of the island is obtained by the opposite of the root of the forcebalance equation, $ v free = \u2212 v root$.
C. Transport ordering
Under the adopted ordering the solutions depend on three ratios of transport parameters $ C / \mu \u22a5$, the Schmidt number $ Sc = \mu \u22a5 / D$ and $ C / D$. This transport ordering facilitates the reduction of the system (3)–(6) to a set of equilibrium (Sec. III C) and transport (Sec. III D) equations in the steadystate limit $ \u2202 / \u2202 t = 0$.
III. DERIVATION OF THE COMPRESSIBLE ISLET TRANSPORT MODEL
A. Introduction
The described transport ordering enables a partial solution of the governing equations (3)–(6) in the steadystate scenario, wherein all time derivatives vanish in the reference frame where the island is stationary. This approach reduces the equations to a set of two equilibrium equations and a set of four 1D transport equations describing the profiles. Subsequently, the IslET (Island Equilibrium and Transport) model is formulated based on these resulting equations. This model was previously developed in Refs. 9 and 11 incorporating electron temperature gradients but under the assumption of incompressibility ( $ v = 0$). Here, we extend the derivation to incorporate parallel ion compressibility, characterized by $ \u2207   v \u223c v E \xb7 \u2207 n$. However, for the sake of focusing on compressible effects, we omit consideration of electron temperature gradients.
B. Mathematical preliminaries
The firstorder field equations take the form of differential equations along the magnetic field (magnetic differential equations) or stream lines (convective differential equations).
C. Lowest order: Equilibrium
To lowestorder ( $ O ( \delta 0 )$), we neglect all the transport terms and set the time derivatives to zero so that the equations of motion describe an equilibrium state to this order.
1. Electron equilibrium
2. Ion equilibrium
3. Summary of equilibrium results
Together, the four profile functions representing the isothermal electron streamfunction $ H ( \psi )$, the potential vorticity $ K ( \phi )$, the canonical azimuthal momentum $ G ( \phi )$, and the inductive current $ I ( \psi )$ completely specify the variation of the four fields n, $\phi $, $\psi $, and v across the island.
D. Firstorder: Electron transport
In this section, we derive transport equations for the profiles by expressing the solubility conditions for the firstorder corrections. Since fluid equations for the loworder moments involve higher moments, we start with the highestorder moment equations and work our way down to the lowestorder moment, the continuity equation, for each particle species, starting with electrons.
1. Electron momentum transport
2. Electron particle transport
E. Firstorder: Ion transport
1. Ion momentum transport
2. Ion particle transport
F. Boundary conditions
Equations (41), (43), (51), and (54), complemented by the boundary conditions (8) and (9), constitute the complete nonlinear IslET model that determines the size and velocity of magnetic islands. In view of the reduced nature of the IslET model compared to the complete dynamical model, however, it is necessary to revisit the question of boundary conditions.
G. Overview of solution
Here, we regroup the transport equations of the IslET model before discussing possible solution methods.
1. Governing equations
2. Constant $\psi $ regime
Up to this point, our analysis has refrained from making any assumptions about the flux function $ \psi ( x , y )$. To further simplify the equilibrium and transport equations, we may use the constant $ \psi \u0303$ approximation, where $ \psi \u0303$ is the amplitude of the flux perturbation. As we can see from $ \u2207 2 \psi = 1 + O ( \beta \u0302 )$, the constant $\psi $ approximation is justified whenever $ \beta \u0302 \u226a 1$. With the abovementioned approximations, the solution is determined by two equilibrium parameters: the island halfwidth $ w = 2 \psi \u0303 1 / 2$ and $ \Delta \u0302 \u2032 = \Delta \u2032 \rho s / \beta \u0302$, as well as by the three ratios of transport coefficients that control the profiles near the island.
IV. ANALYTIC SOLUTION FOR THIN ISLANDS
We have developed an analytical solution of the compressible equilibrium and transport model for unmagnetized islands. This solution neglects electron temperature gradients and considers two distinct regions: the inner region, where $ x \u223c w \u226a 1$, corresponding to the magnetic island itself, and the outer region, where $ w \u226a x \u223c 1$. We also obtain the slip curves for the velocity of the island in the inner region obtained from matching to the lowest order solution in the outer region.
A. Thinisland regime: Outer region solution
Therefore, in the outer region, we have two equilibrium equations: Eq. (68) for the lowest order to find $ \Phi 0 ( x )$ and Eq. (69) for the firstorder to obtain $ \Phi 1 ( x )$. To solve each equilibrium equation, we need $ K ( \phi ) , H ( \psi ) , G ( \phi )$ profiles. Following our zerothorder convention, we denote $ K 0 , H 0 , G 0$ as the zerothorder solutions.
The constants $ \delta H \u221e$, $ \delta G \u221e out$, $ \delta G \u221e out \u2032$, and $ \delta K \u221e$ are integration constants determined by matching with the inner region solutions (Subsection IV C). Moreover, once we establish the relation between $ \delta G \u221e out$ and $ \delta G \u221e in$ in terms of $ \delta G \u221e$, they can be determined by the asymptotic behavior of the parallel velocity: $ \delta G \u221e = v \u221e$.
In the incompressible limit $ \epsilon c = 0$, we recover the isothermal limit of the result found for the outer region in Ref. 11. The appearance of the $ \epsilon c 2$ term in the abovementioned equation shows that compressible effects appear in the electric potential to first order.
B. Thinisland regime: Inner region solution
C. Solution matching
D. Summary of results for the profiles
The final expressions of the profiles are
E. Slip curves
When the force exerted on the island changes, it modifies its velocity. Plotting the velocity as a function of the force gives the slip curves. Varying the compressibility parameter $ \epsilon c 2$ produces different slip curves.
Figure 1 shows the slip curves giving the island velocity $ v 0$ as function of the force $ F y / \mu $, comparing the exact results obtained from the numerical evaluation of the integrals $ h a$ and $ h c$ (solid lines) with the approximations of Eqs. (105) and (106) (dotted lines). Each color corresponds to a different value of $ \epsilon c 2$. It can be seen that, at constant force, an increase in the velocity values is observed as $ \epsilon c 2$ increases. Bifurcations appear at low magnitudes of $  F y / \mu $ for large enough $ \epsilon c$. The convergence point, where all lines converge, corresponds to $ F y / \mu $ when $ v 0 = 0$, indicating the absence of electric drift velocity.
Moreover, the ratio between the normal resistivity C and the particle diffusion coefficient D can be changed; when $ D \u226b C$, the velocity increases much faster than when $ C = D$ and bifurcations appear for smaller values of $ \epsilon c 2$; when $ C / D = 1$ (left side of Fig. 1), bifurcation appears for $ \epsilon c 2$ near 7, but when $ C / D = 3$ (right side), it appears for $ \epsilon c 2 = 3$.
Equilibrium solutions exist when $ F y = 0$, corresponding to the crossing point on the velocity axis (yaxis) in Fig. 1. This is not, however, the natural velocity of the island, as it would be the case if the complete solution to first order had been used. It is just the free rotation that the island would have if there was not a back reaction from the outer region due to the absorption of the onacoustic wave by the plasma. For $ C / D = 1$, there is one equilibrium solution for every $ \epsilon c 2 < 7$; however, when $ \epsilon c 2 = 7$ there are two equilibrium solutions, and a third one appears for every $ \epsilon c 2 > 7$. This behavior is also present in the right side, $ C / D = 3$, but for $ \epsilon c 2 = 3$ instead of 7.
The free rotation velocity, $ v free$, of the magnetic island can be expressed as a function of the ion compressibility parameter $ \epsilon c 2$ by imposing the constraint $ F y = 0$ in Eq. (113).
Plotting the free rotation velocity $ v free$ as a function of the compressibility factor is equivalent to making the bifurcation diagram of the slip curves.^{22} In Fig. 2, these equilibrium points are plotted for $ C / D = 1 , \u2009 2$ and 3. For each $ C / D$ scenario, there are three solution branches; one in the electron diamagnetic direction, $ v free > 0$, and two in the ion diamagnetic direction, $ v free < 0$.
The electron branch corresponds to a branch of sources, indicating unstable points where solutions tend to diverge away from. Additionally, in the ion direction, a saddlenode bifurcation occurs when $ \epsilon c 2$ exceeds a threshold value corresponding to the node point of the respective $ C / D$ case. Consequently, the node splits into an upper branch and a lower branch. The upper branch represents sinks, meaning that stable equilibrium solutions can be reached when $ \epsilon c 2$ is greater than the threshold value for small values of the velocity in the ion direction. However, the velocity decreases and tends toward zero as $ \epsilon c 2$ increases.
On the other hand, the lower branch constitutes a branch of sources, signifying unstable solutions. With increasing $ \epsilon c 2$, the velocity grows in the negative direction and tends toward infinity in the ion direction.
These results were obtained through $ \delta H \u221e$ resulting from the matching of the inner region solution with the lowestorder approximation of $\phi $ in the outer region. This means that the drift acoustic waves radiated by the island do not play any direct role in the inner region, there is no momentum redistribution by the acoustic waves in the inner region for thin islands at lowestorder; they indirectly play a role through the boundary conditions. As a result, the behavior observed in Figs. 1 and 2 is not the expected one according to previous numerical computations,^{7,17} where it was found that the island velocity decreases when $ \epsilon c 2$ is increased. Therefore, it is necessary to go to the next order in the solution for the outer region.
V. DISCUSSION
We have extended the Island Equilibrium and Transport model in Ref. 11 to describe a compressible plasma but restricting to constant temperature, focusing on describing the forces exerted by the plasma on the magnetic island, particularly the region around the magnetic island.
Previous works studied this model by considering small values of the ion compressibility parameter,^{7} found only numerical solutions,^{17} or adopted physically motivated assumptions for profile functions H and G.^{18} Here, we developed an analytical solution for arbitrary $ \epsilon c 2$ and $ v 0$, by calculating the profiles through solving the transport equations in the thinisland limit.
Under our framework, the inner region solution found for $ x \u223c w$ has to be matched to the outer region solution ( $ x \u226b w$) in order to get the complete solution. This is accomplished by an expansion in the island width w and here the matching was done to the lowest order outer region solution. Although matching to higher order is needed for a full solution, the force computed to this order provides the torque exerted by the inner region which is transmitted to the outer region by radiation of driftacoustic waves, which is what we were interested in. We found that this drag force $ F y$ depends on two integrals; one proportional to the particle diffusion coefficient D and $ 1 \u2212 v 0$, and the second one proportional to the velocity $ v 0$ and the normalized resistivity $ C = 0.51 ( \nu e / \omega * e ) ( m e / m i ) / \epsilon c 2$ times the ion compressibility parameter, i.e., the second term in the force is proportional to $ 0.51 \u2009 v 0 ( \nu e / \omega * e ) ( m e / m i )$. When the measure of the collisionality, C,^{23} is larger than the particle diffusion D, the second term in $ F y$ dominates and the velocity becomes more susceptible to the applied force, the velocity increases faster. At constant applied force, the velocity of the island increases when the ion compressibility increases either in the electron direction $ v \u221e > 0$ or in the ion direction $ v \u221e < 0$, and when $ \epsilon c 2$ exceeds a threshold value a saddlenode bifurcation appears.
We studied the slip curves bifurcation diagram, $ v free$, of the magnetic island as a function of the ion compressibility parameter and obtained that there are three equilibrium solution branches; an unstable branch in the electron direction, two branches coming from a bifurcation in the ion direction. The bifurcation node corresponds to the threshold value for $ \epsilon c 2$. Thus, we have three cases, when $ \epsilon c 2$ is less than the threshold value, the velocity increases in the electron direction; when $ \epsilon c 2$ exceeds the threshold value, the upper branch in the ion direction is stable and the velocity approaches to zero as $ \epsilon c 2$ increases, and the lower branch that is an unstable solutions branch which tends to infinity in the ion direction as $ \epsilon c 2$ increases.
In future work, we will solve the electric potential in the exterior region including momentum change. This will complete the description of thin island evolution. Note that the equation for the electrostatic potential in the exterior region is linear and was already studied in the context of the investigation of mode penetration;^{24} it has also been obtained with a semianalytical, iterative approach in Ref. 7.
In the present work, we have removed the necessity for an iterative approach by solving the nonlinear problem in the inner region analytically.
ACKNOWLEDGMENTS
This work was funded by U.S. DoE Contract No. DEFG0396ER54346 and by projects PAPIITUNAM IN110021 and CONACyT A1S24157, Mexico.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
M. S. CancinoEscobar: Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). J. J. Martinell: Formal analysis (equal); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). F. L. Waelbroeck: Conceptualization (lead); Formal analysis (equal); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
Data sharing is not applicable to this article as no new data were created or analyzed in this study.