1-D Modeling of the Screw-Pinch Plasma in PROTO-SPHERA

A simple steady-state model for a 3-species mixture (ions, electrons, and neutrals) in a screw-pinch plasma configuration is developed. The model is applied to the central plasma column of the PROTO-SPHERA experiment. Degree of ionization, azimuthal current density, and azimuthal ion velocity are calculated. Full ionization is found at plasma temperatures above 1.5 eV, with neutrals confined in an outer shell where radial plasma flow develops and drives both azimuthal current and azimuthal flow.


Introduction
Magnetically confined plasmas are commonly formed inside a toroidal vacuum vessel surrounded by a toroidal magnet.Plasma-facing components and conductors are subject to severe thermal, nuclear, and electromechanical loads at the inboard side of the torus.To give a proof of principle of a magnetic confinement configuration which can overcome this problem, the PROTO-SPHERA project [1,2] is aimed at producing hot toroidal plasma in a simply connected machine topology, i.e., without solid elements between the plasma and the symmetry axis.The toroidal magnetic field is generated by current flowing in an axial plasma pinch (called the centerpost plasma discharge).The main project task is to form and sustain a very low aspect ratio (i.e., nearly spherical) toroidal plasma encircling the centerpost discharge (Figure 1).Current in the centerpost discharge has been driven so far up to 10 kA, both in Argon and in Hydrogen.The target current, which will be driven after an upgrade of the power supply, is 60 kA.The centerpost discharge is shaped by an externally applied magnetic field to wet large-area annular electrodes.The voltage difference between electrodes is regulated to drive the preset plasma current, which heats the plasma and generates the toroidal (azimuthal) magnetic field.The circuit is closed by return conductors placed outside the vacuum vessel.From fast-camera images, the plasma appears as a linear discharge with variable diameter for a length of about 1 m and takes the shape of a mushroom in front of both annular electrodes (Figure 2), in agreement with magnetic field lines geometry.The plasma region is surrounded by a larger volume filled by neutral gas.The linear part of the discharge is essentially a screw pinch with variable cross-section.The plasma diameter is about 0.5 m at the equatorial plane and decreases to about 0.13 m at the ends of the linear part, consistent with clearance left by two tight-fitting magnetic coils; plasma luminosity sharply drops at field lines that intercept protections of such coils.The main properties of the centerpost plasma are the capability of conducting the desired current within the voltage limits of the power supply and the uniform spreading of current on the anode surface.Such avoidance of anode arc-anchoring is likely due to spontaneous plasma rotation.The line-average density has been measured at the equatorial plane by a second-harmonic laser interferometer [3].Typical line-average number densities are 4 × 10 20 m −3 and 1.5 × 10 20 m −3 in Ar and H discharges, respectively.Temperatures in the 2-8 eV range have been found in preliminary Langmuir probe measurements in the mushroom-shaped region.Electric potential measurements are available at the electrodes and at the magnetic coil casings.The voltage drop between the coils at the ends of the linear section gives a rough estimate of the axial electric field, a typical value being 50 V/m.Plasma current and electric potentials are in steady state for most of the plasma duration, which, depending on power supply settings, ranges from 0.3 to 0.8 s.
A modeling activity is being carried on to help devising optimization strategies and to forecast plasma properties at higher centerpost current.To keep the problem as simple as possible at this preliminary stage, geometrical complexity is neglected and the centerpost plasma is modeled assuming cylindrical symmetry.This approach allows capturing the essential features of plasma dynamics by means of simple calculations and provides useful insights for setting up more realistic models.The purpose of the analysis presented in this work is to evaluate qualitatively (i) the degree of ionization, (ii) the plasma rotation and (iii) the self-generated azimuthal current.The latter could provide a seed for torus formation.
Density measurements and preliminary temperature measurements indicate that mean free paths of plasma particles are much smaller than the plasma size.A fluid description with three interacting species (ions, electrons, and neutral gas) [4,5] has then be adopted.Neutral gas is referred-to as neutrals in the following.The 1-D analysis for ions and electrons is similar to the one developed for magnetized gas discharge plasmas [6,7], but uniform neutral gas density was assumed in those works, while neutral dynamics is kept into account in this work.The effects of decreasing neutral density inside the plasma (neutrals depletion) were reviewed in [8], but without including the confining force that arises from plasma current.On the other hand, screw-pinch analyses including the full electromagnetic force, see for example [9], usually assume fully ionized plasma.
General governing equations are collected in the Appendix A. Simplifying assumptions on which the pinch model is based and the ensuing simple transport equations are presented in Section 2. Numerical results are presented and discussed in Section 3. Conclusions and perspectives are given in Section 4.
Figure 2. Image of PROTO-SPHERA plasma in Argon as reconstructed from three cameras, one for the anodic region (top), one for the nearly cylindrical region (middle) and one for the cathode (bottom).View at transitions from linear to mushroom-shaped regions is impeded by coils.The mushroom-shaped plasma is clearly visible in the anodic region at the top.The cathode region appears obscure since sensitivity has been adapted to hot electron emitters (visible as an array of brilliant spots).

