Impact of a free normal velocity boundary on the modelling of external MHD modes

Abstract

Free normal-flow (NF) conditions at the plasma boundary are shown to be essential for three-dimensional magnetohydrodynamic (MHD) simulations to agree with linear stability theory. A comparative verification study is presented between two different formulations of the boundary conditions (BCs) for velocity perturbations: (i) fully consistent free NF implementation and (ii) fixed NF formulation, neglecting flow perturbations at the numerical boundary. Numerical results are compared with consolidated figures of merit from the linear theory of external kink modes. We consider two classes of initial equilibria presenting different numerical challenges: a uniform current channel surrounded by pure vacuum and a shaped Wesson-like equilibrium, with smooth (polynomial) radial dependency. Only the fully consistent free NF formulation is invariably accurate in modelling the plasma interface at the numerical boundary, without the need of enforcing a pseudovacuum region at the edge of the simulation domain, as in most analogous past studies. This study employs the cylindrical code SPECYL (Cappello & Biskamp, Nucl. Fusion, vol. 36, no. 5, 1996, p. 571) that solves a full-MHD model without pressure gradients, whose fully consistent resistive wall module with free NF BCs was recently successfully verified against the independent code PIXIE3D (Spinicci et al., AIP Adv., vol. 13, no. 9, 2023, p. 095111).

1 Introduction

Magnetically confined plasmas aimed at thermonuclear fusion best perform close to pressure- and current-driven magnetohydrodynamic (MHD) stability boundaries and are thus subject to a wide spectrum of potentially dangerous large-scale relaxation phenomena, e.g. the tearing modes, the vertical instabilities, the kink modes, the edge localised modes (see, e.g. White 2006; Wesson & Campbell 2011). The effect of such instabilities is generally deleterious (White 2006), unless stable nearby equilibria may be attained by nonlinear saturation or through externally imposed magnetic perturbations. This is the case in hybrid (Piovesan et al. 2017; Burckhart et al. 2023) and advanced (reversed shear) tokamak scenarios (Gormezano et al. 2008; Zhang et al. 2021), or in quasi-single-helicity states in the reversed field pinch (RFP) dynamo/flux-pumping effect (Bonfiglio, Cappello & Escande 2005; Marrelli et al. 2021).

Among the most dangerous MHD instabilities, an important class is represented by the free-boundary modes. Characterised by a finite fluid velocity at the plasma edge, they may lead to disruptions in the tokamak configuration (Granetz et al. 1996; Wesson & Campbell 2011); this is the case of vertical instabilities and external kink modes.

Advanced numerical modelling is crucial for understanding such complex and typically nonlinear dynamics. Starting from the 1980s, an increasing attention has been devoted to the formulation of realistic magnetic boundary conditions (BCs), mostly relying on resistive-wall modules (see, e.g. Gimblett 1986; Hender, Gimblett & Robinson 1989; Haines, Gimblett & Hastie 2013) typically implemented at the numerical boundary of MHD codes (as in Schnack et al. 1987; Strauss et al. 2004; Marx & Lütjens 2017; Bonotto et al. 2020; Bunkers & Sovinec 2020; Isernia et al. 2023), but not necessarily (Ferraro et al. 2016).

The problem of a free (fully consistent) normal velocity boundary has deserved renewed interest, also prompted by some theoretical studies (see Zakharov 2008; Zakharov, Galkin & Gerasimov 2012; Zakharov & Li 2015), suggesting that its absence in most present-day codes may prevent them from reliably modelling the vertical displacement events in tokamak plasmas. Two numerical studies have specifically investigated the impact of the velocity boundary formulation, performed, respectively, with the M3D code (Strauss 2014) and with the NIMROD code (Bunkers & Sovinec 2020), both demonstrating its critical implications for numerical predictions. Still more recently, realistic velocity BCs have been indicated among the unresolved demanding challenges for the modelling of tokamak vertical displacement events (see Artola et al. 2024).

The nonlinear code-to-code verification study of a resistive wall module with a fully consistent velocity boundary has been recently performed (Spinicci et al. 2023) between the SPECYL (Cappello & Biskamp 1996) and the PIXIE3D codes (Chacón 2008), extending a previous nonlinear verification study between the same two codes, where ideal wall BCs were adopted (see Bonfiglio, Chacón & Cappello 2010). The plasma interface was modelled as an axisymmetric thin shell (as in Hender et al. 1989), separated from an outer ideal wall by a vacuum region. The purely perpendicular velocity boundary featured a free normal-flow (NF) and allowed a fully self-consistent coupling within the MHD model. The quantitative reliability of both SPECYL's and PIXIE3D's NF boundary implementations was verified by a remarkably general mutual agreement on several test cases, both in the tokamak and the RFP configurations.

This paper employs the SPECYL code in a comparative study of the fully consistent BCs with a free NF, just outlined in the previous paragraph, against an approximated fixed NF formulation that does not allow finite velocity perturbations at the numerical boundary, often employed in nonlinear codes (analogous BCs are adopted, e.g. in Bonfiglio et al. (2010), McAdams (2014), Ferraro et al. (2016), Marx & Lütjens (2017) and Ramasamy et al. (2022)). The reliability of both is measured against well-known figures of merit from the ideal MHD linear theory of external kink modes (see, e.g. Freidberg 2014). The thin shell of SPECYL's BCs can be set to a high enough resistivity for the magnetic field diffusion through it to be almost instantaneous on the characteristic time scales of the plasma dynamics. In this condition, here dubbed the transparent wall limit, free-boundary modes driven by plasma current are expected to be unstable in the presence of a fully consistent velocity field. We show in this study that only SPECYL's free NF formulation allows proper modelling of the plasma interface through a transparent shell, in optimal agreement with the linear theory.

Remarkably, external kink modes simulations in the transparent-wall limit substantially differ from the typical approach adopted by other codes in the past, that leverages a high-resistivity and low-density pseudovacuum region (PV) to displace the numerical boundary (where BCs are enforced) away from the plasma interface. In the PV, high resistivity and vanishing density strongly damp the plasma flow, so that simplified (even not fully consistent) velocity BCs may be enforced with very low impact on simulations outcomes.

The paper is organised as follows. In § 2 we start with an overview of SPECYL's boundary module. Section 3 presents the numerical set-up, leveraged to pursue both the free NF velocity boundary approach and the simplified velocity boundary with PV. Section 4 contains the main body of our results: we study the external kink mode stability of two axisymmetric equilibria, characterised by a uniform current density and a shaped current density, respectively. Section 5 gives conclusions and outlooks. Four appendices are also present, dedicated to, respectively, an overview of linear kink stability analysis (Appendix A), the analytical profiles and growth rates for the simplified initial equilibrium characterised by a flat current (Appendix B), an overview of the semianalytical tool employed to produce the theoretical figures of merit for the Wesson-like current density profile (Appendix C) and a discussion in relation to the remarkable benchmark illustrated in Ferraro et al. (2016) (Appendix D).

2 The SPECYL code and its BCs

