J. Plasma Fusion Res. SERIES, Vol. 9 (2010)

Model Data Fusion: developing Bayesian inversion to constrain
equilibrium and mode structure
M. J. Hole, G. von Nessi, J. Bertram, J. Svensson †, L. C. Appel ‡, B. D. Blackwell, R. L. Dewar andJ. Howard
Research School of Physics and Engineering, Australian National University, ACT 0200, Australia.
†Max Planck Institute for Plasma Physics, Teilinstitut Greifswald, Germany.
‡Euratom/UKAEA Fusion Association, Culham Science Centre, Abingdon, Oxon OX143DB, UK.
March 1, 2010

(Received: 10 November 2009 / Accepted: 1 March 2010)

Recently, a new probabilistic “data fusion” framework based on Bayesian principles has been developed on
JET and W7-AS. The Bayesian analysis framework folds in uncertainties and inter-dependencies in the diagnostic
data and signal forward-models, together with prior knowledge of the state of the plasma, to yield predictions ofinternal magnetic structure. A feature of the framework, known as MINERVA (J. Svensson, A. Werner, Plasma
Physics and Controlled Fusion 50, 085022, 2008), is the inference of magnetic flux surfaces without the use of aforce balance model. We discuss results from a new project to develop Bayesian inversion tools that aim to (1)distinguish between competing equilibrium theories, which capture diﬀerent physics, using the MAST sphericaltokamak; and (2) test the predictions of MHD theory, particularly mode structure, using the H-1 Heliac.
Keywords: Bayesian inference, equilibrium, tokamak, stellarator

1. Introduction

data D. A widespread technique used is least-square fitting, in which prior assumptions are included via a penalty
term in the fit.
The confluence of higher performance plasmas withdiagnostic development has however led to a dichotomy:
data is often inconsistent with ideal MHD equilibria, sometimes misleading scientists to propose new phenomena to
describe data “artifacts”. An example is filamentation offlux surfaces in the Rijhnuizen Tokamak, which was inferred from fluctuations in measurements of electron temperature [3]. Subsequent analysis of the Thomson scattering detection chain revealed that with the correct photoelectron statistics and 2D instrument profile, similar structures could arise from noise. [4, 5].
Recently, a new integrated data-modelling approach
for inference of fusion plasma parameters has emerged
which oﬀers a natural framework with which to resolve artifacts from model behaviour. In contrast to least square fitting, the Bayesian approach to inference in fusion plasmas,
developed by multiple authors, [6, 7, 8, 9, 10, 11, 12, 13]
involves the specification of an initial prior probability distribution function (pdf), P(I), which is then updated bytaking into account information that the measurements provide through the likelihood pdf P(D|I). The result is theposterior distribution P(I|D) given by Bayes’ formula

With the advent of large scale neutral beam heating, magnetically confined fusion plasmas have driftedfrom the simple picture of ideal magnetohydrodynamics(MHD), which describes the plasma as a single, stationary, isotropic Maxwellian fluid. Due principally to neutral beam heating, several tokamak experiments now boast
plasma toroidal rotation speeds that approach the thermal
Mach speed and have significant stored energy residing in
the energetic particle population produced by charge exchange with thermals. Motivated by these developments, a
range of new descriptions have emerged that include thermal rotation [1] as well as energetic particles [2]. Despite this, ideal MHD is still the foundation of nearly all
analysis. Detailed magnetic reconstruction based on thistreatment ignores the energetic complexity of the plasma,and can result in model-data inconsistencies, such as thermal pressure profiles which are inconsistent with the total
stored kinetic energy of the plasma.
A parallel development has been the improvement in
the diversity, accuracy and resolution of plasma diagnostics. Interpretation, however, often requires a detailed
knowledge of the plasma equilibrium. For example, inference of the toroidal current profile jφ (ψ) from line of sightmeasurements of the polarization angle requires a knowledge of the poloidal flux ψ across the plasma. Formally,
diagnostic forward functions relate the vector of plasma
parameters I to the measurement vector D. For a linear
system, such as toroidal current inference in a double null
configuration, I and D are normally related through a response matrix M with additional contributions C, such that
D = MI+C. Inference involves inverting this relationshipto give plasma parameters I that are consistent with theP(I|D) = P(D|I)P(I)/P(D).

(1)

The advantage of the Bayesian approach over traditional
inversion techniques is two-fold: (i) prior knowledge, including known parameter inter-dependencies is made explicit, and (ii) as the formulation is probabilistic, random
errors, systematic uncertainties and instrumental bias areintegral part of the analysis rather than an afterthought.
In this work we present initial results of forward mod-author’s e-mail: matthew.hole@anu.edu.au479©2010 by The Japan Society of Plasma
Science and Nuclear Fusion Research