The Cylindrical Pinch Model
The model for the centerpost plasma is based on some simplifying assumptions which lead to a description in terms of ordinary differential equations.First, the plasma configuration is represented as a cylindrical screw pinch, with negligible flows (excepting axial current) in the axial direction and with constant quantities in both axial and azimuthal directions.This description seems reasonable in the region between the main constrictions shown in Figure 1; no attempt is done to extend the analysis to the mushroom-shaped regions.Second, steady-state condition is assumed, consistent with the experimental observation of flat-top conditions for most of the plasma duration (0.3-0.8 s, with 5 ms initial transient).Third, charge separation is neglected; charged plasma sheaths at plasma-solid interfaces are out of the scope of this work.The model cylindrical plasma is assumed to carry a uniform current density within a radius (R), which mimics the outer radius of the bright region in Figure 2, which comprises field lines that do not intersect coil protections.The present analysis does not apply beyond such radius, where the plasma flow becomes of the scrape-off type, with substantial axial losses towards solid surfaces.Equations for mass continuity and momentum balance in a mixture of electron, ion and neutral fluids [4,5] are employed.The Argon case is considered, assuming singly charged ions.Energy balance is left out of the present analysis: uniform temperature is assumed and its effect is studied parametrically; modeling energy balance in Ar is a formidable task as it requires accounting for radiation losses from many excited states.
General equations are collected in the Appendix A, while relevant components in cylindrical coordinates (r, φ, z) are considered in the following subsections.The mixture is described by ion velocity V i , plasma density n, electric current density j, plasma temperature T (assumed common to electrons and ions), neutrals velocity V N , neutrals density N and neutrals temperature T N .The z components of velocities are eliminated by the assumption of negligible axial flows.The radial component of current density is eliminated by assumed steady state and 1-D geometry.The other two components are eliminated using Ohm's law.Plasma temperature is regarded as a parameter and neutrals are assumed to remain at room temperature.Azimuthal components of velocities are eliminated using the respective momentum balances.The radial velocity of neutrals is eliminated by the continuity equation.Three non-linear differential equations for radial ion velocity, plasma density, and neutrals density remain, as detailed in the following subsections.Equations are expressed in SI units for all quantities excepting temperature, which is measured in electronvolts.Final equations are also given in dimensionless form.
The main limitation of this 1-D model is the implication of zero radial current density.The presence of this component, which is possible in 2-D provided that z-dependence and axial return currents preserve charge neutrality, could modify the momentum balance and then it could affect results on azimuthal components of both current density and velocity.On the other hand, 1-D results on the degree of ionization should be exportable to more realistic configurations.

Mass Continuity
Continuity Equations (A1) and (A2) in steady-state and cylindrical symmetry reduce to 1 r where the plasma number density n = n e ≈ n i has been introduced and recombination has been neglected.The ionization rate coefficient in m 3 /s as a function of electron temperature in eV units is [10] where i is the first ionization energy ( i =15.76 eV for Ar) and V Te is the electron thermal speed Equations ( 1) and ( 2) imply that regular solutions fulfill