This section aims at a general overview of the SPECYL code and its boundary module, which is described in detail in Spinicci et al. (2023).

2.1 The SPECYL code

The SPECYL code computes the numerical solution of nonlinear three-dimensional (3-D) MHD equations by evolving in time $t$ the magnetic field $\boldsymbol {B}$ and the flow velocity field $\boldsymbol {v}$ according to the following viscoresistive scheme:

(2.1)\begin{gather} \rho(\partial_t\boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{v})=\boldsymbol{J\times B}+\rho\nu\nabla^2\boldsymbol{v} , \end{gather}

(2.2)\begin{gather}\partial_t\boldsymbol{B}={-}\boldsymbol{\nabla}\times\boldsymbol{E}, \end{gather}

(2.3)\begin{gather}\boldsymbol{E}=\eta\boldsymbol{J-v\times B}, \end{gather}

(2.4)\begin{gather}\boldsymbol{J}=\boldsymbol{\nabla}\times\boldsymbol{B}, \end{gather}

(2.5)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B}=0. \end{gather}

In the equations above, $\rho$ is the axisymmetric mass density, $\boldsymbol {J}$ is the current density, $\nu$ and $\eta$ are a uniform viscosity and an axisymmetric resistivity, respectively, and $\boldsymbol {E}$ is the electric field.