M.J. Hole et al., Model Data Fusion: Developing Bayesian Inversion to Constrain Equilibrium and Mode Structure

els of current tomography in the Mega Ampere Spherical Tokamak (MAST), we formulate Bayesian inferenceof force balance, and we identify plans for inference ofmode structure in the H1 heliac. The paper is structured
as follows: Sec. 2 briefly outlines the MAST experimentand present first results for Bayesian inference of poloidal
flux and current profiles from MAST motional Stark eﬀect(MSE) measurements. Section 3 outlines a Bayesian inference framework for force balance, synthesizes a pressure profile posterior given the flux surfaces from MSE
and Thomson Scattering data, and computes an estimateof f (ψ) f  (ψ), with f (ψ) the toroidal flux function. Using
this, we are able to quantify the impact of poloidal currents on the equilibrium configuration. Section 4 outlines
a Bayesian inference model for mode structure and develops a plasma model to compute candidate modes in the H-1
heliac. Finally, Sec. 5 contains concluding remarks.

2.5

2

1.5

1

Z (m)

0.5

0

−0.5

2. MAST and current tomography−1

The Mega Ampere Spherical Tokamak, which firstpublished physics results in 2001 [14], is one of the world’s
largest spherical tokamaks. While MAST physics andtechnology development have contributed across a broad
range of fusion science [15, 16], two properties of MAST
help motivate this research: high performance, and precision diagnostics. It is the combination of these properties in a spherical tokamak, which has a relatively weaktoroidal field and therefore large sensitivity to poloidal
currents compared to conventional tokamaks, which helps
motivate selection of MAST to develop Bayesian inference
for force balance.
Recently, both MAST neutral beam injectors havebeen upgraded to 3.8 MW. This has enabled plasma performance to be routinely lifted above βn ≈ 5, which wasreported in 2005 [17]. MAST is equipped with an array of
precision diagnostics [16], including a high spatial resolution, single-time point, ruby Thomson Scattering system, amulti-time point Nd:YAG Thomson Scattering system witha coarser spatial resolution, Motional Stark eﬀect, Charge
Exchange Recombination, and fast magnetics [18].
Our development of Bayesian inference of the current profile on MAST closely follows the seminal workof Svensson and Werner [11]. In that work, the plasma
was represented as a grid of toroidal axisymmetric current
beams, each with rectangular cross-section and each beam
carrying a uniform current density. In MAST, we haveplaced these beams so as to fill-out the entire limiter region as depicted in Figure 1. The magnetic field generated
is a summation of Biot-Savart’s law over current beams.
A key advantage to using a series of current beams
with finite cross-section to model the plasma current (asopposed to a filamentary model) is that the semi-analytic
expressions for the corresponding magnetic field and vector potential have no singularities, even at points within thecurrent beam itself. Indeed, if one were to use filaments to

−1.5

−2

−2.5

0.5

1
R (m)

1.5

2

Figure 1: Plasma beam cross-sections fill out the MAST limiter
region in the plasma beam current model. Pick-up coils
are shown as solid bars near the poloidal flux coils, flux
loops as a “*” on both the centre column and poloidal
field coils, and MSE as “*” points on a chord through the
midplane.

model the plasma current, there would be many singularpoints in the calculated magnetic field within the plasmathat would make subsequent flux-surface calculations difficult and somewhat questionable [11].

2.1

Diagnostic Signals

In axisymmetric devices, there exists a direct relation
between the poloidal flux, ψ, and the toroidal componentof the current vector. Specifically, defining the cylindricalcoordinate system (R, Z, φ) with R the major axis, Z the vertical axis and φ the geometric toroidal angle, the following
relation holds:

ψ(R, Z) =480



B · dS =

A · dl = 2πRAφ ,(2)

M.J. Hole et al., Model Data Fusion: Developing Bayesian Inversion to Constrain Equilibrium and Mode Structure

where Aφ is related to the toroidal current jφ via BiotSavart’s law:jφ (r )µ0(3)dV  .
Aφ (r) =4π|r − r |and where µ0 is the permeability of free-space. In Eq. (3) ris the position vector r = ReR +Zez +φeφ , which eR , eZ , andeφ unit vectors in the R, Z and φ directions, respectively. By
inspection, the AR and AZ components due to jφ are zero.
Taking the curl of Eq. (2) and utilizing B = ∇ × A it canbe shown
1 ∂ψ
,(4)
BR = −
R ∂Z
1 ∂ψ
BZ =
,
(5)
R ∂R
where BR , BZ and Bφ are the components of the magnetic
field. From Eqs. (2)–(5), it can be seen that if jφ is represented by an axisymmetric beam of rectangular cross section and uniform current distribution, then ψ, BR and BZ
will vary linearly with respect to the total current going
through that beam. The importance of this will become
clear in the next two paragraphs.
The main diagnostics used to infer the current distribution are pickup coils Pi , full flux loops Fi and the polarisation angle γi of the emitted light from neutrally excitedspecies during neutral beam injection due to the motional
Stark eﬀect (MSE). The responses of these diagnostics toa current running through a beam, I, is given byThus, given that Pi (R, Z; I), Fi (R, Z; I) andtan γi (R, Z; I) all have a linear dependence on thecurrent flowing through I, one may write out a generalized
prediction vector P for all the pickup coils, flux loops andpolarization angles, as
P = MI + C.