Azimuthal Momentum Balance
Model assumptions on cylindrical symmetry and steady state imply radial magnetic field B r = 0, radial current density j r = 0, azimuthal electric field E φ = 0 and dp/dφ = 0. Retaining only dominant rate coefficients, i.e., neglecting k ion when summed with k eN or k iN collision coefficients and neglecting m e k eN when summed with m i k iN , the azimuthal component of Ohm's law takes the form: where ω ce = eB z /m e is the electron gyrofrequency and the electron-ion collision rate coefficient e is related to e-i collision frequency ν e = nk ei and to Spitzer resistivity η S = m e k ei /e 2 .Constant values of k eN and k iN coefficients are assumed, as listed in Table 1.The azimuthal plasma momentum balance reads Substituting ( 6) into ( 5) and neglecting m e k eN when summed with m i k iN one obtains If the electron cyclotron frequency is dominant, i.e., then the azimuthal current density takes the simple form The azimuthal ion velocity is estimated assuming zero net azimuthal momentum and neglecting convective inertia in (6): A posteriori numerical checks based on this expression show that ( 7) is widely verified, while neglect of convective inertia in passing from (6) to ( 9) is only marginally justified.

Radial Momentum Balance
The radial component of plasma momentum balance reads where the total plasma pressure p = p i + p e = ne(T i + T e ) (temperatures being in eV units) has been introduced and small coefficients have been neglected as in the previous subsection.The j z B φ term is mainly due to pinch force from the externally imposed plasma current; the plasma pressure it can balance is given by the Bennett relation or its kinetic generalization [11].The j φ B z term arises from the self-generated azimuthal current.The radial component of neutral gas momentum balance reads