All quantities appear in dimensionless units: $\rho =\hat {\rho }/\rho _0$ (where $\rho _0=\hat {\rho }|_{r=0}$); $\boldsymbol {B}=\boldsymbol {\hat {B}}/B_0$ (where $B_0=|\boldsymbol {\hat {B}}|_{r=0{,t=0}}$); $\boldsymbol {v}=\boldsymbol {\hat {v}}/v_A$ ($v_A=B_0/\sqrt {\mu _0\rho _0}$ being the Alfvénic velocity, with $\mu _0$ vacuum permeability); $t=\hat {t}/\tau _A$ ($\tau _A=r_{\text {BC}}/v_A$ being the Alfvénic time and $r_{\text {BC}}$ the radial span of the computational domain, coincident with the plasma radius $a$ unless a PV region is included); $\boldsymbol {J}=\hat {\boldsymbol {J}}\mu _0 r_{\text {BC}}/B_0$; $\boldsymbol {E}=\hat {\boldsymbol {E}}/v_A B_0$. Hatted quantities are in physical units. In dimensionless units, the scalar kinematic viscosity corresponds to the inverse viscous Lundquist (or Reynold's) number, $\nu =\tau _A/\tau _\nu \equiv M^{-1}$ (where $\tau _\nu =r_{\text {BC}}^2/\nu$ is the viscous time scale); and resistivity corresponds to the inverse Lundquist number, $\eta =\tau _A/\tau _R\equiv S^{-1}$ (where $\tau _R=\mu _0 r_{\text {BC}}^2/\eta$ is the resistive-diffusion time scale).

The geometry consists of a periodic cylinder: a spectral approach is adopted for the periodic coordinates, $0\leq \theta \leq 2{\rm \pi}$ and $0\leq z\leq 2{\rm \pi} R$ ($R$ being the plasma major radius), while adopting a finite-differences staggered scheme in the radial direction, $0\leq r\leq r_{\text {BC}}$. The pressure gradients are neglected, and the axisymmetric density $\rho (r)$ and resistivity $\eta (r),$ along with the uniform scalar viscosity $\nu$ are kept constant in time. No pressure or temperature dependence is assumed for $\rho$ and $\eta$, thus presently limiting SPECYL's ability to accurately model highly nonlinear processes.

2.2 The boundary module

A thin cylindrical shell of arbitrary uniform resistivity is modelled at the plasma edge, separated from a coaxial ideal conductor by a tuneable-width vacuum region. In our notation, we call $r_{\text {BC}}\geq a$ the numerical boundary radius (always located at the thin shell position), and $b>r_{\text {BC}}$ the outer ideal wall radial position.

The magnetic boundary leverages the consolidated thin-shell relations (see e.g. Gimblett 1986), descending from the assumption that the radial-$\boldsymbol {B}$ is continuous across the wall, and prescribing

(2.6)\begin{gather} \partial_t B_{r,\text{BC}}=\frac{1}{\tau_w} [\partial_r (r B_r)]_{-}^{+}, \end{gather}

(2.7)\begin{gather}\boldsymbol{E}_w=\frac{r_{\text{BC}}}{\tau_w} [\boldsymbol{\hat{r}}\times\boldsymbol{B}_t]_-^+{+} E_{z}^{0,0}\hat{\boldsymbol{e}}_z, \end{gather}

where $\boldsymbol {B}_{r,\text {BC}}=\boldsymbol {B}_r\vert _{r=r_{\text {BC}}}$, $\tau _w=\mu _0 r_{\text {BC}} \delta _w/\eta _w$ is the resistive time scale of the magnetic field diffusion through the shell (with $\delta _w$ shell thickness and $\eta _w$ shell resistivity), $\boldsymbol {E}_w$ is the 3-D electric field inside the wall, $E_z^{0,0}$ is the uniform electric field representative of the inductive loop voltage and $\boldsymbol {B}_t=B_\theta \hat {\boldsymbol {e }}_\theta+B_z\hat {\boldsymbol {e}}_z$ is the tangential magnetic field, $\hat{\boldsymbol{e}}_\theta$ and $\hat{\boldsymbol{e}}_z$ being the azimuthal and axial unit vectors. Here, square brackets indicate the difference of the enclosed quantity between just outside the wall (in vacuum) and just inside (plasma edge); the vacuum solution is obtained in SPECYL from the analytical solution of Poisson's problem, and from the self-consistent value of $B_{r,\text {BC}}$, as thoroughly described in Spinicci et al. (2023).

Equation (2.7) provides a boundary condition for $\boldsymbol {B}_{t}$ by imposing the continuity of the tangential electric field between wall surface and plasma edge, via Ohm's law,

(2.8)\begin{equation} \boldsymbol{E}_w(\boldsymbol{B}_{t,\text{BC}})= \eta_\text{BC}(\boldsymbol{\nabla}\times\boldsymbol{B}_\text{BC})_t - (\boldsymbol{v}_\text{BC}\times\boldsymbol{B}_{\text{BC}})_t, \end{equation}

where the subscript ‘bc’ indicates that quantities are evaluated at numerical boundary.

This is a general set-up, capable of reproducing diverse experimental conditions, from an ideal wall attached to the plasma (when $\tau _w\gg \tau _{\text {sim}}$, $\tau _{\text {sim}}$ being the simulation duration), to a generic plasma-facing shell of arbitrary resistivity. In addition, the proximity of the outer ideal conductor can freely vary, from attached to the plasma to a virtually infinite distance.

If the magnetic diffusion through the resistive wall is instantaneous with respect to the time scales of plasma dynamics ($\tau _w\ll \tau _{\text {dyn}}$), the wall acts as if transparent to the magnetic field. In this regime, a fully self-consistent boundary module must be capable of reproducing free-boundary current-driven instabilities, including the external kink modes. Good quantitative accuracy is expected, as long as SPECYL's cylindrical geometry can realistically describe the evolution of a large-aspect ratio plasma, including the linear growth of the mode, possibly up to low-amplitude saturating kink modes (see Eriksson & Wahlberg 1997).

A 3-D velocity at the plasma edge, including a free NF, is mandatory for numerical self-consistency of resistive-wall boundary conditions. Moreover, allowing finite flow into the wall is required for a hot fusion plasma, which must be set free to impinge on the wall and annihilate (see Zakharov et al. 2012).

Two sets of velocity BCs hereinafter described are sketched in figure 1; a schematic view of the poloidal cross-section illustrates the plasma (pink), the thin resistive shell at $r=r_{\text {BC}}$ (claret) and the outer ideal wall (grey). For the case of SPECYL's free NF formulation (already described in Spinicci et al. (2023)), $r_{\text {BC}}=a$, and a dotted black contour suggests fully consistent modulation of the velocity boundary in reproducing an $m=2$ external kink mode. The fixed NF implementation enforces Dirichlet boundary conditions on velocity fluctuations and is mostly stable to the free-boundary modes, unless enforcing a suitable PV region (shaded in lighter pink) to displace the plasma interface away from the numerical boundary.

Figure 1. Schematic representation of the two boundary formulations. The plasma is represented in pink, surrounded by a thin resistive shell (claret) at $r=r_{\text {BC}}$ and by an outer ideal wall at $r=b>r_{\text {BC}}$ (in grey). In the free flow formulation (a), the plasma edge corresponds to the numerical boundary ($a=r_{\text {BC}}$), and a 3-D velocity is permitted: a black dotted contour suggests self-consistent edge flow modulation compliant with an $m=2$ external mode. Fixed flow boundary conditions are illustrated in (b): in this case, the numerical boundary is separated from the plasma interface by a PV region (in a lighter pink shade), and no flow perturbation is present at $r=r_\text {BC}$.

The free NF velocity boundary of SPECYL aims at an $\boldsymbol {E}_w\times \boldsymbol {B}$ drift at the plasma edge, thus purely perpendicular to the magnetic boundary,

(2.9)\begin{equation} \boldsymbol{v}_a=\frac{\boldsymbol{E}_w\times \boldsymbol{B}_a}{(\boldsymbol{B}_a)^2},\quad v_{{\parallel},a}=0. \end{equation}

A finite parallel velocity could be enforced by means of, e.g. a Bohm speed at the plasma boundary, even if its value would still be irrelevant to self-consistency with (2.6)–(2.8), as they only depend on $\boldsymbol {v}_{\perp }$.

For the purpose of this comparative study, we will also consider an approximated fixed normal flow formulation, enforcing Dirichlet boundary conditions for all non-axisymmetric fluctuations, on the top of a purely one-dimensional (radial and axisymmetric) pinch flow velocity (analogous to the one in Bonfiglio et al. (2010)), defined as

(2.10)\begin{equation} \boldsymbol{v}_\text{BC}^{0,0}=\frac{(\boldsymbol{E}_w\times \boldsymbol{B}_\text{BC}\boldsymbol{\cdot} \hat{\boldsymbol{e}}_r)^{0,0}}{(|\boldsymbol{B}_{t, \text{BC}}|^2)^{0,0}}\hat{\boldsymbol{e}}_r, \quad \boldsymbol{v}_\text{BC}^{m,n}=0\quad \text{for }m,n\neq 0, \end{equation}

where $\hat {\boldsymbol {e}}_r$ is the radial unit vector. Equation (2.10) describes the leading term of (2.9) in the standard RFP operation regime, while it is typically small in tokamak configuration; for the sake of this study, (2.10) represents almost a Dirichlet BC also for the axisymmetric harmonic of NF.

The theoretical figures of merit used throughout this work reproduce well-known results from the linear stability theory of external kink modes, and are derived leveraging the energy principle, briefly discussed in appendix A (see also Freidberg 2014). The resulting linear differential equation for the plasma displacement $\boldsymbol {\xi }$ is numerically solved, unless an analytical solution can be attained (as for the case reported in Appendix B), for any given initial axisymmetric equilibrium utilising the semianalytical LENS code (i.e. linear Euler and Newcomb solver, from the names of the two main differential equations). The LENS is implemented in the IDL (interactive data language) coding language and has been widely tested, both against analytical figures of merit and manufactured solutions (see Spinicci 2023). An overview of LENS is reported in Appendix C.

To comply with SPECYL's implementation, we will enforce the so-called straight tokamak limit theory (corresponding to a very large plasma aspect ratio), assuming a cylindrical plasma with a circular cross-section.

3 Numerical set-up

The external kink mode stability will be studied for the two velocity boundary formulations already reported in figure 1, and for two alternative initial equilibria, each defined by a specific initial current density distribution,

(3.1)\begin{equation} J_{z}^{0,0}(r) = \begin{cases} J_{z}^{0,0}(0)\, \left(1-\left(\dfrac{r}{a}\right)^2\right)^\lambda & \text{if }r\leq a \\ 0 & \text{if }r> a \end{cases}. \end{equation}

We will refer to the case of $\lambda =0$ as the flat current case, and to the $\lambda =1$ case as the Wesson current case, since this is the typical equilibrium employed in Wesson's theoretical works (see, e.g. Wesson 1978).

Throughout this study we will leverage an initial axisymmetric Ohmic equilibrium (as in Delzanno, Chacón & Finn 2008), defined as the solution of the force-free condition $\boldsymbol {J}\times \boldsymbol {B}=0$ and the parallel Ohm's law $\boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {B}=\eta \boldsymbol {J}\boldsymbol {\cdot }\boldsymbol {B}$. The use of these equations is meant to provide initial equilibrium profiles that are approximate solutions of the full system of equations solved by SPECYL, apart from the small convection and viscous terms in the momentum balance equation. The Ohmic equilibrium equations, as well as (2.1)–(2.5), assume mass-density and a resistivity only depending on $r$,

(3.2)\begin{gather} \eta(r) = \min{\{\eta_0 (1+A r^B)^C ;\ \eta_{\max}\}} , \end{gather}

(3.3)\begin{gather}\rho(r) = (1-\rho_\text{BC}) \left[\frac{1}{2}-\frac{1}{\rm \pi} \arctan{(\alpha(r-a))}\right] + \rho_\text{BC}. \end{gather}

The $\eta _{\max }$ in (3.2) is a truncation value preventing numerically demanding edge values of $\eta$. The two profiles in (3.2) and (3.3), along with the uniform scalar viscosity, are set as a specific input to each simulation. A uniform loop voltage in the plasma volume yields the initial equilibrium current density, roughly proportional to plasma conductivity $1/\eta$; the parameters $A, B$ and $C$ are chosen accordingly, to reproduce the desired initial current profiles as in (3.1). As for the plasma density, the model parameters $\alpha$ and $\rho _\text {BC}$ are chosen so to approximate a sharp-step-like profile across the plasma interface; for the case of the free flow formulation, where the interface is at numerical boundary, we adopt a uniform density profile ($\alpha =0,\ \rho _\text {BC}=1$).

Figure 2 illustrates the hereinafter considered Ohmic equilibria. Figure 2(a,b) correspond to the two alternative choices of BCs: the free and fixed NF formulations, respectively. Figures 2(a i) and 2(b i) correspond to the equilibria with a flat current density ($\lambda =0$), while figures 2(a ii) and 2(b ii) represent the Wesson-like initial current distribution ($\lambda =1$). For each of the four combinations, two diagrams are displayed: in the upper one the SPECYL's profiles of current density and safety factor (in black) are compared with their respective reference profiles (in pink and teal, respectively), whereas in the lower diagram the two corresponding input profiles of resistivity and mass density are also represented (solid lines), along with the theoretical sharp-step-like density profile (orange dashed).

Figure 2. Initial equilibria employed in this study. The panels correspond to different velocity boundary formulations: (a) free NF and (b) fixed NF with PV (shaded). Subpanels (i) and (ii) report equilibria corresponding to the two alternative initial current profiles enforced. For each of the four combinations two diagrams are reported: the superior one illustrates SPECYL's initial current density and safety factor (in black), in comparison with the analytical model (coloured dashed lines), whereas the lower one displays the two equilibrium-relevant input profiles of SPECYL's density and resistivity (solid lines). The theoretical (infinitely sharp) mass density profile is also reported as an orange colour dashed line. Unlike elsewhere, the radial coordinate on horizontal axes is normalised to the plasma radius $a$ for better graphical clarity.

As anticipated, the fully consistent free normal flow implementation allows full identification of the plasma boundary with the numerical boundary ($a=r_{\text {BC}}$), whereas the simplified fixed NF boundary requires that we displace the plasma interface away from the numerical boundary ($a< r_{\text {BC}}$) in order to get it linearly unstable: the resulting high-resistivity, low-density PV region is shaded.

All along, the displayed level of confidence in reproducing the theoretical reference profiles appears to be widely acceptable, as the only visible discrepancy – the smooth transition of the current density profile across $r=a$ for the flat current model with a fixed NF – is insufficient to produce sensible discrepancies, e.g. in the safety factor profile.

Table 1 presents the main simulation parameters, corresponding to the initial equilibria represented in figure 2, always enforced in the following unless otherwise specified.

Table 1. Numerical set-up for the initial equilibria represented in figure 2. For each one we report the coefficients defining the input profiles of plasma resistivity ($A,\ B,\ C, \eta _{\max }/\eta _0$) and of plasma density ($\alpha,\ \rho _\text {BC}/\rho _0$), the number of radial points ($N_r$), yielding comparable resolution in the plasma region, depending on its width ($a/r_{\text {BC}}$).

The edge current in Wesson's model, as much as both $J_z^{0,0}$ and $\rho$ in the PV region, should be null; this is numerically unfeasible in SPECYL. Instead, table 1 reports the most extreme values for $\eta _{\max }/\eta _0,\ \alpha \ \text {and}\ \rho _\text {BC}/\rho _0$ still compliant with numerical stability of the code time stepping scheme. To recreate more reliably the numerically challenging Wesson-like equilibrium, we evolve the axisymmetric initial state for a short time (approximately spanning the initial $1 \tau _A$) with magnified viscosity and resistivity ($\nu =\eta _0=10^{-2}$, $\eta _{\max }=5$) to substantially damp unwanted residual flow and current in the PV region, before starting the non-axisymmetric dynamics.

The plasma region width ($a/r_{\text {BC}}$) for the fixed flow case balances two opposing and quantitatively relevant requirements. On the one hand, as will be manifest in § 4, the Dirichlet condition for the velocity fluctuations, enforced by (2.10), produces a non-physical boundary layer in which the internal plasma flow unnaturally adapts to meet the imposed value in $r_{\text {BC}}$; this takes approximately the outermost $10\,\%$ of the numerical domain, and should be included in the PV region. On the other hand, since the plasma density and current in the PV cannot be exactly null, the PV width should remain significantly thinner than the plasma region width, to provide a negligible contribution to the integral mass and current (further amplified by the cylindrical Jacobian). Also, in equally spaced meshes like in SPECYL, a finite PV width proportionally enlarges the computational cost, by subtracting radial resolution from the plasma core region. For both initial equilibria, the optimal value of $a/r_{\text {BC}}$ had to be determined a posteriori, maximising the agreement with the theoretical expectations; in both cases, this coincides with $a=0.88 r_{\text {BC}}$, as reported in table 1. The number $N_r$ of radial mesh points is fixed so to have approximately the same resolution (${\sim }900$ points) for the plasma region with both fluid boundary formulations.

The dynamical evolution of such equilibria is here described with a two-dimensional (2-D) spectrum of modes, composed by the axisymmetric one, with the addition of several non-axisymmetric harmonics with a fixed helical pitch, such that $m=2n$. This configuration is adequate to study the linear stability to the external kink $m=2, n=1$, already considering just four non-axisymmetric modes with $m>0$, plus their conjugates. The initial destabilising contribution, driving the system towards relaxation, is applied in the form of a velocity perturbation,

(3.4)\begin{equation} \delta v_r^{2,1}=\frac{\mathcal{A}}{r} \sin{\left(\frac{{\rm \pi} r}{r_{\text{BC}}}\right)^4}, \end{equation}

where $\mathcal {A}=10^{-6}$. In all simulations we operate with a very fast resistive wall time scale, $\tau _w/\tau _A=10^{-6}$, in compliance with the transparent wall assumption, and enforce a confidently large aspect ratio $R/a=20$, as prescribed for the straight tokamak limit.

Uniform plasma viscosity, $\nu =M^{-1}$, and on-axis resistivity, $\eta _0=S_0^{-1}$, are kept mutually identical throughout this study (with the exception of Appendix D). Their specific value will be reported for each case study, individually. In all cases, the enforced values should tend to ideal plasma conditions (very low dissipation); this, along with the challenging equilibrium configurations already discussed and the inherently fast dynamics of external kink modes, requires very small time steps

Table 2 summarises the hitherto described input parameters underlying the dynamical evolution in our simulations. For some of them (namely, $S_0,~M\ \text {and}~b/a$) a range of enforced values is provided; for the ideal wall proximity, $b/a$, we have underlined the most commonly adopted value.

Table 2. Numerical set-up for the dynamical evolution of SPECYL's simulations in this study. We report the adopted values for the 2-D helical pitch ($m/n$), the number of modes with positive $m$ ($N_{m>0}$), the aspect ratio ($R/a$), the ideal wall proximity to the plasma edge ($b/a$), whose typical value is underlined, the transparent shell time scale ($\tau _w/\tau _A$), the on-axis Lundquist number and uniform viscous Lundquist number ($S_0$ and $M$) and the simulations time step ($\Delta t/\tau _A$).

4 Numerical results

We present in this section SPECYL's predictions of linear instability to the external kink mode 2/1 of the two equilibria introduced in § 3, and leveraging both sets of velocity BCs. It will be apparent that not only the free NF velocity boundary provides a better match for the analytical theory, but it also achieves more reliable asymptotic convergence to the expected outcomes.

4.1 The flat current model

The flat current model is widely present in numerical verification studies (e.g. in Hoelzl et al. 2014; Ferraro et al. 2016; Marx & Lütjens 2017), owing to its reduced complexity and to the existence of an exact analytical solution of the Energy principle (see Shafranov 1970; Wesson 1978). In fact, for a 2/1 external kink mode in the straight tokamak limit, the linear eigenfunction is $\boldsymbol {\xi }^{2,1}(r)=\xi _{r,a}^{2,1} (\hat {\boldsymbol {e}}_r+\text {i}\,\hat {\boldsymbol {e}}_\theta )r/a$ and the dispersion relation reads

(4.1)\begin{equation} \gamma^2\left(\frac{\tau_A}{\varepsilon}\right)^2 = \max{\left\{0;\ \frac{2}{q_a^2}(2-q_a) \left[1-\frac{2-q_a}{1-(a/b)^4}\right]\right\}}, \end{equation}

where the characteristic time scale is the toroidal Alfvénic time $\tau _A/\varepsilon =R/v_A$, $\varepsilon =a/R$ being the inverse aspect ratio. A complete analytical derivation is available in Appendix B. For the confidently large aspect ratio enforced and for $b/a\gg 1$, the characteristic dynamical time scale $\tau _{\text {dyn}}=1/\gamma$ predicted by (4.1) ranges in the tens of $\tau _A$; hence, our choice of $\tau _w/\tau _A=10^{-6}$ reliably fulfils the transparent-wall condition.

Figure 3 reports a convergence study to the theoretical growth rate from ideal MHD theory for a specific value of the uniform safety factor ($q(r)=q_a=1.5$) and considering various levels of plasma dissipation. As SPECYL plasma becomes increasingly ideal, we see that the fully consistent free NF boundary produces an asymptotic convergence to the expectation. A significantly different trend can instead be observed for the fixed flow boundary with PV, which non-monotonically oscillates around the target value for rather dissipative SPECYL plasmas, crossing it before $S_0=M=10^6$ and asymptotically plummeting towards stability as the plasma gets increasingly ideal (no visible mode growth can be detected in our simulations above $S_0=M=10^{12}$). Indeed, a very ideal plasma corresponds to a very ineffective PV, owing to the numerical limitations on $\eta _{\max }/\eta _0$ that keeps vacuum resistivity unphysically finite.

Figure 3. Comparison of linear growth rates for the flat current case study for $q_0=q_a=1.5$ and for various values of $S_0=M$. Free NF BCs asymptotically reach the ideal MHD theoretical value (dotted horizontal line) for vanishing plasma dissipation. The fixed NF BCs cross the target value around $S_0=M=10^{6}$ and decrease asymptotically towards stability (which is almost reached above $S_0=M=10^{12}$).

Figure 4 extends the growth rate analysis to a wider range of initial flat current equilibria, for various values of $q(a)$. The additional normalisation factor $q(a)/\varepsilon$ for $\gamma$ is widely present in the literature (e.g. in Wesson 1978) and it is therefore employed here as well. For the free NF formulation we employ $S_0=M=10^7$ (but similar results could be displayed for $10^8$ or higher), whereas for fixed flow BCs the dissipation level that results in the best agreement with linear theory has been chosen, corresponding to the target-crossing simulation (i.e. $S_0=M=10^6$; see figure 3). For the latter set of BCs, slightly more ideal simulations with $S_0=M=10^7$ are also shown in green, for comparison. Again, the theoretical dispersion relation (see (4.1)) is reported as a black line. Linear growth rates from SPECYL are numerically estimated with a logarithmic fit, except for $q(a)\leq 1$ and $q(a)\geq 2$ where a fast-oscillating low-amplitude signal is observed, lacking any detectable global trend.

Figure 4. Comparison of linear growth rates for the flat current case study, for various values of $q(a)$. The theoretical dispersion relation is represented as a black contour. (a) The free NF BCs (blue) stick robustly to the theoretical expectation; (b) the fixed NF BCs with PV are mostly consistent with the ideal MHD theory in the target-crossing set-up (red), besides slightly lower precision for $q(a)\leq 1.4$ and $\gamma <0$ in the ideally stable domain (due to sensible dissipation levels). Such an overall positive result is, however, frail, as choosing a slightly more ideal set-up (green) visibly spoils the previous agreement.

The simulations performed with a 3-D fluid boundary are in robust agreement with theoretical values (within $1\,\%$). The SPECYL simulations performed with a fixed velocity boundary and PV at the target-crossing dissipation level achieve a mostly comparable agreement with theory in the mode instability window with $1< q(a)<2$, despite somewhat lower precision on the left-hand edge of the hill. The effect of sensible plasma dissipation is, however, manifest in some simulations with $q\gtrsim 2$, significantly departing from the ideal mode stability condition $\gamma =0$ with a negative damping rate. Reducing the plasma dissipation level, the ideal stability condition is recovered, at the cost of a much weakened agreement with the theoretical expectation in the unstable region.

The radial velocity and magnetic field eigenfunctions for the three set-ups employed in figure 4 are benchmarked with the ideal MHD predictions, in figure 5, along with the radial profile of the axial current oscillations. The radial velocity is confronted with the radial component of the linear plasma displacement: $\xi _{r}^{2,1}/\xi _{r,*}^{2,1}=v_{r}^{2,1}/v_{r,*}^{2,1}$ (being $v_{r,*}=v_r\vert _{r=r_*}$, and $r_*=0.9a$ an arbitrarily chosen position for normalisation). Once again, the reliability of the free NF is remarkable. The same consideration applies to the fixed flow BCs for what concerns the plasma domain. Yet its behaviour in the PV is manifestly ill-behaved, in two main aspects: there is finite (and considerable) flow in the vacuum region, and the radial magnetic field pitch variation across the plasma edge reveals exaggerated surface currents. The strong deviation of SPECYL magnetic field in the PV from the vacuum linear eigenfunction (see (B9)) deteriorates with decreasing plasma dissipation and it is both sustained by the finite velocity, through a dynamo contribution, and by finite conductivity effects, allowing stray current fluctuations $\boldsymbol {J}^{2,1}=\boldsymbol {\nabla }\times \boldsymbol {B}^{2,1}$.

Figure 5. Benchmark of eigenfunctions profiles between SPECYL simulations with different velocity BCs and the analytical theory (black, dashed curves), for the flat current equilibrium. (i,ii) The SPECYL 2/1 radial velocity and magnetic field are compared with the radial displacement $\xi _r^{2,1}$ and the radial magnetic eigenfunction of ideal MHD: (a) the free NF formulation reproduces remarkably the theoretical curves; (b) the target-crossing simulation with the fixed NF BCs and (c) the slightly more ideal simulation with the same BCs are mostly reliable in the plasma region, with significant discrepancies in the PV. All profiles are normalised at $r_*=0.9a$, indicated by vertical dotted lines. (iii) We also report the (axial) current density oscillations for a unitary magnetic perturbation, compared with the analytical $J_z^{2,1}(r)\approx J_{\text {surf }}^{2,1}\delta (r-a)$: exceeding currents across plasma interface sustain exaggerated pitch variation of $B_r^{2,1}$ for the fixed flow formulation.

Figure 5(a/,iii,b/,iii,c/,iii) reports the dominant (axial) contribution to the current density fluctuations: for each subpanel, the linear current density is normalised to the amplitude of the magnetic perturbation. The analytical expectation for this case study is nearly proportional to a Dirac's delta distribution, centred across the plasma interface: $J_z^{2,1}(r)\approx J_{\text {surf }}^{2,1}\delta (r-a)$ (as discussed in Appendix B), which is quite well reproduced by the free NF implementation. The fixed NF formulation instead presents a much broader current distribution across plasma edge, much of which is on the PV side. Despite the current peak being lower than for the free NF case, the total surface current is larger, motivating the exaggerated pitch variation of $B_r^{2,1}$ across interface, as we already commented.

Up to this point the outer ideal wall in SPECYL's BCs formulations has always been kept largely distant from the plasma edge ($b/a=100$; see table 2). In figure 6 we explore the stabilising effect of bringing the ideal shell closer to the plasma edge, for two fixed flat current equilibria with $q(a)\approx 1.1$ and $q(a)\approx 1.5$ (qualitatively equivalent results would be obtained for any other value of $1< q(a)<2$). The set-up with dissipation values that mostly agree with the analytical theory ($S_0=10^{7}$ and $S_0=10^{6}$ for the free and fixed NF boundaries, respectively), are benchmarked with (4.1) (black line) showing remarkable qualitative agreement; until the ideal wall is as far as $b=3a$, its effect is mostly negligible and all numerical datasets correctly saturate to the far-wall limit growth rates already discussed for figures 3 and 4. As the ideal wall closes in, SPECYL's modes 2/1 progressively damp, up to the stability threshold proximity ($b/a\approx 1.78$ for $q_0=1.1$, $b/a\approx 1.19$ for $q_0=1.5$), which is accurately reproduced by both sets of BCs. As the ideal wall approaches further, the time traces of SPECYL's modes become fast oscillating noisy signals without any measurable global trend.

Figure 6. Scan on the outer ideal wall proximity $b/a$, within the flat current model and for two different equilibria ($q(0)=1.1$ in (a), $q(0)=1.5$ in (b)): ideal wall stabilisation of MHD activity becomes increasingly effective as it closes in to the plasma edge (shaded). The best performing set-ups for both the free and the fixed NF boundary are compared with the theoretical prescription (black). Two vertical dotted lines highlight the tightest possible proximity with either approach: owing to finite PV width, the fixed NF formulation cannot reproduce a closer proximity than $r_{\text {BC}}/a\approx 1.136$.

As for the quantitative agreement, the reliability of the free NF BCs lies within $6\,\%$ (ever decreasing towards $1\,\%$ as the ideal wall moves apart), whereas the level of compatibility of the fixed NF boundary lies within $20\,\%$ for $q_0=1.1$ and within $6\,\%$ for $q_0=1.5$. We highlight with two vertical dotted lines the tightest achievable proximity with either set of BCs, as $b=r_{\text {BC}}$; the presence of a PV region in the case of the fixed NF prevents simulating an ideal wall closer that $b\approx 1.136a$. For the present case study this inaccessible parametric region presents reduced interest, but would prevent simulating a portion of the unstable region for any flat current equilibrium having $q(a)\geq 1.6$, unless picking a larger, less favourable, $a/r_{\text {BC}}$ ratio, or using free NF conditions. Furthermore, this limitation could in principle restrict the ability of the fixed flow boundary to reproduce, e.g. real experimental conditions when a narrow vacuum interval separates the plasma edge from the first wall.

The displayed level of agreement for the fixed NF formulation in figure 6(a) is in line with what previously obtained with the JOREK-STARWALL code enforcing an analogous set-up for the same case study of a flat current, at $q(0)=1.1$ (see McAdams 2014). Substantially improved agreement with the analytical theory was achieved in Ferraro et al. (2016) with the M3D-C$^{1}$ code, by adopting a resistivity profile strongly discontinuous across the plasma interface; such results can be reproduced also in SPECYL, however, finding reduced reliability as the resonance radius approaches the plasma interface, as documented in Appendix D.

4.2 Wesson's model

Wesson's current model presents no known analytical solution for its eigenfunction $\boldsymbol {\xi }$ and eigenvalue $\gamma$. Alongside, its self-consistent modelling requires a numerically cumbersome fast-ramping resistivity profile at the plasma edge, corresponding to a vanishing current density at plasma interface, as already visible in figure 2. These are the reasons why it finds almost no application in verification studies of MHD codes.

Conversely, while certainly more realistic than a uniform current channel surrounded by pure vacuum, its numerical complexity makes it an important and challenging test for SPECYL's new BCs.

Also for this case study, the chosen numerical set-up enforces $\tau _w/\tau _A=10^{-6}\ll 1/\gamma$, and $R/a=20$. Theoretical expectations are obtained through the LENS code.

Figure 7 presents an analogous convergence study to the one presented in figure 3; the free and fixed NF BCs with PV (in blue and red, respectively) are compared with the expected theoretical growth rate (black dotted line), for a specific value of $q(a)$ and at various levels of plasma dissipation. The qualitative trends closely resemble what is already observed for the flat current equilibrium; only a 3-D and fully consistent free NF boundary can robustly reproduce the expected ideal MHD dynamical behaviour in the limit of low plasma dissipation, whereas the fixed NF presents just a (this time narrow) crossing of the target growth rate and becomes asymptotically stable as the PV becomes too inviscid and conductive to suppress MHD activity.

Figure 7. Convergence study to the theoretical growth rate predicted by the LENS code (black dotted line) of SPECYL's simulations with free and fixed NF boundaries, respectively, employing the Wesson-like equilibrium with $q(a)\approx 1.8$, for various levels of plasma dissipation.

The quantitative aspect is yet rather different from figure 3, since even the asymptotic convergence of the free NF is manifestly less effective, reaching a $17\,\%$ discrepancy only around $S_0=M=10^{10}$. The rationale for this is to be sought in the effect of strong and fast-ramping resistivity profiles at plasma edge. In fact, combining (2.3) and (2.4) from SPECYL's MHD model we get

(4.2)\begin{equation} \partial_t\boldsymbol{B}=\boldsymbol{\nabla}\times (\boldsymbol{v}\times\boldsymbol{B})- \frac{\text{d}\eta}{\text{d}r}\hat{\boldsymbol{e}}_r \times\boldsymbol{J}-\eta(\boldsymbol{\nabla}\times\boldsymbol{J}). \end{equation}

The first term on the right-hand side of (4.2) is the only one prescribed by the ideal MHD and corresponds to the dynamo action in the magnetic field variation. Of the other two terms, the one proportional to the resistivity gradient was vastly negligible in the plasma region in the previous case study with uniform $\eta$, as it still is now in the plasma core, where $\eta$ is mostly uniform. Conversely, at the plasma edge, both for the free and for the fixed NF BCs, this same term is now roughly one to two orders of magnitude larger than the other resistive term, for any assigned value of $S_0$.

The two velocity boundaries produce comparable predictions at high plasma dissipation. In fact, in such a regime the very high resistivity produces the effect of a PV region at the plasma edge even for the free NF case, strongly damping the edge velocity and thus making the two models quantitatively alike.

Looking at the fixed flow boundary simulations in figures 3 and 7 we notice a pretty similar qualitative trend, where numerical growth rates exceed the theoretical one at high dissipation and converge to stability for an ideal plasma. This is, however, not the case for the 3-D flow boundary simulations, that asymptotically converge to the expected growth rate in the ideal MHD limit, but they do it from below in figure 3 and from above in figure 7. This effect can be linked to the above mentioned tendency of large edge resistivity gradients (only present in Wesson's current modelling) to behave as PV. Notably, the PV has the function of displacing the effective plasma radius to a more internal position with respect to $r_{\text {BC}}$. The most direct consequence of this is that the effective edge safety factor must be evaluated at a more internal position and will thus be lower than the nominal value of $1.8$; as we will see in figure 8, $q(a)$ values slightly below $1.8$ are associated with larger growth rates according to the specific dispersion relation of Wesson's model.

Figure 8. The linear growth rates of several Wesson-like equilibria, corresponding to diverse values of $q(a)$, are illustrated. We compare the theoretical expectation obtained by the LENS code with (a) SPECYL's simulations with the free NF boundary (blue), at a very reduced plasma dissipation, and (b) SPECYL's simulations with fixed NF, corresponding to the target-crossing dissipation level (red), as discussed in figure 7, and for a mildly more ideal plasma set-up (green). In this challenging numerical test, the free NF formulation retains quantitative and qualitative reliability, unlike the fixed NF formulation.

In figure 8 we extend the verification benchmark to a range of initial Wesson-like equilibria, each defined by a different value of $q(a)$. The theoretical dispersion relation is obtained through the LENS code and is represented as a black contour. The least dissipative set-up of figure 7 is leveraged in SPECYL's simulations with free NF BCs (blue circles), whereas for the fixed NF BCs with PV the target-crossing dissipation level is represented in red, while a slightly more ideal set-up is reported in green. Also in this quite challenging case, the SPECYL code with its full 3-D velocity boundary displays a remarkable capability of reproducing the correct dispersion relation. Quite at the opposite, the fixed NF implementation seems drastically incapable of achieving even qualitative agreement with the theoretical expectations; not only are single growth rates strongly incorrect (by up to $190\,\%$ for $q(a)<2$ for the best performing set-up), but the instability window is completely inaccurate, predicting as significantly unstable equilibria that should instead be completely stable to the external kink modes. Internal resistive tearing modes for $q(a)>2$ cannot explain such behaviour, as they should present drastically lower growth rates for the given set-ups ($\gamma \tau _A/\varepsilon \sim 10^{-4}$, as estimated by LENS). Instead, the high resistivity at the plasma edge is once again blurring the interface between plasma and vacuum, as will appear evident from figure 9; the choice of a more conductive set-up (as represented by the green circles) importantly reduces the overstretching of the instability window, however, at the cost of further spoiling the agreement with theory for $q(a)<2$. Interestingly, SPECYL's numerical results with a fixed NF boundary are very sensitive to the specific value of $a/r_{\text {BC}}$: slight increase or decrease of such parameter significantly extend the unstable region beyond $q(a)=2$ and may even produce unphysical kink stability where $q(a)\lesssim 1.6$. In our perspective, such a strong dependence on a set-up parameter that retains no physical meaning would in itself already drastically undermine the reliability of fixed NF BCs in modelling free-boundary instabilities for the large class of equilibria with a vanishing edge plasma current.

Figure 9. Comparison of radial profiles of $v_r^{2,1}$ (i), $B_r^{2,1}$ (ii) and $J_z^{2,1}$ (iii), for free NF (in blue) and fixed NF BCs with PV (target-crossing case in red). A dashed black line represents the theoretical profiles of ${\xi }_r^{2,1}/{\xi }_{r,*}^{2,1}$, $B_r^{2,1}/B_{r,*}^{2,1}$ and $J_z^{2,1}/J_{z,*}^{2,1}$ ($r_*$ being marked by vertical dotted lines), as computed by the LENS solver. The large resistivity gradient at plasma edge produces deviations of SPECYL's eigenfunctions from the theoretical model, especially visible for the flow and the current density profiles. For the fixed NF case, the PV behaves as a higher-resistivity plasma rather than a vacuum: $v_r^{2,1}$ presents the vertical asymptote at resonance ($r/a\approx 1.054$, indicated by a cyan dotted line) typical of tearing modes, while a variation in the pitch of $B_r^{2,1}$ reveals finite surface currents $\boldsymbol {J}^{2,1}$ at the same resonance layer, also clearly visible in (b iii).

A more complete picture can be drawn by looking at the fluid and magnetic radial eigenfunctions in figure 9. Here, SPECYL's outcomes with the free and fixed NF boundary formulations, in blue and red, respectively, are compared with the radial displacement $\xi _r^{2,1}$ and the radial magnetic fluctuation $B_r^{2,1}$, as obtained from the LENS code. Despite the good agreement with theoretical radial profiles of both formulations within the plasma core region, high resistivity effects can be seen in all plots at the plasma edge, in the form of a mild reduction of edge radial velocity and magnetic field, correlated to an excess of edge current oscillations. It should be noted that the good performance of fixed NF is in this case, as well as in figure 5, strongly cherry-picked; we already commented on how frail the displayed agreement is upon small variations in physical conditions ($q(a),\ S,\ M$) or in $a/r_{\text {BC}}$.

Looking at the PV, shaded in figure 9(b), we may observe some apparent marks of resistive plasma dynamics in a region that we would rather regard as vacuum. In particular, $v_r^{2,1}$ presents a sharp change of sign across $r_0\approx 1.054a$, paired with a pitch variation of $B_r^{2,1}$ across the same radial position, revealing the presence of a (mostly axial, due to straight tokamak ordering) finite slab current $\boldsymbol {J}^{2,1}$, well visible in figure 9(b iii). This is the signature of a resistive internal (tearing-like) mode across a resonance layer: $q(r_0)\approx 2$.

5 Final remarks and future perspectives

New resistive wall BCs with a fully consistent 3-D velocity boundary were recently implemented in the nonlinear MHD codes SPECYL and PIXIE3D. Such a module consists of a thin shell of tuneable resistivity at plasma edge, separated from an outer ideal conductor by an arbitrary-width vacuum region. The magnetic boundary enforced at the plasma edge leverages the resistive Ohm's law across the thin shell and is self-consistently coupled with an $\boldsymbol {E}\times \boldsymbol {B}$ velocity field, dictated by the internal MHD and by Ohm's law at boundary.

In this paper we presented a challenging verification study against ideal MHD external kink modes theory. The conceptually simplest set-up to perform this study consists in modelling a free interface at plasma boundary, possibly also including the effect of an ideally conductive wall at finite distance. In practice, this is achieved by setting SPECYL's thin shell resistivity so high that it becomes transparent for the magnetic field dynamics; such a numerical set-up should be appropriate for the study of linear free-boundary MHD instabilities driven by plasma current.

With the fully consistent free NF formulation recently implemented in the SPECYL code, the linear evolution of external kink modes can be properly reproduced already looking at the transparent shell as a free interface in itself. This is, however, not the case with a not fully consistent set of velocity BCs, imposing null edge velocity perturbations. As the latter is mostly stable with a rigid wall attached to the plasma boundary, we leveraged a high-resistivity, low-density PV region to displace the plasma edge from the numerical boundary; such an approach is routinely employed in most analogous studies already present in the literature (see, e.g. Merkel & Strumberger 2015; Marx & Lütjens 2017).

The study considered the instability to external kink modes of two classes of initial tokamak equilibria characterised by circular cross-section and by a very large aspect ratio. The first one is defined by a uniform current channel, coincident with the plasma itself, and the second one features a Wesson-like parabolic current density that vanishes at plasma edge.

The first case study presents a well-known analytical solution for the linear instability eigenfunctions and growth rates; along with its reduced numerical complexity, this is why it is used in most numerical verification studies already present in the literature. This case study could be remarkably reproduced on a wide range of equilibrium current intensities and of external ideal wall positions by both sets of boundary conditions, each for a selected numerical set-up. However, the free NF formulation achieved robust and asymptotic convergence to the theoretical growth rates and eigenfunctions in the correct limit of an increasingly ideal plasma set-up, whereas the fixed NF formulation best fitting set-up was found for intermediate plasma viscoresistive values, while showing asymptotic stability in the ideal MHD limit. We identify the unavoidable numerical limitations in the modelling of the PV region as the most likely candidates to explain such behaviour.

The second case study presents no known analytical solution and needs enforcing a specific semianalytical solver; along with its increased numerical complexity, this makes its employment in this work rather exceptional, with respect to past analogous verification studies of nonlinear MHD codes. Also in this case, a robust and asymptotic convergence of SPECYL's free NF simulations was illustrated, even if reached at considerably lower dissipation levels than in the previous case, owing to increased numerical challenges. On the other hand, the fixed NF formulation proved itself extremely unreliable even in predicting bare modes stability. Growth rates analysis for a fixed equilibrium current intensity found the fixed NF formulation in SPECYL to be substantially incompatible with the expectation, both in the resistive plasma limit and in the ideal plasma limit (where it predicts once more full modes stability). Particularly concerning appears the strong and non-monotonic dependence of simulations outcomes on set-up parameters retaining no physical meaning, namely the distance of the numerical boundary from the plasma interface.

The main findings of this work are the following.

1. (i) A free NF formulation is essential for complete self-consistency of resistive wall BCs; indeed, it allows reliable modelling of the linear growth of current-driven external modes identifying the plasma interface with the numerical boundary. As far as we know, this result is unprecedented for a nonlinear MHD code.

2. (ii) The enforcement of a PV region at plasma edge is unavoidable in presence of BCs bearing a velocity formulation in partial consistency with the rest of the leveraged model. Nonetheless, its usage has general applicability in predictive MHD modelling, and it allows to describe the nonlinear evolution of external instabilities (especially in more advanced codes where $\rho$ and $\eta$ can evolve through time); its presence can indeed be justified to account for a scrape-off layer in realistic experimental conditions. Our study has, however, highlighted reasons for concern when its enforcement aims to account for an actual vacuum region, especially in the face of physical configurations that are challenging from the numerical point of view. In particular, the more or less narrow parametric region, where the theoretical models are matched with acceptable fidelity, appears to be hardly predictable a priori (unlike the asymptotic convergence to the right theoretical value, displayed by the fixed NF formulation). Even if a suitable set-up for fixed NF with PV formulations was eventually found in optimal agreement with Wesson's model predictions (possibly, extending the approach explored in Appendix D), an a posteriori determination of such a set-up might not be completely satisfying for a numerical tool that aims at performing predictive modelling.

3. (iii) Albeit the specific set-up featuring a transparent wall, coincident with the plasma interface, has reduced applicability in predictive simulations, we hope to provide with it a remarkably challenging and repeatable test for the self-consistency of free NF BCs in MHD nonlinear codes. The SPECYL's documented reliability on this critical set-up is expected to hold in the general case of a finitely conductive wall, or even considering a cold and sparse scrape-off layer around the plasma core.

Along with the already published verification study against PIXIE3D, this paper completes the numerical verification process of SPECYL's resistive wall module. The code will now be employed in both first-principle studies and predictive numerical modelling in the tokamak and the RFP configurations. For the RFP, preliminary numerical studies performed with the SPECYL code enforcing a simplified formulation of its boundary module with a fixed NF velocity like in (2.10) are reported in Marrelli et al. (2021). Analogous modelling work is now being carried out with the free NF boundary formulation, to extend the validity of previously achieved results and in view of anticipating the next device RFX-mod2 (Marrelli et al. 2019; Terranova et al. 2024). A central point of interest for the RFP community concerns the self-organised helical plasma states at high current regimes, with periodicity $m=1,\ n=7$ in the RFX-mod device (see Bonfiglio et al. 2013; Marrelli et al. 2021), which still lack self-consistent numerical predictability. Stimulated by the promising results presented in this paper we are confident in SPECYL's resistive wall module future findings also in this merit.

Free NF BCs are expected to play a determinant role especially in plasma-limiter configurations, namely including RFP routine operation and tokamak disruptions. The impact of a fully consistent fluid-magnetic formulation could be of especial interest in the modelling of externally imposed resonant magnetic perturbations at the plasma edge, widely used in numerical simulations for the modelling of disparate physical conditions, e.g. including edge localised modes triggering (Cathey et al. 2020), hybrid tokamak scenarios (Piovesan et al. 2017) and stimulated QSH states in RFP configuration (Veranda et al. 2017).

First-principle studies are also being pursued for the numerical validation of Zakharov's theory on vertical displacement events in the tokamak configuration. Such a model predicts the existence of external kink-like instabilities flowing into the wall and sustaining surface currents in its passive structures; the combination of these currents with the plasma-confining magnetic field would be at the origin of experimentally observed intense horizontal loads on device support structures (see Zakharov et al. 2012). Conceptually similar simulations to the ones presented in this work, but considering more realistic wall resistivities, will prove if a fully consistent velocity boundary with a normal velocity perpendicular to the plasma edge is indeed sufficient to produce the so-called wall-touching kink modes and sustain quantitatively comparable currents with respect to the model predictions.