Here, M is the response matrix of the current vector I corresponding to all the plasma beams modelling the plasma,and C represents various other contributions to the prediction vector, which are constant relative to I.

2.2

tan γi (R, Z; I) =A0 BZ (R, Z; I) + A1 BR (R, Z; I) + A2 Bφ (R, Z; I),
A3 BZ (R, Z; I) + A4 BR (R, Z; I) + A5 Bφ (R, Z; I)

CAR Prior and Bayesian InferenceIn [11], Svennson and Werner used a Conditional
Auto-regressive (CAR) prior to perform Bayesian inference for the current beam model of the plasma, and thisis the model adopted in this paper. The advantage of thischoice in prior is that plasma beam currents are correlated
to each other but in a spatially localised way. Thus, thisprior has the eﬀect of enforcing some smoothness between
adjacent current beams, while minimising spatially longranged eﬀects due to manipulating the current in a particular beam. This is desirable in that such behaviour allows one to get a clear relation between individual current beams and diagnostic measurements. To construct thisprior distribution, the following distribution over all current beams is considered


1
(11)
P(I) ∝ exp − IT QI ,2

Pi (R, Z; I) = BR (R, Z; I) cos(θi ) + BZ (R, Z; I) sin(θi ), (6)Fi (R, Z; I) = ψ(R, Z; I),

(10)

(7)

where the superscript T denotes the transpose, and Q is
the precision matrix. Equation (11) is proportional to azero-mean normal distribution, and so the conditional distribution of one current Ii given all others I−i satisfies