Transport Equations
Equations presented in previous subsections could be closed with energy balance equations [5] accounting for heat transport and radiation.A simpler approach is followed in this work: equations are closed assuming a common constant temperature T e = T i as a parameter and the character of solutions is studied for different values of the parameter.Neutrals are assumed to remain at room temperature T N .
From the axial component of Ohm's law, uniform electron temperature implies uniform Spitzer resistivity and then nearly uniform j z , in fact the the last two terms in (12) are negligible, as verified in numerical calculations.The pinch term is then easily calculated: where I p is total axial current and R is the pinch outer radius.Centrifugal terms and radial variations of B z are neglected, as justified a posteriori.Combining ( 1), ( 4), ( 8), ( 10), ( 11) and (13) a system of three ordinary differential equations results for plasma density, ions radial velocity, and neutrals density.
where f p = µ 0 I 2 p /(2π 2 R 4 m i ) and plasma and neutral gas squared sound speeds have been introduced, c 2 s = e(T e + T i )/m i and c 2 N = eT N /m i respectively.A convenient dimensionless form can be obtained introducing the normalized radius x = r/R, the Mach number M = V ir /c s and normalized densities g = n/n B and G = N/n B , where is the radial excursion of density in the absence of radial flow (alike to a Bennett pinch). ( where p 1 = ω ce ω ci R/(c s n B k ei ) is related to slowing of perpendicular diffusion by magnetization, s is a normalized ionization rate, and p 5 = c 2 N /c 2 s .Zeros at the left-hand sides of (17) and (18) correspond to the well-known sonic singularity, which is usually met at plasma-solid interfaces, where it is resolved keeping into account substantial charge separation.
Another important characteristic of the system is introduced by neutrals dynamics.In fact, plasma velocity has a source term proportional to neutrals density, namely the 4th at the RHS of (17).On the other hand, neutrals density can steeply increase with radius since, according to (1), the radial plasma flux nV ir (or gM) is an increasing function of radius, and then the second term on the RHS of (19) can give rise to super-exponential radial growth.If this happens, the source term in (17) can bring radial velocity to hit the sonic limit.
In practice, since physically acceptable solutions must have no singularities inside the simulation domain , this effect sets an upper limit to the central density of neutrals.The limit changes dramatically with plasma temperature, as discussed in Section 3, the main reason being the strong temperature dependence of the ionization rate coefficient, see Equation (3).
The system is solved in the r ≤ R domain, with boundary conditions at the plasma center V ir (0) = 0, n(0) = n 0 and N(0) = N 0 .The central plasma density n 0 is adjusted to fit the measured line-average density (excepting particular cases as detailed in Section 3.1).The central neutrals density N 0 is adjusted to have N(R) values close to the surrounding gas density, as discussed in next section.

Numerical Results
Equations ( 14)-( 16) are solved numerically employing a standard python ode solver.To avoid the geometric singularity at r = 0, solutions are advanced starting from a slightly off-axis boundary r 0 = 10 −3 , with boundary conditions V ir (r 0 ) = 0, n(r 0 ) = n 0 and N(r 0 ) = N 0 .The number density of the surrounding gas is N out ≈ 10 21 m −3 ; the boundary condition on central gas density N 0 is adjusted to obtain N(R) slightly below N out .The boundary condition on central plasma density n 0 is chosen to reproduce the measured line-average density whenever possible (see next subsection).System parameters are temperature (through k ion , k ei and c s dependencies), plasma radius (being the pinch force coefficient f p ∝ R −4 ), plasma current (being f p ∝ I 2 p ) and axial magnetic field (being ω ce ω ci ∝ B 2 z ).Plasma current is fixed at I p = 10 kA.Different values of temperature, plasma radius, and plasma density are considered in the following.Values of the parameters for each case, including the dimensionless ones used in Equations ( 17)-( 19) are listed in Table 2.

Dependence on Temperature
The plasma radius varies along the pinch axis from 0.27 m to 0.065 and the external magnetic field varies from 0.02 T to 0.6 T; intermediate values of R = 0.12 m and B z = 0.17 T are assumed in this subsection.The character of solutions changes dramatically with temperature; temperature values are chosen in the range 1.0-2.5 eV, which is adequate to capture the different possible regimes.
For T e = 1.0 eV, the minimum possible n 0 , corresponding to g(0) ≈ 1 has been assumed.The resulting neutrals density varies by less than 20 % across the plasma radius, while the radial profile of plasma density has nearly parabolic shape (Figure 3).The ion radial velocity tends to diverge towards the edge while its value is only 4 × 10 −3 times the sonic limit.The reason is the small n at denominator in the first term at the RHS of Equation ( 14), which arises from finite ion radial flux nV ir required in steady state to drain fresh ions produced by ionization.The line-average density is 4.7 × 10 20 m −3 , i.e., about 20% above the measured value.Lower densities give rise to physically unacceptable solutions as, at this low temperature, plasma pressure becomes insufficient to balance the pinch force and then density decreases to zero inside the simulation domain.The estimated axial electric field E z = ηj z is 235 V/m in this case, much higher than the experimental value of about 50 V/m.
For T e = 1.5 eV, there is sufficient room to decrease n 0 and match the experimental line-average density; the chosen value n 0 = 5.5 × 10 20 m −3 corresponds to g(0) = 1.2.At this temperature, neutrals density varies by almost two orders of magnitude from the plasma center to the edge, while the plasma density profile remains parabolic with low edge value (Figure 4).The line-average density fulfills the experimental constraint in this case.Radial velocity tends to diverge towards the edge; the reason is the same as in the 1.0 eV case.The axial electric field is 128 V/m, still substantially higher than the experimental estimate.
For T e = 2.0 eV the experimental density is matched with n 0 = 5 × 10 20 m −3 , which corresponds to g(0) = 1.5.The neutrals density is very small across most of the plasma, and steeply increases in an outer shell (Figure 5).The plasma density profile is flatter than a parabola and the edge value becomes significant.The axial electric field is 83 V/m.Radial velocity steeply increases towards the edge; the reason in this case is not a small denominator but it is the cumulative effect of the source term proportional to neutrals density, as discussed at the end of Section 2.4.Very low N 0 is required in this case because the radial growth of neutrals density is steep at this temperature, and with larger N 0 the sonic singularity would occur in the simulation domain.The slope of N(r) keeps increasing with radius, so the profile cannot smoothly match the outer gas density N out ; this matching must occur in the scrape-off layer, where field lines intercept coils protections and the radial plasma flow nV ir decreases because plasma escapes in the axial direction and recombines at solid surfaces.The scrape-off layer is outside the scope of the present analysis, as discussed above.Profiles at higher temperature are qualitatively similar to the 2.0 eV case, but with flatter plasma density and steeper increase of neutrals density in the outer shell.Results at T e = 2.5 eV are shown in Figure 6.The estimated axial electric field is 60 V/m in this case, in reasonable agreement with the experimental estimate.The dramatic change of profiles at T e > 1.5 eV (far below the 15.76 eV first ionization energy of Ar) is due to reduced neutrals penetration depth.Neutrals penetration is pushed by external gas pressure and contrasted by collisions with outgoing ion flux, which in turn depends on ionization.
The normalized neutrals penetration depth δ can be estimated from ( 1) and (19).In dimensionless variables the former reads 1 x Keeping only terms that according to numerical results, are dominant at the plasma edge, in particular the second term at the RHS of (19), one obtains 1 The depth decreases with temperature as k ion strongly increases.At T e = 1 eV, δ = 15, i.e., there is substantial neutrals density throughout the plasma radius, in agreement with simulations.At T e = 1.5, 2.0 and 2.5 eV one finds δ = 0.51, 0.01 and 0.004, respectively.The last two values are in fair agreement with the relative width across which neutrals density increases from 20 % to 100 % of the edge value.Uniform neutrals temperature has been assumed for simplicity in the model; if neutrals temperature increases inside the plasma, then δ is to be interpreted as a scale of pressure variation and the scale of density variation becomes even shorter than δ.
Results on degree of ionization, azimuthal current density from ( 8) and azimuthal ion velocity from ( 9) are shown in Figure 7 for the lowest temperature case and in Figure 8 for the highest temperature case.The degree of ionization in the plasma core is 40 % at T e = 1 eV and 100 % at T e = 2.5 eV, while absolute values of self-generated azimuthal current density and velocity increase substantially at higher temperature.

Dependence on Plasma Radius and Density
Substantially different profiles are found when changing the plasma radius.Profiles obtained with R = 0.14 m (and B z scaled to preserve magnetic flux) are shown in Figure 9.The ratio between edge and central plasma density (or pressure) is 0.5 in this case, whereas a lower value of 0.14 has been found in the simulation with R = 0.12 m (Figure 6).Broadening of density profiles with increasing R is due to lower current density and then to weaker pinch force, total current being fixed.Lower edge density is recovered in simulations at lower central density (3.3 × 10 20 instead of 4.5 × 10 20 ), as shown in Figure 10.The edge density depends non-linearly on central density, and no acceptable solutions are found (at the chosen radius and temperature) for n 0 < 3.2 × 10 20 m −3 .
Consideration of different plasma radii reveals an issue for any future 2D analysis.Envisaging simulations shown in Figures 6 and 9 as representative of different cross sections at different z of the actual pinch, it follows that pressure cannot be constant along field lines, in fact, if it is constant along the field line at the plasma center, then it changes along the field line running at the plasma edge.The reason is that from (13), the pinch force decreases with increasing plasma cross-section.Constant edge density at constant axial density can be restored in the presence of significant azimuthal current density throughout the plasma and not only in the outer region.Driving such current against resistive drag requires radial current j r = 0, which would be acceptable in 2-D since zero divergence can be preserved by variations of the axial current along z.
However, pressure variations along magnetic field lines remain a relevant issue, in fact the core plasma pressure has to build-up somewhere between the electrodes and the core region.Collisions with neutrals can give at most the pressure of the surrounding gas, which is smaller than the core plasma pressure.Other possible reasons of pressure variation along field lines are pressure anisotropy, inertial forces, and electrostatic forces.The second reason has been invoked in [12] for plasmas with similar local parameters but far from steady state.Electrostatic forces are likely dominant in axial pressure variation close to electrodes.The role of convective inertia and of charge separation will be investigated in future work.

Discussion
The model developed in Section 2 allows accounting for three species dynamics with simple and easily controllable calculations.The main result is that owing to the high PROTO-SPHERA densities, the plasma is fully ionized already at T e = 2 eV, well below the Ar first ionization energy.Neutrals penetration is limited to an external layer, in which a substantial radial plasma velocity develops to drain fresh charges produced by ionization.This radial flow across the applied magnetic field drives both azimuthal current and azimuthal ion flow.The latter is a likely cause of the good plasma spreading on the anode surface.The self-generated azimuthal current could be a seed for the formation of a torus encircling the pinch, but this subject requires further work, in particular the uniform B z assumption should be removed and different applied magnetic fields should be considered, including the presence of x-points, where magnetic reconnection can occur and magnetic flux can be transferred from the pinch to the torus.Another issue that calls for further work is force balance along magnetic field lines, in the presence of the radial pinch force which varies by an order of magnitude in the z direction.Sorting out this item requires consideration of electrostatic and inertial forces, as well as axially resolved diagnostic information on plasma density, temperature, and electrostatic potential.where k iN , k eN and k ei are ion-neutral, electron-neutral and electron-ion collision rate coefficients respectively; V e is electron velocity, E and B are electric and magnetic field respectively, and p denote pressures.Isotropic pressures are assumed.Using ∇ • (nVV) = V∇ • (nV) + n(V • ∇)V more familiar forms result − m e n i n e k ei (V i − V e ), (A6) + m e n e N k eN (V e − V N ) + n i n e k rec (m i (V i − V N ) + m e V e ), (A7) where electron inertia has been neglected.Please note that Spitzer resistivity η S = m e k ei /e 2 .Substituting the resistive term from (A9), the ion fluid momentum balance can be cast as − N (m i n i k iN + m e n e k eN + (m i + m e )n e k ion )(V i − V N ) + m e e N(k eN + k ion )(j − ρ e V i ).(A10) The combination of plasma-neutrals collision coefficients in the second line corresponds to the ambipolar diffusion coefficient of gas discharges with uniform temperature and uniform neutrals density D a = e(T e + T i ) N(m i k iN + m e k eN ) .

Figure 1 .
Figure 1.Schematic of PROTO-SPHERA target configuration, including annular electrodes, centerpost discharge in red color and quasi-spherical toroidal plasma (ST) in brown color.

TeFigure 3 .
Figure 3. T e = 1.0 eV.Upper frame: plasma (solid line) and neutrals (dashed line) number density profiles in 10 20 m −3 units.Lower frame: ion velocity in percent of the sound speed.

Figure 4 .
Figure 4. T e = 1.5 eV.Upper frame: plasma (solid line) and neutrals (dashed line) number density profiles in 10 20 m −3 units.Lower frame: ion velocity in percent of the sound speed.

TeFigure 5 .
Figure5.T e = 2.0 eV.Upper frame: plasma (solid line) and neutrals (dashed line) number density profiles in 10 20 m −3 units.Lower frame: ion velocity in percent of the sound speed.

TeFigure 6 .
Figure 6.T e = 2.5 eV.Upper frame: plasma (solid line) and neutrals (dashed line) number density profiles in 10 20 m −3 units.Lower frame: ion velocity in percent of the sound speed.

TeFigure 9 .
Figure9.T e = 2.5 eV and R = 0.14 m.Upper frame: plasma (solid line) and neutrals (dashed line) numer density profiles in 10 20 m −3 units.Lower frame: ion velocity in percent of the sound speed.

TeFigure 10 .
Figure10.T e = 2.5 eV, R = 0.14 m and lower density.Upper frame: plasma (solid line) and neutrals (dashed line) numer density profiles in 10 20 m −3 units.Lower frame: ion velocity in percent of the sound speed.
N ) + ∇ • (NV N V N ) = −∇p N + m i n i N k iN (V i − V N ) + m e n e N k eN (V e − V N ) − m i n e Nk ion V N + n i n e k rec (m i V i + m e V e ), (A4) m e ∂ ∂t (n e V e ) + ∇ • (n e V e V e ) = −en e (E + V e × B) − ∇p e− m e n e N k eN (V e − V N ) + m e n i n e k ei (V i − V e ) + m e n e (Nk ion V N − n i k rec V e ), (A5) m e n e ∂ ∂t V e + (V e • ∇)V e = −en e (E + V e × B) − ∇p e − m e n e N (k eN + k ion )(V e − V N ) + m e n i n e k ei (V i − V e ).(A8)Introducing charge density ρ e = e(n i − n e ) and current density j = e(n i V i − n e V e ), (A8) can be cast as a generalized Ohm's lawen e E = (j − en i V i ) × B − ∇p e − m e n e N (k eN + k ion )(V i − V N ) +m e e (N(k eN + k ion ) + n i k ei )(j − ρ e V i ), (A9)