Inference of field reversed configuration topology and dynamics during Alfvenic transients

Active control of field reversed configuration (FRC) devices requires a method to determine the flux surface geometry and dynamic properties of the plasma during both transient and steady-state conditions. The current tomography (CT) method uses Bayesian inference to determine the plasma current density distribution using both the information from magnetic measurements and a physics model in the prior. Here we show that, from the inferred current sources, the FRC topology and its axial stability properties are readily obtained. When Gaussian process priors are used and the forward model is linear, the CT solution involves non-iterative matrix operations and is then ideally suited for deterministic real-time applications. Because no equilibrium assumptions are used in this case, inference of plasma topology and dynamics up to Alfvenic frequencies then becomes possible. Inference results for the C-2U device exhibit self-consistency of motions and forces during Alfvenic transients, as well as good agreement with plasma imaging diagnostics.

A field reversed configuration (FRC) has no externally imposed toroidal field, belonging to the category of compact tori 1,2 . The poloidal field in an FRC has one component arising from magnets arranged on a common linear axis and another component generated by a toroidal plasma current flowing in opposite direction to the magnet currents. Under transient conditions, an additional magnetic field component arises from toroidal currents flowing in the vessel (the flux conserver (FC) current), which are induced by changes in plasma current distribution and/or transients in the external magnet currents. When the plasma current is strong enough to reverse the externally imposed magnetic field, a closed field structure topologically similar to a torus is formed (Fig. 1). Closed field lines in an FRC enclose a point of maximum kinetic pressure and null magnetic field called the o-point. The separatrix is the flux surface with null poloidal flux and separates the internal closed field region from the open field line scrape-off layer (SOL) region.
The C-2U device, built and operated by Tri Alpha Energy (TAE), is the first device to demonstrate that FRCs can be sustained in near steady state using neutral beam injection 3,4 . TAE's C-2U device relied largely on FC effects to stabilize plasma displacements, so the discharge lifetime was of the same order as the time constant of the vessel. TAE's C-2W device, presently in its initial operational phase, will extend the discharge duration over this limit, so plasma control will become necessary to stabilize the separatrix shape and position 5 . A method to determine the magnetic field structure and related control variables in real time (with sampling frequency in the range 10-100 kHz) is then required. Some first order approximations for FRC geometry parameters are available from the excluded flux radius. However, these cannot distinguish an FRC from a high beta mirror 6 , so they are not particularly useful if the FRC state itself is uncertain.
Determination of the magnetic field structure in an FRC is a challenging problem. While magnetic field structure inside the plasma can be measured by inserting probes inside the plasma 7 , this cannot be done in high temperature plasmas without severely disrupting the plasma confinement. In FRC plasmas, the magnetic field structure must be determined indirectly from external magnetic probes 8 , laser polarimetry systems 9 , etc.
Determination of the internal current sources from external measurements is termed the inverse problem. The technique used to determine the current sources from the sensor data is the inference technique. When Lorentz forces are balanced by plasma pressure, there is no net acceleration of the plasma, and the plasma is said to be in equilibrium. The determination of the magnetic structure corresponding to plasma in equilibrium is referred to in the literature as 'equilibrium reconstruction'. This work departs from the standard equilibrium reconstruction approach 10 and use instead the current tomography (CT) method 11,12 , a well-validated alternative already studied in connection with real-time control of tokamaks 13 . The CT method uses Bayesian inference 14 of Gaussian processes (GPs) 15 to solve the inverse problem. The GP modelling used by the CT method can be tailored to a multiplicity of related tomography problems 16 , in particular to the specifics of the FRC magnetics. There are several advantages of the CT method that make it ideal for plasma control. First of all, when the relationship between current sources and sensors is linear (such as is the case with magnetic probes), and the physics assumptions can be reduced to linear relationships among current sources and measurements, the solution depends on the sensor data through non-iterative matrix operations and, for this reason, is deterministic and suitable for real time. A version of the algorithm for C-2W device has already been implemented in a field-programmable gate array and verified to run under 10 μs. Second, as no equilibrium restrictions are necessarily required, the CT method can infer Alfvenic oscillations from magnetic sensor data. Fast transients can then be resolved accurately and with very low latency, both factors known to have an impact on control systems performance. Third, the CT method is able to fuse information from multiple sensor data sets and boundary conditions using a unified inference approach. This allows straightforward scalability should other magnetic sensors become available at a later stage. Sensors based on Polarimetry 17 and Hanle effect 18 , for instance, are both planned for TAE's C-2W device. Finally, the CT method provides uncertainty measures on all inferred outputs. This is interesting information on its own, but it has also an interest for advanced control applications, since the uncertainty information can be factored in as part of a robust control scheme 19 .