2
βi j I j ) /(2τi ) ,(12)p(Ii |I−i ) ∝ exp (Ii −(8)where θi is the angle between a pickup coil’s normal andthe outboard midplane, and A0 , A1 , A2 , A3 , A4 and A5 areconstants for the particular MSE viewing geometry. The
notation “;” denotes that the subsequent argument, in thiscase current I, is held constant. (By convention, the notation “|” is reserved for probability theory.) Both Pi and Fi
have a linear dependence on ψ, BR and BZ and hence they
are linearly dependent on the total current going througheach beam of the current plasma beam model. The function tan γi has a non-linear dependence on BR and BZ ; however given that both neutral beams and MSE viewing optics
all are on the mid-plane for MAST, A1 , A2 , A3 , A4 are allapproximately zero. The following approximation is hence
used for tan γi :A0 BZ (R, Z; I) + A1 BR (R, Z; I).
(9)tan γ(R, Z; I) =A5 Bφ (R, Z; I)jfor some βi j . Matching terms between Eq. (11) and Eq.
(12) shows that βii = 0, βi j = −Qi j /Qii and τi = 1/Qii .
Moreover, it turns out the symmetry of Q implies thatβi j τ j = β ji τi [19]. We have set all variances equal. By
setting βi j to correspond to Ii having a mean that is simply
the average of the currents of the beams horizontally andvertically adjacent to it, the precision matrix takes the form1
1
1− W ,(13)Q=τ
4

where 1 is the identity matrix and W the adjacency ma¯trix with Wi j = 1 if i and j are horizontally or verticallyadjacent to current element Ii , and Wi j = 0 otherwise.
From Eq. (12), it can be seen that the CAR prior isconstructed by exactly specifying the conditional distribution p(Ii |I−i ) and subsequently transforming this expression
into a zero-mean normal distribution for I. CAR distributions are discussed in detail in [19] and have the advantage

where we have retained the term A1 . Since the current
beam plasma model does not take into account poloidal
currents, vacuum field values for Bφ are used in Eq. (9) toperform the current tomography calculations. The correction to Bφ due to poloidal currents is investigated in Sec.
3.1.

481

that hidden parameter inter-dependencies associated withdirect manipulation of the covariance matrix in a normal
distribution are totally avoided in the CAR construction(see [11] for a discussion on this point).
We have used a normal distribution for every diagnostic measurement D as it is the least informative, or maximum entropy distribution when only the standard deviation (“error”) and mean value is available for each diagnostic signal [20], as is the case for MAST diagnostic data
signals. Under this assumption, the likelihood P(D|I) becomes a multivariate normal distribution of the data D, andis given explicitly by


1
1
T −1
(P
−
D)
exp
−
Σ
(P
−
D)
(14)
p(D|I) =2(2π)Nd /2 |ΣD |1/2

−1.5

−1

−0.5

psi (webbers)

M.J. Hole et al., Model Data Fusion: Developing Bayesian Inversion to Constrain Equilibrium and Mode Structure

0

2.5

2

1.5

1

Z (meters)

0.5

where ND is the number of measurements, and Σ is the covariance matrix of the measurements, which is determined
experimentally. Finally, substituting Eq. (11) and Eq. (14)into Bayes’ formula shows that the posterior distributionp(I|D) satisfies the following relation:p(I|D) ∝ exp −(P − D)T Σ−1 (P − D) − IT QI . (15)

0

−0.5

−1

−1.5

Recalling Eq. (10), it is known that P is linearly dependent on I; and thus, it can be shown that Eq. (15)is proportional to a multivariate normal distribution in I
(see [11] for details). Hence, p(I|D) has a simple analytic
representation that be directly analysed without resorting
to using Markov-chain Monte-Carlo (MCMC) algorithms.
The speed aﬀorded by the analytic nature of the posterior
makes this analysis amenable to real-time plasma control
applications.
With the maximum and variance of p(I|D) determined, Eqs. (2) and (3) can be used to construct thepoloidal flux surfaces corresponding to the contours of
ψ(R, Z). Figure 2 shows poloidal flux surfaces from MAST
discharge # 22254 at 320 ms using pickup coils, flux loops
and MSE. Discharge #22254 is a deuterium plasma in a
double-null configuration, which was heated with 3.1 MWof neutral beam heating and a plasma current of I p =800 kA. The time of 320ms is analysed here as it corresponds to a firing of the high-resolution TS system closest
to the peak β for this shot. The figure shows a contour
plot of ψ(R, Z) which is calculated from I corresponding
to the maximum of the p(I|D) distribution. Overlaid on
the contours are traces of the poloidal field coil cross sections and conducting surface cross sections for the MAST
experiment. The last closed flux surface calculated fromthe plasma beam model is outlined in bold with the corresponding EFIT last closed flux surface overlaid in purple
for comparison. One outcome of the Bayesian approach is
generation of pdfs from which the error to the fit can beinferred. Figure 3 shows the inferred safety factor, q, profile and distribution obtained by sampling the posterior 200
times. The width of the fits is a characteristic measure of
the error in q.

−2

−2.5
0

0.5

1
R (meters)

1.5

2

Figure 2: Poloidal flux surfaces inferred for MAST shot # 22254
at 320 ms using pickup coils, flux loops and MSE.

3. Inference of force balance

The physics goal of our work is to exploit the improved resolution of diagnostics to infer the validity of different force balance descriptions. Specifically, we have inmind development of a framework that will validate energetic plasma resolved force balance models. To build
towards this goal, we commence by applying Bayesian inference to the axisymmetric Grad-Shafranov equilibrium
description of ideal MHD force balance, J × B = ∇P.
The well known Grad-Shafranov equation [21] can be
written
F(R, Z) = 0

(16)

where
F(R, Z) = −µ0 R jφ (ψ) + µ0 R2 p (ψ) + µ20 f (ψ) f  (ψ), (17)−µ0 R jφ (ψ) = R∂ 1 ∂ψ ∂2 ψ.
+
∂R R ∂R ∂z2(18)

For a real plasma, the presence of nonideal eﬀects will
mean F(R, Z) is nonzero. Our long term aim is to compute P(F |D) using Bayes formula, such that P(D|F ) =P(D|F )P(F )/P(D). Here,the vector D comprises all the482

M.J. Hole et al., Model Data Fusion: Developing Bayesian Inversion to Constrain Equilibrium and Mode Structureof the magnetic axis and q profile from the recomputed solution, we are able to quantify the impact of poloidal currents on current tomography without the need for a separate magnetic reconstruction by EFIT. This folds p (ψ) intojφ (ψ) through the influence of poloidal currents.
For an ideal gas, the pressure is given by

12

10

q

8

6

p = ne kB T e + ni kB T i

4with ne , ni the electron and ion density, and T e , T i the electron and ion thermal temperature. At 320 ms, high resolution Thomson scattering measurements of T e and ne areavailable. The data is distributed normally about each data
point, and a mean and standard deviation provided. To obtain a pressure estimate, we have assumed that the meanvalues of density and temperature satisfy ni /ne = 0.8, andT i /T e = 1.1, as typically extracted from a charge exchange
recombination measurement. We have also assumed thesame distribution of data as Thomson Scattering.
We have next computed a fit for the pressure p(ψ) using Monte Carlo simulation [22]. The procedure is as follows: at each radial grid point, the inverse transformation
method is used to generate ne , T e , ni , T i samples that satisfying the prescribed pdfs P(ne ), P(T e ), P(ni ), P(T i ). To
compute p(ψ) and its distribution we have mapped eachsample p(r) to p(ψ) using ψ(r) across the inboard chord
determined from MSE. To find a smoothed p(ψ) we havethen fitted a fourth order polynomial in normalized flux.
The sampling process is repeated until the pdf for p(ψ) nolonger changes.
Figure 4 shows the sampled pressure and fit for thepressure using N s = 2000 samples. Overlaid is the polynomial fit for each sample across the midplane. The widthof the polynomial fit (shown in black) yields a lower caseestimate to the standard deviation in p(ψ).

2

0
0

0.2

0.4
0.6
Normalised Flux

0.8

1

Figure 3: Safety factor, q profile as a function of normalised
poloidal flux found by sampling the posterior 200 times
for shot #22254 at 320ms. The poloidal flux is normalised such that ψn = 0 is the magnetic axis and ψn = 1
is the edge.

magnetics, motional Stark eﬀect, Thomson scattering andcharge exchange recombination data, while the column
vector F contains F(R, Z) evaluated from Eq. (17) at different (R, Z) across the plasma.
An integral representation for the prior P(F ) can beobtained by using the transformation property ∂(F , p (ψ), f (ψ) f  (ψ)) =P(F , p (ψ), f (ψ) f  (ψ) ×  ∂( jφ (ψ), p (ψ), f (ψ) f  (ψ)) P( jφ (ψ), p (ψ), f (ψ) f  (ψ)) (19)and integrating both sides across p and f. That is,P( jφ (ψ), p (ψ), f (ψ) f  (ψ))  d p d f f  (ψ)
P(F ) =f (ψ) f  (ψ))  ∂(F ,p (ψ),∂( jφ (ψ),p (ψ), f (ψ) f  (ψ)) 

(21)

(20)

where the integrand d p d f f  = d p (ψ)d f (ψ) f  (ψ). An
ideal MHD plasma satisfies Eq. (16), such that the residual force is zero, and so P(F ) = δ(F ). Unfortunately,
the joint distribution function P( jφ (ψ), p (ψ), f (ψ) f  (ψ))
is non-separable, preventing direct integration of Eq.
(20). In future work we will fold interdependency ofjφ , p and f (ψ) f  (ψ) together to enable calculation ofP( jφ (ψ), p (ψ), f (ψ) f  (ψ).

3.1

Inference of toroidal flux function

Some progress can be made if we assume the plasma
obeys ideal force balance, and p (ψ) is assumed to be independent of f (ψ) f  (ψ) and jφ (ψ). In this instance, andproviding we are able to estimate p (ψ) and its distribution from measurements, then f (ψ) f  (ψ) across the midplane can be computed through Eq. (17). As a first step
to inference of force balance in a real plasma, we computef (ψ) f  (ψ), integrate to find f (ψ), and substitute this backinto MINERVA. By examining the changes in the position

Figure 4: Inferred pressure profile p(r) across inboard midplaneradial chord at 320 ms in # 22254. The sampled data areshown as points, and the black lines are fits to p(ψ) for
each sample set using a fourth order polynomial in ψ.

Using the pressure profile fit of Figure 4, we havesampled p (ψ) and jφ across the midplane and computed

483

M.J. Hole et al., Model Data Fusion: Developing Bayesian Inversion to Constrain Equilibrium and Mode Structuref (ψ) f  (ψ). Figure 5 shows the mean and standard deviation for p (ψ), jφ and f (ψ) f  (ψ). For reference, the scheduler EFIT solution for this time slice is over-plotted in bold.
While the MSE current profile is qualitatively similar to
EFIT, the pressure profile is very diﬀerent, particularly inthe edge and core regions. This diﬀerence is principally
responsible for the diﬀerence in f (ψ) f  (ψ) at the core andedge seen in panel (c).
Next, we have integrated f (ψ) f  (ψ) to obtain f (ψ)
from panel (c). This yields a toroidal flux profile thatvaries approximately linearly in poloidal flux from f (ψ) =−0.407 at the edge to f (ψ) = −0.496 at the core. That is,poloidal plasma currents are paramagnetic, and increase
Bφ by 22% at the core. Qualitatively, this poloidal current
should lift the on-axis safety factor ≈ 22%, increase thechange in poloidal flux across the plasma, but not change
the geometry of flux surfaces or the Shafranov shift. Workis in progress to investigate the impact of the correction tof (ψ) in MINERVA.

corresponding to a given plasma mode we obtain the forward model for D given I. Given some prior P(I), thelikelihood P(D|I) can then be constructed.
In the remainder of this section we derive the forward model for Global Alfvén Eigenmodes (GAEs) in H-1
plasmas. As a starting point to this sub-task, we compute
GAEs in a cylindrical plasma model (coordinates (r, θ, z))incorporating a vacuum region in the region r p < r < rw
which is encased by a perfectly conducting wall at r = rw .

Figure 6: Rotational transform as a function of average minor
radius for the typical range of κh in H-1NF, a three fieldperiod heliac with average minor radius r p = 0.2m andmajor radius R = 1m.

GAEs are discrete modes that accumulate at the minimum ω2min of the Alfvén continuum ω2A = k2 v2A [23] wherek denotes the component of the wave vector parallel tothe equilibrium magnetic field and vA denotes the Alfvénspeed. Expanding k = k · B and substituting into ω2A =k2 v2A yields ω2A ∝ (n − -ιm)2 /ρ where m, n are the poloidal
and toroidal mode numbers respectively, and ρ is
 the mass
density profile, assumed to take the form ρ = ρ0 1 − r2 /r2p
with ρ0 the density on-axis. We define ω2min to be the global
minimum of ω2A over the interval (0, r p ), if such a minimum exists. That is, ω2min is the smallest ω2A (r0 ) with both
ω2A (r0 ) = 0 and ω2A (r0 ) > 0 for some radius r0 suchthat 0 < r0 < r p . Figure 7 shows the dependency ofωmin on κh for diﬀerent m, n modes. We have restricted
our attention to the frequency interval 0 < f < 40 kHz,
where f = ω/(2π), since we are only interested in stable modes, and current H-1 diagnostics cannot resolve frequencies above 40 kHz. Theoretically, for a given κh ,GAEs lie at frequencies just below the lines in Figure 7,
accumulating at ωmin , with an increasing number of radial
nodes the closer the eigenfrequency lies to ωmin .
We have modelled the plasma region by a stellarator
normal-mode formulation in a cylindrical plasma [24, 25],
corrected for non constant density profile. The eigenvalue
equation for this problem can be written

Figure 5: Inferred parameters across the midplane inboard radialchord using Monte Carlo simulation. Panel (a) shows thepressure gradient p (ψ), panel (b) the toroidal current andpanel (c) the toroidal flux flux function f (ψ) f  (ψ).

4. Inference of mode structure

A second objective of the project is to test predictions of MHD mode theory, particularly mode structure.
To achieve this goal we intend to develop a Bayesian inference tool to identify oscillations using the mode structure of a collection of candidate modes. The low beta heliac plasma configuration of H-1NF is ideally suited to thistask, as there is little uncertainty in the equilibrium, andprecision control of the helical coil currents enables accessto a wide range of magnetic rotational transform ( -ι) profiles, as shown in Fig. 6. The combination thus oﬀers arich testing ground for Bayesian mode analysis.
For inference of mode structure, the relevant plasma
parameters I are the set of eigenmode-specifying parameters (frequency, mode numbers, radial structure moments)
and D is a set of Mirnov array measurements made outside the plasma. By solving the eigenmode problem in theplasma, and computing the vacuum-region magnetic fieldLϕ = ω2 Mϕ,(22)

where ω is given in units of inverse Alfvén time τ−1A =√B0 /R µ0 ρ0 , with B0 the field strength at the magnetic axis,484

M.J. Hole et al., Model Data Fusion: Developing Bayesian Inversion to Constrain Equilibrium and Mode Structure

If ω2min appears in the plasma and is positive, a shooting method is used to find solutions to Eqs. (22), (28),(31), and Eq. (29) or Eq. (30) with frequency 0 < ω <ωmin . More explicitly, we select the appropriate boundary condition for ξr at r = 0, and then ’shoot’ out to theplasma/vacuum interface, searching frequency space in theregion 0 < ω < ωmin for eigenfrequencies that are consistent with Eq. (28). Figure 8 shows an example of an(m, n) = (3, 4) eigenmode found in this manner, locatedat the point (κh , ω) indicated by the cross-hairs in Figure
7. The global structure seen has promising features whencompared to recent spectral measurements of the mode:the mode structure is global in radial extent, and the radial
position of the peaks nodes broadly matches observations.

Figure 7: A plot of the minimum in the continuum ωmin /(2π) vsκh for a selection of mode numbers (m, n). The cross-hairsmark the location of the eigenmode shown in Fig. 8.

and where ϕ is defined in terms of the radial element of thefluid displacement

n 
(23)
rξr = exp i(mθ − z) ϕ(r),R√with i = −1. The operators L and M are give by:
L

M

=

=

d
1 d
(n − m -ι)2 r +
r dr
dr

-ϊ
m2 
2
−
D
(n
−
m
ι
)
,
+
(n
−
m
ι
)
S
m
r2
ρ(r) m2
1 d ρ(r) dr +.
−r dr ρ0 dr
ρ0 r 2

−

(24)
(25)

Figure 8: The radial fluid element displacement ξr as a function of minor radius for (m, n) = (3, 4) and κh = 0.73
at the eigenfrequency indicated This eigenmode has angular frequency ω = 12.1 k rad s−1

Equation (22) is derived by averaging over helical ripple
and assuming large aspect ratio. The Suydam stability parameter is given byβ0 
p (r)Ω (r)2 2
where the average field line curvature is



2
2
Ω(r) =  N r -ι + 2 r -ιdr
DS = −

5. Conclusions

In this work we have described a Bayesian inversion
framework for inference of force-balance and mode structure. Progress has been reported in three areas: development of Bayesian inversion for current profiles; inferenceof the toroidal flux profile assuming the plasma satisfied
ideal MHD; and calculation of a Global Alfvén eigenmnodes in a helical-ripple averaged cylindrical plasma
model.
Based on working by Svensson et al. we have developed Bayesian inversion for current profiles and poloidal
flux contours in MAST. The model comprises a grid ofrectangular toroidally extended current beams, and is constrained by pickup coils, flux looops and MSE data. Owing
to the placement of the MSE viewing optics the tangentof the polarisation angle of emitted light is linear in BR
and BZ . Signals from the pickup coils and flux loops arealso linear in BR and BZ . This feature means the prediction
vector for the data is linear in the currents. We have useda conditional autoregresssive prior to describe the initial
current distribution, which has the advantage of enforcing
some smoothness between adjacent current beams while
minimizing spatially long ranged eﬀects. Posterior calculations of the poloidal flux contours of MAST #22254 at
320 ms illustrate good boundary agreement to EFIT.

(26)

(27)

√
Here j = −1, ’prime’ denotes radial derivatives, -ι̇ = r -ι ,p(r) is the normalized pressure,  is the inverse aspect ratio,
N is the number of turns made by the helical windings andβ0 = 2µ0 p0 /B20 is the ratio of plasma pressure to magnetic
pressure at magnetic axis.
Turning now to the vacuum region, the magnetic field
B must satisfy ∇×B = 0 and ∇·B = 0. In cylindrical geometry we can write B in terms of modified Bessel functions,
up to a scale factor.[26] The boundary conditions[[Br ]] = [[B ]] = 0(28)where [[B]] = Bvacuum − B plasma , together withξr (0) =ξr (0)=Br (rw ) =0, |m|  1,

0, |m| = 1,
0,

(29)
(30)
(31)

relate the vacuum solution to that in the plasma.

485

M.J. Hole et al., Model Data Fusion: Developing Bayesian Inversion to Constrain Equilibrium and Mode Structure

In Sec. 3 we outlined an inference technique forthe extraction of force balance, as applied to the Grad
Shafranov equation. Using Thomson scattering data and aMonte Carlo simulation technique, we fitted pressure profiles to the inboard midplane chord of MAST, and computed p (ψ). Assuming ideal MHD force balance we thencalculated f (ψ) f  (ψ). While the f (ψ) f  (ψ) profile is similar to EFIT, the change in f (ψ) from the core to the edge is22%, which is slightly larger than the EFIT value. Qualitatively, this poloidal current should lift the on-axis safety
factor by ≈ 22%, and increase the change in poloidal flux
across the plasma, but not change the geometry of flux
surfaces or the Shafranov shift. At present, we are substituting the corrected Bφ profile inferred from f (ψ) into
MINERVA to compute the eﬀect of poloidal currents on
the calculation of jφ and ψ surfaces. This constitutes a firststep towards combining jφ , p (ψ) and f (ψ) f  (ψ) to validate
diﬀerent equilibrium descriptions.
Finally, we have developed an analysis approach forthe inference of mode structures, and constructed a GAE
model. The model, which is in cylindrical geometry, accounts for a non constant density profile, features helicalripple averaging to compute the Suydam criterion and generate -ι, and includes the vacuum region. The model has
been applied to the range of -ι profiles accessible to H-1
plasmas, and the minimum in the continuum computed asa function of helical winding current ratio κh . The eigenfunction of a candidate (m, n) = (3, 4) GAE mode was determined at an angular frequency of 12.1 krad s−1 . The
eigenfunction has similar structure to recent observationsin H1, including the same number of radial nodes which
occur at similar radial positions.
In ongoing work we are developing forward modelsof Thomson scattering and charge exchange recombination
spectroscopy for inference of force balance in MAST plasmas. We will be exploiting this information not only to improve equilibrium reconstruction, but to use Bayesian inference as a tool with which to resolve competing equilibrium models. In particular, we intend to compare force balance descriptions of ideal MHD, ideal MHD with flow [1],two fluid models [27], and energetic fluid resolved equilibria [2]. The Bayesian model oﬀers a rigorous frameworkwith which to quantify the fit and thereby elucidate the important underlying physics. With regards mode structure,
we intend to transform the cylindrical eigenfunction to abeam cross-section equilibria [28] and identify a candidate
mode set which spans the diﬀerent type of mode structures
observed. The next stage of this work involves deploymentof forward models for mode structures into MINERVA tocompute the posterior for the mode class and parameters.
On a longer time scale we envisage the cylindrical plasma
model will be replaced by a fully 3D ideal MHD wave
mode code, CAS3D [29]. As with equilibrium modelling,
the importance of mode analysis by Bayesian inversion isthat it oﬀers a rigorous framework in which to identify
wave mode structures, some of which have deleterious ef-

fects on plasma performance.

Acknowledgments

This work was jointly funded by the Australian Government through International Science Linkages Grant
CG130047, the Australian National University, the United
Kingdom Engineering and Physical Sciences ResearchCouncil, and by the European Communities under the contract of Association between EURATOM and UKAEA.
The views and opinions expressed herein do not necessarily reflect those of the European Commission.
[1] L. Guazzotto, R. Betti, J. Manickam, and S. Kaye, Phys.
Plasmas 11, 604 (2004).
[2] M. J. Hole and G. Dennis, Plas. Phys. Con. Fus. 51, 035014(2009).
[3] N. J. L. Cardozo et al., Phys. Rev. Lett. 73, 256 (1994).
[4] B. P. van Milligen, I. G. J. Classen, and C. J. Barth, Rev.
Sci. Instrum. 74, 3998 (2003).
[5] N. J. L. Cardozo et al., Phys. Rev. Lett. 90, 0031 (2003).
[6] G. A. Cottrell, in Maximum Entropy in Action, edited by
B. Buck and V. A. Macaulay (Oxford Science, Oxford,1990).
[7] A. P. Millar, D. C. McDonald, and D. A. Diver, Plas. Phys.
Cont. Fus. 42, 337 (2000).
[8] R. Fischer, A. Dinklage, and E. Pasch, Plas. Phys. Cont.
Fus. 45, 1095 (2003).
[9] A. Dinklage, R. Fischer, and J. Svensson, Fus. Sci. Techn.
46, 355 (2004).
[10] J. Svensson et al., Rev. Sci. Instrum. 75, 4219 (2004).
[11] J. Svensson and A. Werner, Plas. Phys. Cont. Fus. 50,
085002 (2008).
[12] O. Ford et al., Rev. Sci. Instrum. 79, 10F324 (2008).
[13] M. Reginatto and A. Zimbal, Rev. Sci. Instrum. 79, 023505(2008).
[14] A. Sykes et al., Phys. Plas. 8, 2101 (2001).
[15] B. Lloyd and et al, Nuc. Fus. 47, S658 (2007).
[16] H. Meyer and et al, Nuc. Fus. 49, 104017 (2009).
[17] M. Hole et al., Plas.Phys. Con. Fus. 47, 581 (2005).
[18] M. J. Hole, L. C. Appel, and R. Martin, Rev. Sci. Instrum.
80, 123507 (2009).
[19] J. Besag, Journal of the Royal Statistical Society. Series B
36, 192 (1974).
[20] D.
S.
Sivia
and
J.
Skilling,
Data Analysis: A Bayesian Tutorial, 2nd ed. (Oxford
Science Publications, Oxford, 2008).
[21] J. Wesson, Tokamaks, 2nd ed. (Oxford Univ. Press, Oxford,1997).
[22] W.
H.
Press,S.
A.
Teukolsky,
W.
T.
Vetterling,
and
B.
P.
Flannery,Numerical Recipes in Fortran 77: The Art of Scientific Computing,
2nd ed. (Univ. of Cambridge Press, Cambridge, 1997).
[23] K. Appert, R. Gruber, F. Troyon, and J. Vaclavik, Plas.
Phys. 24, 1147 (1982).
[24] T. Tatsuno and M. Wakatani, Nuc. Fus. 39, 1391 (1999).
[25] R. L. Dewar et al., Phys. Rev. E 70, 0066409 (2004).
[26] D. W. Ross, G. L. Chen, and S. M. Mahajan, Phys. Fluids
25, 652 (1982).
[27] K. G. McClements and A. Thyagaraja, Mon. Not. R. Astron. Soc. 323, 733 (2001).
[28] M. G. Davidson, R. L. Dewar, H. J. Gardner, and J.
Howard, Aus. J. Phys. 48, 871 (1995).
[29] C. Schwab, Phys. Fluids B - Plas. Phys. 9, 3195 (1993).

486