Results
Inference of Alfvenic transients in FRC. In the C-2U device 3 , two individual toroidal current rings are produced inductively (θpinch technique) in two opposing quartz formation sections placed at both ends of a stainless steel vacuum vessel (Fig. 2). These are produced simultaneously using pulsed power, fast magnetic field transients, and then accelerated out of their respective formation sections at supersonic speeds v z~3 00 km/s. Collisions of both FRCs take place inside the confinement vessel near or at the mid-plane z = 0.
The merging process occurs during the first few 10s of μs of the discharge, transforming the kinetic energy of the two initial compact tori into thermal energy of a single, static FRC 20 . Neutral beam heating is then applied to this initial FRC to provide the necessary heating and current drive to sustain the discharge against thermal and resistive flux losses.
When the accelerations in both formation sections are slightly different with respect to each other, FRCs do not collide exactly in the middle of the confinement section, leading to a merged FRC with a residual velocity. The resulting FRC is bounced back and forth in the axially stabilizing external mirror field until its position is stabilized around the machine mid-plane. Analysis of these oscillations provides a way to test the compliance of the inferred forces and accelerations with Newton's second law, using estimations of the plasma mass obtained by other diagnostics, as it will be shown. A plasma discharge exhibiting Alfvenic plasma oscillations around the midplane is chosen for the study, as illustrated in Fig. 3. Contours of poloidal flux and forward prediction of the actual magnetic measurements are shown every 10 μs. The frequency of the oscillations is about~20 kHz.
The axial position of the o-point is shown in Fig. 4 with high time resolution (every 10 μs), along with other geometric descriptors and plasma variables related to those. The o-point position is strongly correlated with the vessel current imbalance I ζ V , defined as: where I z>0 V is the net toroidal current flowing in one half of the vessel with z > 0, and I z<0 V is the net toroidal current flowing the other half of the vessel with z < 0. For a static plasma, the vessel current imbalance is zero. As plasma moves back and forth, midplane antisymmetric current components are induced in the vacuum vessel, which eventually dissipate ohmically. These are in the direction to oppose and slow down the plasma movements.
As a result, a strong correlation between the o-point axial position and the vessel current imbalance exists, as shown.
The separatrix radius is found to be proportional to the o-point radius, which is in agreement with Eq. (7). The trapped flux ψ 0 also matches approximately the approximation (9) for an elongated FRC. These approximations are not used in the inference process but as a check for consistency of the final results with these limiting cases.
The total number of deuterons in the plasma is estimated from line integrated density measurements integrated over the excluded flux radius. Plasma mass is estimated to be m p = 1.3 × 10 −7 kg from the deuteron mass times the deuteron inventory. The acceleration € z of the o-point can be determined from its position z (see Fig. 4). The net Lorentz force F z exerted over the whole plasma current distribution can be determined from the inferred current distribution and derived magnetic field. It turns out the product of the plasma mass and acceleration € z is consistent with the inferred electromagnetic force Formation section 1.2 within one standard deviation, as illustrated in Fig. 5. So Newton's second law is recovered from the inference results. The algorithm, however, is not very accurate during the first 50 μs or so of the discharge (right after formation/merging) presumably because the smoothing prior used cannot adequately describe the abrupt profiles resulting from shock waves or violations of other prior assumptions. Another test of relevance is to check whether the axial forces are proportional to some measure of plasma position z If Eq. (3) is valid for some axial range, a Hooke's constant can be defined. For a rigid plasma current distribution subjected to an infinitesimal displacement, the Hooke's constant can be evaluated from the plasma current distribution and the externally applied flux ψ ext (from magnets and FC currents) as an integral extending over the plasma domain Ω 21 Note that when taking derivatives the flux created by the plasma does not change with z, as the plasma is considered a rigid object; only the external flux does change due to the relative motion.
A positive Hooke's constant corresponds with a magnetic configuration that is axially stable and vice versa. The evolution  Axial force and displacement are linearly dependent in a range of +/−1 m around the mid-plane, so plasma dynamics can be approximated by a linear partial differential equation in this range. This is interesting for the future plasma control goals, since control theory and practice are well established for linear systems 22 . The Hooke's constant is positive due to the axially stabilizing external field and therefore consistent with an axially stable magnetic configuration that reaches the mid-plane after a few oscillations, as observed. The inferred value of about 1000 N/m is in agreement with the results obtained using the Lamy Ridge code 23 . The inference method can also provide the axial stability properties of the magnetic configuration. This is an important information for plasma control of future devices, which will require to establish and sustain an axially unstable plasma in equilibrium around the mid-plane z = 0 24 . A method to determine the axial stability properties of the magnetic configuration will be therefore required.
Comparison with plasma imaging. High-speed imaging of visible plasma emission is an independent technique that can yield information about the plasma dimensions. In this study, qualitative agreement between visible light emission from intrinsic oxygen impurity ions and the dynamics of the inferred poloidal flux contours serve as additional validation of the proposed inference method. Photons emitted from the 3d→3p transition (at 650.0 nm) of O 4+ were measured using a filtered high-speed camera with a radial view of the plasma 25 .
Emissivity of this spectral line was reconstructed (assuming axis symmetry) using the Simultaneous Algebraic Reconstruction Technique 26 . The core FRC electron temperature and density are more than sufficient to ionize the O 4+ charge state and populate higher charge states via electron impact excitation; therefore, minimal emission from this spectral line is found in the core. Instead, emission is peaked in the SOL where the electron temperature and density are lower and diffusive transport from boundary sources competes with ionization to higher charge states.
An example comparing the results of the magnetic inference method with the emissivity reconstruction for this spectral line is shown in Fig. 6. A relatively high-density plasma discharge (#48269) was chosen so that the emission measurement had good signal. Good agreement in the temporal dynamics of the reconstructed poloidal flux and emissivity is observed. This agreement provides further validation of the proposed inference method and is all the more encouraging since the two reconstructed quantities are derived from independent measurements (magnetics vs. photons) and analysis procedures.
The overall consistency of inferred results (Fig. 7) is also good, with the following highlights: (a) R s ¼ p 2R 0 is really a very good approximation, within one sigma. (b) The long FRC trapped flux (Eq. (9)) ψ 0 ¼ B w R 3 s =R w is also a very good approximation, being its overall magnitude in agreement with similar results obtained The approximation (14) does not reproduce well the inferred FRC length. (e) There is a correlation between FRC length and plasma current, as expected from Eq. (11). However, the approximation (11) does not reproduce well the inferred results, partly because this approximation does not consider a current distribution flowing outside the separatrix. (f) Vessel current decays in about~5 ms, comparable with the characteristic time over which the FRC is passively stabilized.

Discussion
We have used the CT method to provide a direct inference of the internal FRC magnetic topology, both during steady state and fast Alfvenic transients. The viability of the approach has been verified in a number of ways, including comparisons with approximate results from a long FRC approximation, recovering of a force balance dynamic equation, and comparison with imaging of visible plasma emission. All current sources have been modelled as GPs and inferred from external magnetic measurements using Bayesian analysis.
Smoothing priors (for plasma current and vessel current distributions) and a flux-conserving prior derived from Lenz's law (for the magnet currents) have been used in the inference. From all the inferred current sources, FRC topology and dynamic properties have been obtained. This includes the separatrix geometry and the axial stability properties of the magnetic configuration, among others.
When GP priors are used, and linear relationships among current sources and measurements can be established, the CT solution involves non-iterative matrix operations and is then ideally suited for deterministic real-time applications. Because no equilibrium assumptions are used in this case, inference of plasma topology and dynamics up to Alfvenic frequencies then becomes possible. The FRC topology and dynamics have been determined during Alfvenic oscillations, with excellent self-consistency of results.
Methods FRC approximations. The inference results of experimental data presented have been compared with first-order approximations for FRC parameters, which are summarized below. These are valid for an elongated FRC inside a FC of constant radius R w with negligible field line curvature at the mid-plane, termed the long FRC approximation.
The radial pressure balance condition relates the magnetic field component in the axial direction right beneath the inner vessel walls at the o-point plane B w with the average kinetic pressure of the plasma: So the average kinetic pressure of the plasma at the o-point (where the magnetic pressure is null) must necessarily be equal to the magnetic pressure at the confinement vessel walls When the plasma pressure is a flux function, and Eq. (5) is fulfilled, then the opoint radius is proportional to the separatrix radius 1.
In addition to the former, the plasma is in axial force balance, and the maximum beta achievable by an ideal FRC surrounded by a perfect FC of constant radius R w is given by the Barnes' average β condition 29 , which depends solely on the ratio of separatrix radius R s to FC wall radius R w . When both axial and radial pressure balance are fulfilled, the flux at the o-point (trapped flux) is given by The plane perpendicular to the machine axis that contains the o-point is termed the o-point plane. The intersection of the o-point plane with the machine axis determines the point were the axial component of the magnetic field is minimum. From Eqs. (5) and (9), the magnitude of the field reversal B ax at this point is From Ampere's law and Eq. (10), the total plasma current in a very elongated FRC of length L can be approximated by The plasma elongation is defined as A common approximation to separatrix radius and length comes from the excluded flux radius axial profile, which can be derived directly from magnetic sensors 30 , as explained below. If ψ i ; B i z are the flux and field profiles determined at positions z i along the internal wall of the vacuum vessel with radius r = R w , the excluded flux radius profile is defined as The position of the plasma mid-plane can then be estimated from the position where the excluded flux radius has its maximum.
A first-order approximation for the FRC separatrix radius is to consider it equal to the excluded flux radius at the mid-plane. A first-order approximation for the FRC x-point position is taken as the point along the axis Z 2/3 where the excluded flux radius has fallen to 2/3 of its value at the mid-plane 28 . The FRC length is then approximated by Bayesian inference of GPs. Bayesian inference is used in this paper to calculate the posterior distribution of currents given the magnetic measurements. The method, however, is generic enough to be used in a variety of related tomographic problems, which can be stated as follows. Given a forward model D = H(X) relating a continuous variable X(r) function of location r = (r 1 ,r 2 ,r 3 ) with some discrete set of measurements in the data vector D, the objective is to obtain all the solutions for X(r) that can explain the measurements in D. These are arranged into a probability distribution p(X|D) termed the posterior. A likelihood probability distribution p(D|X) measures the misfit between the model predictions H(X)and the measurements D. The probability of the spatial variable p(X) prior to taking any measurements is termed the prior probability distribution. According to Bayes theorem 14 , the posterior can be obtained from the likelihood distribution and the prior as The term in the denominator p(D) is called the evidence (or marginal likelihood) and normalizes the volume of the posterior distribution to 1.
Given prior and likelihood, the most likely solution is the maximum a posteriori (MAP) estimate, the solution in the posterior with the highest probability.
In the particular case where the forward model is linear, the spatial variable X(r) can always be discretized on a fine grid of dimension k, and a matrix K 2 R n k can be used to relate the discretized variable X 2 R k with a set of n measurements in Assuming additive Gaussian measurement noise ε = N(0,Σ D ) independent of X, the likelihood function can be modelled by an n-dimensional Gaussian distribution where Σ D 2 R n n is the data covariance matrix. The prior distribution can also be approximated by a multivariate probability distribution over X where Σ X 2 R k k is the prior covariance kernel and μ X 2 R k is the prior mean. It is usually convenient (but by no means necessary) to consider a zero mean μ X = 0 on the prior. The posterior distribution can likewise be approximated by a k-dimensional Gaussian probability distribution.
Since all probability distributions are Gaussian, the posterior distribution can be obtained analytically, since Gaussian distributions are transformed into Gaussian distributions through linear operations. The posterior mean (MAP estimate) and covariance are given in this case by 11 : As the dimension k of the multivariate normal distribution is made increasingly large, the multivariate normal distribution approaches a continuous distribution, and at this limit a GP is obtained 15 . In our case, the vector X becomes a continuous function X(r) of the spatial location. All possible solutions for X(r) can then be thought of as being generated by a stochastic process, described by the corresponding GP.
For a large number of situations in plasma physics, transport processes will work in the direction to reduce the spatial gradients of X(r). In other words, our prior belief about X(r) is that it must be a smooth function of r. The prior covariance kernel required for the inference can then be parameterized using the Squared exponential (SE) function, which is one of many options available 15 to model the spatial correlations between the values of a smooth profile variable at two points r and r': In the Bayesian context, σ and λ i are termed the prior hyper-parameters. The standard deviation σ controls the spread of values of X. The scale length λ i determines how quickly the plasma variable can change with the coordinate r i . A large length scale will give a large covariance between the values of the variable X at different r i coordinates, so the prior probability (Eq. (19)) for large differences between the values of the plasma variable X at neighbouring positions r i ; r′ i will be low. In other words, if the plasma profiles are smooth, the corresponding scale lengths will be large, and vice versa.
The SE kernel is by no means the only choice for a prior covariance kernel. A good review of GPs and the most common covariance kernels used can be found in ref. 15 . In general, the prior covariance kernel will have a set of hyper-parameters, arranged in a vector θ for simplicity.
Determination of the prior hyper-parameters can be considered as a continuous model selection problem, where the more likely hyper-parameters are obtained directly from the data 31 . The posterior for the hyper-parameters is p(θ|D), which from Bayes theorem is where p(D|θ) is the likelihood and p(θ) is the hyper-prior (prior for the hyperparameters).
In Bayesian model selection, the optimum set of hyper-parameters θ opt is selected to maximize this probability.
The prior over the hyper-parameters p(θ) in Eq. (24) is usually taken to be flat, since there is no indication of what are the best hyper-parameters before seeing the data. In this case, the optimal set of hyper-parameters that maximizes likelihood of the data with respect to the hyper-parameters is Given a set of hyper-parameters θ, there is an infinite class of plasma profiles X (r) that can be generated by the corresponding prior covariance p(X|θ) through the corresponding GP. The quality of the data fit must be evaluated not just for one particular solution but for all the solutions that can be obtained for a given set of hyper-parameters. The likelihood should be integrated out (marginalized) with respect to all these possible profiles generated by a single set of hyper-parameters, so it becomes a marginal likelihood.
In the particular case at hand where p(X) is a GP, the likelihood is normal and the model linear, the marginal likelihood can be calculated analytically. The expression for its logarithm is 16 For any given prior kernel, the maximum of the expression (28) with respect to the hyper-parameters gives the optimal set of hyper-parameters that explain D.
Inference model for the C-2U device. The C-2U magnetic model used for the analysis comprises a total of 42 magnets, 31 vacuum vessel (passive) segments and a current distribution made of 734 discrete plasma current elements modelled as block coils (Fig. 8). Of special relevance are 8 equilibrium (EQ) magnets in the confinement vessel and 6 FC magnets, which can be used as passive FCs or be connected to power supplies. The magnetic measurement system on C-2U 32, 33 comprises a set of 19 magnetic pick-up probes placed inside the confining vessel and 8 external pick-ups located right underneath the 8 EQ magnets (Fig. 8). There are also Rogowski-based current measurements for all the FC magnets currents I EQ and also for some of the EQ magnets currents I EQ . For the rest of the magnets, only the set point used for its control is known.
The inference problem at hand requires finding the most likely solution for the elements of the total plasma current distribution arranged in a vector I P , along with the most likely solution for current induced in the confining vessel I V and all the magnets I M . A diagram illustrating the magnet location and grid used for the current distribution is shown in Fig. 8.
All the currents to be inferred are arranged into a single current vector All current sources are modelled as GPs as described earlier. The information available to perform the inference comes from (i) set points for all I M , (ii) current measurements for I FC and a few I EQ , (iii) measurements of magnetic field at several locations outside the plasma region, both inside and outside the confining vessel, (iv) null boundary conditions for plasma current distribution, and (v) null boundary conditions for the flux change underneath the equilibrium magnets ∂ψ ∂t ffi 0, which behave as perfect FCs on the timescale of the discharge.
The boundary conditions (iv) and (v) are built directly into the prior, to obtain solutions where the plasma current distribution drops to zero at the domain, and the flux is conserved at the magnet locations (flux-conserving prior).
From the inferred currents in I is then straightforward to calculate the poloidal flux and magnetic field components on the domain grid using the matrix representations M,G R ,G Z of the Biot-Savart operator 34 : Main plasma shape and position variables of interest for control such as xpoint, o-point and separatrix radius can then be obtained directly by searching for nulls on the magnetic field and flux along the axis and mid-plane. Low-order moments of the plasma current distribution of interest for control such as total plasma current or the axial position of current centroid can likewise be obtained from linear operations.
Data availability. All relevant data supporting the findings of this study are available from the authors on request.