D Estève

Y Sarazin
email: yanick.sarazin@cea.fr

X Garbet

V Grandgirard

S Breton

P Donnel

Y Asahi

C Bourdelle

G Dif-Pradalier

C Ehrlacher

C Emeriau

Ph Ghendrih

C Gillot

G Latu

C Passeron

Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport

Trace impurity transport is studied with the flux-driven gyrokinetic GYSELA code [V. Grandgirard et al., Comp. Phys. Commun. 207, 35 (2016)]. A reduced and linearized multi-species collision operator has been recently implemented, so that both neoclassical and turbulent transport channels can be treated self-consistently on an equal footing. In the Pfirsch-Schlüter regime likely relevant for tungsten, the standard expression of the neoclassical impurity flux is shown to be recovered from gyrokinetics with the employed collision operator. Purely neoclassical simulations of deuterium plasma with trace impurities of helium, carbon and tungsten lead to impurity diffusion coefficients, inward pinch velocities due to density peaking, and thermo-diffusion terms which quantitatively agree with neoclassical predictions and NEO simulations [E. Belli et al., Plasma Phys. Control. Fusion 54, 015015 (2012)]. The thermal screening factor appears to be less than predicted analytically in the Pfirsch-Schlüter regime, which can be detrimental to fusion performance. Finally, self-consistent nonlinear simulations have revealed that the tungsten impurity flux is not the sum of turbulent and neoclassical fluxes computed separately, as usually assumed. The synergy mostly results from the turbulence-driven in-out poloidal asymmetry of tungsten density. This result puts forward the need for self-consistent simulations of impurity transport, i.e. including both turbulence and neoclassical physics, in view of quantitative predictions for ITER.

Introduction

Impurity transport studies have recently gained a renewed interest with the perspective of using tungsten divertors in ITER [1] and WEST [2,3] tokamaks. High-Z materials such as tungsten (Z W = 74) are only partially ionized even in the plasma core despite high temperatures, so that they can already lead to prohibitive radiative losses at low concentrations [START_REF] Pütterich | Summers and the ASDEX Upgrade Team[END_REF]. These losses then impact dramatically plasma performance and stability. Actually, on-axis accumulation of tungsten has been reported in several tokamaks, including H-mode plasmas in JET [5] and ASDEX Upgrade [6]. Conversely, in Alcator-Cmod, relatively flat electron density profiles combined with peaked electron temperature in the central region appear to be efficient to prevent tungsten accumulation in the core [7]. Tungsten is not the only impurity to be present in ITER plasmas. In addition to helium ashes, extrinsic impurities will also be present, due to medium-Z (such as neon Z N e = 10 and argon Z Ar = 18) impurity seeding aiming at radiating the power at the edge, and low-Z beryllium (Z Be = 4) due to plasma wall interactions.

Impurity transport is actually a complex issue, even in the relevant limit of trace impurities where their strength parameter is negligible (α = n Z Z 2 /n i Z 2 i 1, with Z the impurity and i the main ion species). The radial impurity flux contains both a diffusive and a convective contribution: Γ Z = -D∇n Z + V n Z . In the absence of transport barriers, in L-mode plasmas, impurity transport is usually found to be dominated by turbulence, at least in the gradient region. Indeed, the experimentally deduced transport coefficients exceed the ones predicted by the neoclassical theory by factors or even orders of magnitudes in tokamak plasmas [8,9,10,11] and spherical tokamaks [12,13]. Conversely, in the very core region [14,15,16] and at the position of transport barriers, either at the edge or at internal transport barriers, impurity transport is often found consistent with neoclassical predictions [17,18,19]. In such regimes, the impurity confinement time sometimes appears to scale with the edge density gradient of the main ion species, as one would expect from neoclassical transport theory [20]. Should neoclassical transport dominate tungsten transport in ITER H-mode edge plasmas, tungsten is then expected to exhibit hollow density profiles in the pedestal region [21].

In addition, it has been recently realized that standard neoclassical predictions for impurity transport can be substantially modified in the highly collisional regime (Pfirsch-Schlüter), relevant for high-Z impurities, in the case where their density profile exhibits in-out asymmetry [22,23,24,25]. Such poloidal asymmetries can emerge as a result of centrifugal force [26] or of anisotropic (i.e. such that T ⊥ = T ) radio-frequency heating, especially Ion Cyclotron Resonance Heating [27,28,29].

All in all, it appears that the turbulent and neoclassical coefficients for impurity transport, both diffusivity and convection velocity, can be of the same order of magnitude, especially for high-Z impurities such as tungsten. At present, first principles simulations of these transports are performed with different dedicated codes, explicitly assuming that both transport channels are additive (see e.g. [7]). One of the key questions is whether this assumption is valid, or whether neoclassical and turbulent transports exhibit synergetic effects which would modify our understanding and predictions of impurity concentration in tokamaks.

In this framework, a multi-species collision operator has been developed for fullf gyrokinetic codes and successfully implemented [30] in the flux-driven GYSELA code [31]. In this paper, we first show analytically that this operator allows one to recover neoclassical predictions of impurity transport in the Pfirsch-Schlüter regime (section 2). The numerical benchmark of the code with respect to neoclassical theory for trace impurity transport is then successfully performed in section 3 in all three regimes of collisionality (banana, plateau and Pfirsch-Schlüter) by considering different types of impurities. Last but not least, section 4 presents the first evidence of synergy between turbulent and neoclassical transport of tungsten. Turbulence driven poloidal asymmetries are suspected to play a key role in this synergy.

Recovering neoclassical fluxes from gyrokinetics

Radial flux and toroidal momentum conservation

The spatial coordinates are (ψ, θ, ϕ) where ψ is the minus poloidal magnetic flux normalized to 2π, θ is the poloidal angle and ϕ the toroidal angle. The Jacobian of this metric is 1/B • ∇θ where the unperturbed magnetic field is

B = I∇ϕ + ∇ϕ × ∇ψ (1)
We consider the gyrocenter distribution function F (z) where z = (ψ, θ, ϕ, v , µ) are the gyrocenter coordinates. Here v is the gyrocenter parallel velocity and µ the magnetic moment. To derive the conservation laws, it is convenient to use the gyrokinetic Fokker-Planck equation [32,33]

∂F ∂t + 1 B * || ∂ ∂z • żB * || F = C(F ) (2)
where ż = dz dt and

B * || = B + mv e b • (∇ × b) (3)
is the Jacobian of the gyrocenter transformation (m is the mass, e the charge and b = B B the unit vector along the magnetic field). For an electrostatic case, the gyrocenter equations of motion are

B * || dx dt = -B * || [H, x] = v B * + 1 e b × ∇ µB + e φ (4)
B * || m dv dt = -B * || [H, mv ] = -B * • ∇ µB + e φ (5)
where

B * = B + mv e ∇ × b (6)
and

H = 1 2 mv 2 + µB + e φ ( 7
)
is the Hamiltonian, φ = J • φ are the gyro-average electric potential. It was shown by several authors [34,35,36,37] that the toroidal momentum conservation equation reads

∂ t L ϕ + ∂ ψ Π ψ ϕ = e Γ ψ -Γ ψ E + I B R ψ (8)
Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport 4

where

L ϕ = m dτ u ϕ F (9) Π ψ ϕ = m dτ F u ϕ v ψ (10)
Γ

ψ E = dτ F ∂ ϕ φ (11) Γ ψ = dτ v ψ F = Γ ⊥ • ∇ψ ψ ( 12
)
where R is the collisional drag force

R = d 3 vmv C(F ) ( 13
)
with

d 3 v = 2π m B * || dµdv . Γ ⊥ is the perpendicular flux. Here u ϕ = I B v is the toroidal covariant component of the gyrocenter velocity (dimensionally: [u ϕ ] = L 2 .T -1 ), while v ψ = v G •
∇ψ is the radial projection (or the radial contra-variant component) of the gyrocenter velocity v G . The phase space volume element in between two magnetic surfaces ψ and ψ + dψ is

dτ = d 3 v dθdϕ B • ∇θ (14)
The bracket ... ψ indicates a flux surface average

... ψ = dθdϕ B•∇θ ... dθdϕ B•∇θ (15)
The collision friction force (last term in the r.h.s. of Eq.( 8)) did not appear in earlier derivations of the momentum conservation equation, either because the collisionless problem was treated, or because an unique species was considered. It is however easy to add. The flux Γ ψ E is the radial flux due to the E × B drift velocity. The sum over all species is a polarization flux, which can be neglected for now. The divergence of the Reynolds stress ∂ ψ Π ψ ϕ contains a turbulent contribution, which is expected to be ρ * smaller than the flux in the gyrokinetic ordering k R o(ρ * ). The neoclassical contribution to Π ψ ϕ is predicted to be even smaller. In steady-state, the flux average of Eq.( 8) then gives

Γ ψ = Γ ψ E + Γ ψ coll (16) Γ ψ coll = - I e R B ψ (17)
This equation is often recast as

Γ ψ coll = Γ ψ BP + Γ ψ P S (18)
where

Γ ψ P S = - I e R 1 B - B B 2 ψ ψ (19)
Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport 5 and

Γ ψ BP = - I e R B ψ B 2 ψ ( 20
)
are respectively the Pfirsch-Schlüter and Banana-Plateau fluxes. This splitting anticipates the fact that, in the Pfirsch-Schlüter regime, the collisional drag force R is governed by the compressibility of both parallel velocity and parallel heat flux (cf. eqs.( 35),( 37)) due to the inhomogeneity of the magnetic field. In absence of turbulence, the flux Γ ψ E due to the E × B drift velocity is small since the poloidal asymmetries of the electric potential are small. The flux Γ ψ coll then coincides with the conventional definition of the neoclassical flux.

Collision operator

We address now the multi-species case. The collision operator for a species "a" that is implemented in the GYSELA code [31] is a simplified one [30], namely C a (F a ) = b C ab (F a ) where the contribution from collisions species "a" colliding with species "b" is

C ab (F a ) = 1 B * || ∂ ∂v B * || D d,ab F M 0a ∂ ∂v f a - m a v U d,a T a -ν s,ab m a v T a U d,a -U ba F M 0a + 2 3
q ab N a T a m a v 2 2T a - 3 2 F M 0a (21)
where f a = Fa F M 0a . Expression of the deflection diffusion D d,ab and the collisional friction rate ν s,ab are given in Appendix A. The distribution F M 0a is the unshifted Maxwellian built with the density N a and the temperature T a

F M 0a (x, v, t) = N a (x, t) m a 2πT a (x, t) 3/2 exp - E T a (x, t) (22)
The single species restoring momentum coefficient is

m a v T a U d,a (v) = 3 4π dΩξf a (v, ξ, γ)) (23)
where ξ = v /v and dΩ = 2πdξ is the element of solid angle in the velocity space. The inter-species momentum restoring coefficients are

U ab = ν s,ab v f a a ν s,ab mav 2 Ta a (24)
where

... a = 1 N a 2π B * || B v ⊥ dv ⊥ dv F M 0a ... = 2 √ π ∞ 0 dx a x 2 a e -x 2 a +1 -1 dξ... (25)
The energy restoring coefficient is

q ab = N a T a σ ab f a a + Q M ab ( 26
)
where

Q M ab = -3N a ma ma+m b ν ab (T a -T b )
is the energy exchange between two centered Maxwellians and σ ab (v) is given by

σ ab (v) = -ν d,ab m a v 2 2T a v ∂ ∂v ln D d,ab F M 0a v (27)
Standard GYSELA simulations use a limited number of mesh points in the µ space (typically 16). Therefore, so as to minimize the error when computing the friction force, the distribution function is projected on the two first Sonine polynomials (generalized Laguerre L (m) j polynomials of order m = 3/2) in the velocity modulus space (cf. reference [36], appendix D). The friction part of the field particle collision operator (2nd line of Eq.( 21)) then reads

C s,ab (F a ) = -ν s,ab m a vξ T a F M 0a V a -V b + 2 5 m a v 2 2T a - 5 2
q a N a T a + 3 5 1 1 + x 2 ba q b N b T b (28)
where V a and q a are the parallel velocity and heat flux. The friction force can then be written as

R ab = -N a m a ν ab V a -V b - 3 5 1 1 + x 2 ab q a N a T a + 3 5 1 1 + x 2 ba q b N b T b (29)
where

x ab = v T b v T a .

Neoclassical impurity flux

We consider now the case of two ion species only: the main ion "i" and a heavy impurity "Z", which is supposed to be in a collisional (Pfirsch-Schlüter) regime. We consider here a neoclassical case, i.e. we ignore the flux due to the E ×B drift velocity in Eq.( 16). The impurity flux is calculated by using the relationship between the flux and the parallel collisional drag force Eq.( 17)

Γ Z • ∇ψ ψ = - I Ze R Zi B ψ (30)
Since

x Zi = v T i /v T Z
1, the contribution from the parallel impurity flux in R Zi Eq.( 29) is negligible so that

R Zi = -N Z m Z ν Zi V Z -V i + 3 5
q i N i T i (31)
In principle, the parallel flows should be calculated by solving a kinetic equation for each species. The calculation is demanding for ions in weakly collisional regime because of the boundary layer near the transition from passing to trapped particles [START_REF] Hirshman | [END_REF]43,23,24]. The complete calculation is delicate and we give here an approximate derivation, which allows recovering the Hirshman expression of the thermal screening coefficient. For passing particles, an acceptable solution of the Fokker-Planck equation is

F a = F M a + δF a ,
where δF a is of the form

δF a = F M a m a v T a,eq W a B - IT a,eq e ∂ ln F M 0a ∂ψ 1 B - B B 2 ψ ( 32
)
where

W a = V a B ψ B 2 ψ + 2 5 E T a,eq - 5 2 1 P a,eq q a B ψ B 2 ψ (33)
and

T a,eq ∂ ln F M 0a ∂ψ = 1 N a,eq ∂P a,eq ∂ψ + e a ∂Φ eq ∂ψ + E T eq - 5 2
∂T a,eq ∂ψ

The approximate distribution function Eq.( 32) yields the following expressions of the parallel velocities and heat fluxes

V a = BV a ψ B B 2 ψ -Ω a I 1 B - B B 2 ψ ( 35
)
where Ω a (ψ) = 1 N a,eq e a ∂P a,eq ∂ψ + ∂Φ eq ∂ψ

and

q a = Bq a ψ B B 2 ψ - 5 2 P a,eq 1 e a ∂T a,eq ∂ψ I 1 B - B B 2 ψ ( 37
)
where N a,eq , T a,eq and P a,eq are the unperturbed density, temperature and pressure, Φ eq the unperturbed electric potential. These expressions can then be plugged into Eq. (31).

However one problem remains: the flux surface average of the parallel velocity and heat flux are unknown. Since heavy impurities are supposed to be in the fluid regime (i.e. such that the pressure anisotropy can be neglected: ∇ : Π ∇P , with Π the pressure tensor and P the scalar pressure), the parallel force balance equation reads

R Zi = ∇ P Z + N Z,eq Ze∇ Φ ( 38
)
which implies BR Zi ψ = 0 [START_REF] Hinton | Theory of plasma transport in toroidal confinement systems[END_REF] and therefore

BV Z ψ = BV i ψ - 3 5 Bq i ψ (40)
Hirshman et al. [START_REF] Hirshman | [END_REF] showed that

q B ψ 1/2 q * B 0 (41)
where q * = q N eq T eq Teq eBLp is the parallel diamagnetic heat flux . So q B ψ can be neglected in the limit of small inverse aspect ratio 1. This leads to the well known constraint of equal average velocities of ions and impurities

BV Z ψ = BV i ψ (42)
This is actually the weak point of the demonstration. The limit of large aspect ratio is indeed a poor one, and never fulfilled in actual numerical simulations. We nevertheless pursue the calculation in that framework to find an explicit expression of the friction force (accounting for the fact that ∂ ψ ln T Z,eq Z∂ ψ ln T i,eq at equal temperature)

R Zi = ν Zi m Z T Z,eq Ze I 1 B - B B 2 ψ N Z,eq ∂ ln N Z,eq ∂ψ -Z ∂ ln N i,eq ∂ψ + Z 2 ∂ ln T i,eq ∂ψ (43)
and consequently the impurity flux

Γ Z • ∇ψ ψ = -ν Zi m Z T Z,eq Z 2 e 2 I 2 1 B 2 ψ - 1 B 2 ψ N Z,eq ∂ ln N Z,eq ∂ψ -Z ∂ ln N i,eq ∂ψ + Z 2 ∂ ln T i,eq ∂ψ (44)
In circular concentric geometry, large aspect ratio, dψ dr = r q and

I 2 1 B 2 ψ - 1 B 2 ψ 2r 2 (45)
This formula then reads

Γ neo Z = Γ Z • ∇r ψ = -N Z,eq D neo z ∂ ln N Z,eq ∂r -Z z ∂ ln N i,eq ∂r -Z z H neo z ∂ ln T i,eq ∂r (46) The neoclassical diffusion coefficient is D neo z = 2q 2 ν Zi ρ 2 Z in the Pfirsch-Sclüter regime, where ν iZ = √ 2ν i N Z,eq Z 2 N i,eq
and ρ Z = m Z v T Z ZeB 0 . We recall that N Z,eq m Z ν Zi = N i,eq m i ν iZ . The last two terms of eq. ( 46) are the main off-diagonal terms of the transport matrix. The pinch velocity often refers to

V pinch z = D neo z Z z ∂ r (ln N i,eq
). The thermal screening factor H neo z is expected to be equal to -1/2 in the Pfirsch-Schlüter regime [40,47]. This expression eq. ( 46) agrees with Hirshman-Sigmar [40,47]. It must be noted that the rule BV Z ψ = BV i ψ is not very accurate without external sources of momentum. Fülöp and Helander [23] showed that diamagnetic corrections are important. It turns out that the Hirshman-Sigmar result still holds in the case where impurity poloidal asymmetries are small. When they are accounted for, large differences are found in the coefficients (diffusion and screening factor) of the impurity flux [22,23,24].

Neoclassical transport in the 3 collisionality regimes

The objective of this section is to benchmark the implemented collision operator against the predictions of neoclassical theory regarding the impurity flux, eq. ( 46). Three neoclassical predictions will be checked hereafter:

• the dependence of D neo z with respect to the collisionality ν * z , • the dependence of the pinch velocity V pinch z with respect to the impurity effective charge Z,

• the thermal screening factor H neo z in the Pfirsch-Schlüter regime.

To this end, purely neoclassical impurity transport is studied with the GYSELA code by performing axisymmetric simulations. This means that all non-axisymmetric modes (n = 0, with n the toroidal mode number) are filtered out at each time step, thus preventing any turbulence to develop, whatever the temperature and density gradients. The main ion species is deuterium. Three different impurities are considered: helium (A He = 4), carbon (A C = 12) and tungsten (A W = 184). A single ionization state is considered for each of them: helium and carbon are assumed fully ionized (Z He = 2, Z C = 6), and tungsten only partially (Z W = 40). The collisionality of deuterium is chosen equal to ν * D = 0.1 at the center of the simulation domain, denoted r p . Assuming equal temperatures, the collisionality of the various impurities scales like

ν * z ≈ √ 2 (Z z /Z D ) 2 (A D /A z ) 1/2 ν * D = 2ν * D Z 2 z /A 1/2
z , with z = {He, C, W }. The collisionalities at mid-radius are given in Table 1. With this choice of ν * D , it appears that they cover the entire range of the collisionality regimes: while deuterium and helium are in the banana regime, carbon is in the plateau regime, and tungsten in the Pfirsch-Schlüter regime. The collisionality profiles are plotted on Fig. 1. These impurities are kept at a trace level, i.e. such that the ratio α = (n z Z 2 z /n D Z 2 D ) is much smaller than unity. Taking this ratio of the order of α ≈ 10 -3 leads to the concentrations listed in Table 1.

A simulation with deuterium only (without any impurity) is first run on about one collision time so as to reach neoclassical equilibrium. It is then restarted with the addition of the chosen impurity. These simulations are performed at the reference (i.e. for hydrogen) normalized gyroradius value equal to ρ * 0 = ρ H /a = 1/150 ≈ 6.67 10 -3 . Assuming equal temperatures, the ρ * values of deuterium and of the chosen impurities then scale like ρ * z = (A 1/2 z /Z z ) ρ * 0 . The corresponding values are given in Table 1. The aspect ratio is 3.2, and the safety factor typically ranges from 1.5 to 2.5. The number of grid points is the following: (N r , N θ , N ϕ ) = (256, 256, 32), N v = 128 → 256, and N µ = 16 → 32. Note that N ϕ is increased up to 128 in the case of turbulent simulations, discussed in section 4. Table 1. Collisionality, concentration and ρ * value of the 4 species of the various simulations discussed in the paper. The values are given at the center of the radial domain.

Diffusion and pinch velocity

The first component of the impurity flux to be studied is the diffusion coefficient D z . The only non-vanishing gradient is the impurity density gradient ∂ r N Z,eq . A single simulation would typically require to run on a particle confinement time to allow one to extract accurately the D z value. So as to save computation time, three simulations have been run at three different values of R/L nz = {0.7; 1.4; 2.1}, measured at r p (with 1/L nz = -∂ r ln N Z,eq ). The profiles are plotted on Fig. 2 (left). Once the impurity has been injected, and after a short reorganization time, the impurity flux reaches a well defined value, which slowly decreases in time with R/L nz . The value recorded at r p is plotted on Fig. 2 (right) for the 3 impurities and the 3 gradients. The impurity flux contains both contributions of the E × B drift and of the magnetic (curvature and ∇B) drift. Each series of points defines a line that passes by zero, as it should. The slope, equal to D z , is obtained by a least squares fit.

On Fig. 3, the D z values obtained with GYSELA in this way are successfully compared with the ones predicted by Hirshman-Sigmar [45] and the ones obtained with the NEO code [48] (adiabatic electrons, Hirshman-Sigmar type of collision operator).

So as to estimate the pinch velocity, a non-vanishing value of the deuterium density gradient is considered in addition to ∂ r N Z,eq . The deuterium density profile is homothetic to the one of impurities, with R/L n D = 0.7 at r p . The steady state impurity density profiles would be obtained when the flux vanishes (points on the x-axis on Fig. 4). In this case, they can reveal extremely peaked, as a result of ambipolarity. Indeed, they scale as N Z,eq /N Z,eq (r = a) = (N i,eq /N i,eq (r = a)) Zz . In the case of tungsten (Z W = 40), such a strongly peaked profile would require prohibitive numerical resources, namely large number of radial grid points. Also, too short density gradient lengths L nz can conflict with one of the gyrokinetic assumptions, stating that ρ z /L nz 1. Conversely, using the same procedure as the one used for the diffusion coefficient does not suffer such difficulties. Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport 13

Thermal screening in the P-S regime

When the main ion species exhibits a non-vanishing temperature gradient, there is an additional contribution to the convective part of the flux, due to thermodiffusion. It is expressed in terms of the neoclassical thermal screening factor H neo z (cf. eq. ( 46)):

V thermo z = Z z D neo z H neo z ∇ r log T .
The associated velocity is expected to be outwards (positive) in the Pfirsch-Schlüter regime, for which the value H neo z = -1/2 is commonly admitted (limit of large aspect ratio, provided ν * is not too large [48]). Should tungsten transport be governed by neoclassical physics in the Pfirsch-Schlüter regime, this term could reveal essential to prevent its accumulation in the core of tokamak plasmas, where its strong radiation critically limits fusion performance. Evidence of impurity outward flow due to thermal screening has already been reported experimentally for low-to medium-Z (helium, nitrogen, neon) materials [49]. Some of the observed outflow of impurities when applying Electron Cyclotron Resonant Heating could well proceed from the same physics [50].

A deuterium-tungsten simulation has been performed with equal initial density and temperature gradient lengths: R/L nD = R/L nW = 0.1 and R/L T D = R/L T W = 1. Also, temperatures are the same. In this case, the tungsten neoclassical flux is mostly governed by the density and temperature gradients of deuterium only, because they are weighted by the charge of tungsten Z W = 40 (cf. eq. ( 46)). The following procedure is used to retrieve the thermal screening factor H neo W from the GYSELA simulation. The neoclassical diffusion coefficient is assumed to be equal to the Hirshman-Sigmar theoretical prediction [45]. Given the actual density and temperature gradients, H neo W is then such that the right hand side of eq. ( 46) matches the measured tungsten flux. Notice that the finite discrepancy between D neo W found in GYSELA simulations and the one predicted by the theory (cf. Fig. 3) impacts the evaluation of H neo W with the adopted method. The theoretical screening factor integrates the contributions from the various regimes, namely "banana-plateau" (BP) and Pfirsch-Schlüter (PS):

H W = D BP W H BP W + D P S W H P S W D BP W + D P S W (47)
It is plotted on Fig. 6 (red line), where it appears to be dominated by the PS contribution. Its slight departure from -1/2 is mainly due to the finite value of ( 0.13 at r p ). At r/a = 0.45, it is compared a with 2 NEO simulations using two different versions of the collision operator (cf. Tab. 2): the most advanced one (Fokker-Planck), and the one retaining only the pitch-angle part. The latter one corresponds to the limit in which the analytical predictions have been made. It already appears that there is a mismatch between both (pitch-angle NEO and analytic theory), of the order of 18%. Part of the explanation likely comes from the → 0 limit which is considered in the analytical derivation: corrections of order 1/2 , which is equal to √ 0.13 ∼ 0.36 in the present case, would be expected. Also, the F-P operator leads to a less effective screening than analytically predicted. This is detrimental to fusion performance: temperature gradients might well be less efficient in screening the natural peaking effect of high Z impurities.

Overall, as expected, the screening factor obtained with GYSELA is negative, i.e. leads effectively to an outward transport of tungsten. It has also the right order of magnitude, although it is about 34% smaller than the NEO prediction with the F-P operator. As already stated, part of this discrepancy is due to the uncertainty regarding the diffusion coefficient D neo W .

GYSELA NEO (Fokker-Planck) NEO (pitch-angle) Theory (H-S) H W value:

-0.211 -0.321 -0.377 -0.461

Table 2. Thermal screening factor H W for the D-W case at r/a = 0.45. We review hereafter some of the foundations of the neoclassical theory in the P-S regime, on which rely the derivation of H neo z , among others. One of the important results of neoclassical theory in the P-S regime is that the parallel velocities and heat fluxes exhibit a poloidal left-right asymmetry due to the magnetic field inhomogeneity. This asymmetry, due the inhomogeneity of the magnetic field, is responsible for their parallel compressibility. Besides, the asymmetry scales with the pressure and temperature gradients, as evident in eqs. ( 35), (37). The poloidally fluctuating part of the parallel deuterium (left) and tungsten (middle) velocities are plotted on Fig. 7, as well as that of the deuterium parallel flux (right). The top row shows the instantaneous 2-dimensional values of these fields at the end of the simulation. The bottom row corresponds to the expected asymmetry, proportional to (1/B -B/ B 2 ψ ) and weighted by the radial gradients. The agreement between the modeling (top row) and theoretical results (bottom row) looks satisfactory, both in terms of magnitude and in shape. With the adopted collision operator, parallel flows in GYSELA turn out to exhibit the right structure. In summary, this set of purely neoclassical simulations -where the sole axisymmetric modes n = 0 have been retained -has shown that the reduced collision operator implemented in GYSELA, eq. ( 21), allows one to recover the main results of the neoclassical theory regarding impurity transport in all three collisionality regimes, banana, plateau and Pfirsch-Schlüter. This includes the diffusion coefficient, the pinch velocity due to the main ion density gradient, and the thermal screening factor due to the main ion temperature gradient. This latter term has been found to slightly depart from the analytical estimate, which is close to H Z = -1/2 in the Pfirsch-Schlüter regime. The reason is twofold. Firstly, numerical simulations with the dedicated neoclassical code NEO also predict a negative value for H Z , but closer to zero. This strongly suggests that the analytical estimate is probably too optimistic. Secondly, the uncertainties in the determination of the diffusion coefficient in GYSELA do affect the computed value of H Z . The next section explores the interplay between neoclassical and turbulent channels for impurity transport.

Evidence of synergy between neoclassical and turbulent impurity transport

In tokamak plasmas, the radial transport of matter and energy, governed by collisions and turbulence, is intrinsically multi scales. While neoclassical transport results from stationary large scale structures, namely static (m, n) = (1, 0) modes (m, n being the poloidal and toroidal Fourier wave numbers, respectively), turbulence develops fluctuating small scale modes m, n 1. On the basis of this scale separation, it is usually assumed that both contributions are additive. In turn, these two transport channels are modeled with different dedicated codes (see e.g. [5,6,7,52,12]). One of the key questions is whether this assumption is valid, or whether neoclassical and turbulent transports exhibit synergistic effects.

This section gives the first evidence of synergistic effects between neoclassical physics and ITG (Ion Temperature Gradient) turbulence with respect to impurity transport. Previous evidences of such a synergy have been reported in the literature. They deal with poloidal flow [51], heat transport [53,54] and parallel momentum transport [36,55]. Also, drift wave turbulence has been shown analytically to be able to significantly contribute to the bootstrap current through direct modification of the axisymmetric plasma current distribution via turbulence-induced electron detrapping [56]. This mechanism could reveal particularly vigorous during transient bursts of turbulent transport. In all of these studies, identifying any synergy between turbulent and neoclassical transports depends on the definition which is given to both transport channels. Obviously, in experiments -like in self consistent simulations including both turbulence and neoclassical physics -the total flux is simply the sum of turbulence and neoclassical contributions. The question is then the following: can this self-consistent flux be recovered by adding up the fluxes coming from two independent simulations, one dealing with turbulence physics only, the other one with neoclassical physics only? If the answer is "No", this means that both channels then exhibit some kind of synergy, which would remain being understood. The adopted GYSELA set up which led to the conclusion that there exists synergy between these transport channels is detailed in section 4.2.

Impurity transport in collisional turbulence

The transport of 3 impurities (He, N e and W ) is studied in the ITG turbulent regime. They find themselves in the three regimes of neoclassical transport, namely banana, plateau and Pfirsch-Schlüter, respectively. The total impurity flux is then the sum of two contributions, carried by the electric and magnetic drifts:

Γ tot z ψ = d 3 v F z (v D,z + v E,z ) • ∇r ψ (48)
Here, the magnetic v D,z and electric v E,z drifts are defined by:

v D,z = m z v 2 G + µ z B ZeB * || b × ∇B B ( 49
) v E,z = b × ∇(J z .φ) B * || (50)
with b = B/B and J z the gyro-average operator. There is not a one-to-one correspondence between these two contributions, those of the magnetic and electric drifts on the one hand, and the definitions of the neoclassical and turbulent fluxes on the other hand. Indeed, part of the neoclassical flux is carried by the electric drift, more precisely by its axi-symmetric component n = 0. This point will be taken into account when looking for possible synergies between neoclassical and turbulent transports, in section 4.2.

The simulations are performed as follows. They consist in two consecutive phases. First, GYSELA is run in the single-species mode, evolving the distribution function of the main ion (deuterium) only, until the saturated regime of turbulence is reached. This phase typically lasts several tens of turbulence auto-correlation times. In a second phase, the simulation is restarted in the two-species mode, where the full distribution functions of deuterium f D and impurity f z are self-consistently advanced in time. In the present cases, there is no source of particles. Impurities are initialized with a density profile n z homothetic to the one of deuterium n D , but at low concentration (c z = n z /n e0 1, with n e0 the constant-in-time electron density) so that the impurity remains at a trace level. Although this is not critical for trace impurities, the deuterium concentration is reduced in the quasi-neutrality equation, so that the sum (c D + Zc z ) remains equal to one. By doing so, one ensures that no charge is injected in the system. The flux-surface average of the radial impurity flux is plotted in Fig. 8 for each of the 3 simulations. It is time-averaged on several turbulence correlation times during the turbulence saturation phase. Two main observations can be made. First, the fluxes carried by the electric and magnetic drifts are of the same order of magnitude, and tend to compensate each other. This trend seems to be more pronounced for the lighter impurity, helium. Such a behavior has also been reported for parallel momentum transport [36,55], although no clear explanation has been proposed so far. Second, this compensation is however not complete so that, for this set of parameters, the total flux is outward for helium, and inward for neon and tungsten. This latter observation is potentially detrimental to fusion performance in tungsten divertor machines like WEST and ITER, but a more systematic study as function of plasma parameters is clearly needed. This is left for future work. Interestingly, as evident on Fig. 9, the magnetic (a) and electric (b) drift components exhibit a rich dynamics: the former due to the turbulent fluctuations of the electric potential, the latter due to the dynamics of both the corrugated equilibrium gradients and of the m = 1 modes (see below). As a result, the total flux Fig. 9(c) features complex dynamical patterns, with large deviations from the mean: small scale avalanche-like events can even lead, locally and temporarily, to the sign reversal of the flux.

Evidence of synergies

The level -and existence -of synergy between neoclassical and turbulent transports largely relies on the definition of these two terms. Obviously, in experiments as well as in self-consistent simulations, the total flux is the sum of neoclassical and turbulent contributions. Then how to capture any possible synergy? This question makes sense when realizing that present estimates of these two transports are performed with different dedicated numerical codes. Then, our strategy aims at assessing whether, and in which regimes, the total flux can effectively be simply obtained by adding up the numerical calculations which are currently performed in the literature. GYSELA allows one to run reduced sets of simulations, either retaining the sole neoclassical contribution or the turbulent one, to be compared to self-consistent simulations where neoclassical and turbulent transports are treated on an equal footing.

In that spirit, for each of the 3 considered trace impurities He, N e and W , 3 simulations have been performed: (i) Purely neoclassical (called "neoclassical" hereafter): in this case, all nonaxisymmetric toroidal modes are filtered out at each time step (i.e. all Fourier modes with n = 0 are set to zero, with n the toroidal mode number). The neoclassical flux is then the sum of both E ×B and magnetic drift contributions (the need to account for both terms is consistent with neoclassical theory, as discussed in [57,58]):

Γ neo z ψ = d 3 v F z (v D,z + v n=0 E,z ) • ∇r) ψ (51)
This definition departs from the more common understanding where neoclassical transport is carried by the magnetic drift v D,z only.

(ii) Mainly turbulent (called "turbulent" hereafter): in this case, single-species collisions only are retained (namely ν ii and ν zz ) , so that momentum or energy exchange between species is not taken into account. The fact that ν Zi = 0 in this simulation implies that the parallel friction force is also vanishing: R Zi = 0. As a result, in virtue of eq. ( 30), the radial neoclassical flux is exactly zero. Retaining intra-species collisions is however important and sufficient to account, among others, for the collisional damping of zonal flows, which are known to contribute efficiently to turbulence saturation. In the case of trace impurities, the collision operator for the main ions only is important. Indeed, impurity-impurity collisions are negligible in this case. The turbulent flux is then governed by the electric drift:

Γ turb z ψ = d 3 v F z (v E,z • ∇r) ψ (52)
(iii) Full (called "total" hereafter): in this case, no simplification is made to the collision operator, nor any filtering applied to the electric potential. More precisely, all terms of the collision operator are retained, involving both intra-and inter-species collisions, in the turbulent regime.

As an example, the case of tungsten is considered hereafter. The same trend is observed for the two other impurities. In Fig. 10, the total tungsten flux (red) from the self-consistent simulation is compared to the sum of the neoclassical and turbulent fluxes (black), coming from the reduced simulations. First, conversely to the magnetic and electric drift components of the flux, the neoclassical and turbulent contributions do not exhibit any clear sign of partial compensation. Second, the total flux of tungsten turns out to be inward in this case. However, this result should not be considered as a general statement. A thorough exploration of the parameter space (including density and temperature gradients, collisionality and the source of tungsten) should be performed before reaching more general conclusions. This is left for future work. Third and most importantly, it appears that the total flux differs from the sum of the neoclassical and turbulent fluxes, by more than a factor two at some radial locations. This result shows that neoclassical and turbulent impurity transports are not additive, in the sense that they cannot be computed separately and then simply added up.

One of the explanations of this synergy comes from the impact of turbulence on the magnitude and radial shape of poloidal asymmetries, i.e. of the m = 1 Fourier mode. Indeed, in/out asymmetries of the impurity density profile are able to greatly modify neoclassical coefficients, up to one or two orders of magnitude, as predicted theoretically [22,23,24,25]. The time-averaged (m, n) = (1, 0) component of the electric potential (actually, this is not exactly the n = 0 component: the Fourier transform is performed at a given toroidal location. However, the m = 1 spectrum exhibits a large peak at n = 0, so that the axisymmetric mode is dominant) is plotted on Fig. 11, for the three simulations with tungsten: neoclassical, turbulent and total. The magnitude of the (1, 0) mode is larger in the presence of turbulence, and exhibits smaller radial scales (b,c) than in the neoclassical case (a). Its magnitude turns out to be of the order of ρ * = 1/150.

In these simulations, the parallel Mach number of the main ion M D is fairly small, so that the one of tungsten does not exceed 10% (M W ≈ (A W /A D ) 1/2 M D ). In this case, centrifugal effects are sub-dominant. One thus expects tungsten density fluctuations to have an adiabatic response, namely δN W / N W F S ≡ N W / N W F S -1 ∼ -Z(e δφ/T Z ) [6]. As a matter of fact, this is exactly the order of magnitude which is obtained in this simulation, as evident on Fig. 12(b). There, the large scale m = 1 mode emerges from the small scale turbulent structures. Interestingly, δN W peaks at about 4.5 10 -8 , while the flux-surface averaged tungsten density (not shown here) is about 1.1 10 -7 . The ratio δN W / N W F S is well in the expected range Z(e δφ/T Z ) ≈ 0.4. Further notice that the poloidal asymmetry is larger in the turbulent regime than in the purely neoclassical simulation, Fig. 12(a), the extrema being larger by 40% to 80%.

In summary, the large scale poloidal structures of the impurity density, which are generated by turbulence (via the turbulence-induced m = 1 mode of the electric potential), add-up and/or compete with those governed by neoclassical physics, so that neoclassical and turbulent transports are not additive (in the sense already defined). A more systematic exploration of the parameter space, especially in terms of temperature and density gradients, as well as the charge and mass of the impurity, is currently performed to assess the conditions for and magnitude of this synergy. This will be the subject of a forthcoming paper, including a more detailed analysis of the mechanism for this synergy (including possible kinetic effects).

Conclusions

Impurity transport is likely to play a key role in limiting ITER performance. ITER plasmas will face intrinsic impurities with helium ashes, and extrinsic impurity influxes due to either edge impurity seeding (neon, argon) which aim at radiating most of the power or due to plasma-wall interaction (tungsten, beryllium). The case of tungsten is particularly severe, since core accumulation leads to strong radiation that can prevent the access to high performance discharges, and can even lead to radiative collapses. Quantitative predictions of impurity transport is therefore critical, as well as the search for possible control mechanisms to prevent too large impurity concentration in the plasma core.

Both collisional (neoclassical) and turbulent transport contribute to the impurity flux. Their respective magnitude depends mainly on the mass and charge state of the impurity, on the main ion density and temperature gradients, on the level (with or without transport barrier) and nature (e.g. ITG versus TEM) of turbulence, and on the heating scheme (isotropic or not, with or without torque injection). It is suspected that both transport channels could exhibit similar magnitudes in ITER plasmas for high-Z impurities like tungsten.

The flux-driven gyrokinetic GYSELA code is used to study impurity transport without any scale separation assumption between neoclassical and turbulent transport channels. Conversely to the usual hypothesis that they are uncorrelated and additive, both are treated on an equal footing in GYSELA. The implemented multi-species collision operator has been shown analytically to allow one to recover neoclassical predictions for impurity transport in the Pfirsch-Schlüter regime. The impurity flux of purely neoclassical simulations (performed by filtering out non-axisymmetric modes) of deuterium plasmas with trace impurities of helium, carbon and tungsten has been successfully benchmarked with respect to both neoclassical theory and NEO simulations in the three collisionality regimes. Especially, the neoclassical diffusion coefficient, the inward pinch velocity due to density peaking, and the thermo-diffusive contribution have been recovered. Interestingly, it appears that the thermal screening coefficient departs from the analytical estimate of -1/2 in the Pfirsch-Schlüter regime. The computed screening efficiency is less, which is detrimental to fusion performance. Finally, first nonlinear simulations in the turbulent regime show that the self-consistently computed impurity flux is not the sum of turbulent and neoclassical contributions computed independently, as usually assumed in the literature. The mismatch can reach up to a factor of two at some radial locations. The discrepancy comes from the excitation of (m, n) = (±1, 0) modes of the electric potential by turbulence, resulting in turbulencedriven in-out poloidal asymmetries of the impurity density profile. Forthcoming simulations will explore and quantify this synergy in a more systematic way by scanning the parameter space. Such results put already forward the need for self-consistent simulations of impurity transport in view of providing quantitative predictions for ITER. necessarily reflect those of the European Commission.

Appendix A. Collision operator: details

The diffusion coefficients in the velocity D v,ab is related to the deflection collision frequencies ν d,ab via the relation

D d,ab (v) = 1 2 ν d,ab (v)v 2 = v 2 T a x 2 a ν d,ab (x a ) (A.1)
where v is the velocity modulus (i.e. E = 1 2 m a v 2 ) and x a is a normalized velocity modulus

x a = E T a = v √ 2v T a (A.2)
We introduce a normalizing self-collision frequency ν aa for the species 'a' x a = v √ 2v T a (A.7)

ν aa = 4 √ π 3 ln Λ (4πε 0 ) 2 N a Z 4 a e 4 m 2 a v 3
x b = x ba x a = v √ 2v T b (A.8) x ba = v T a v T b (A.9)
The functions Ψ, G, Φ depend on the velocity modulus only and are defined as The function Φ(x) is the error function, and the function G(x) is the Chandrasekhar function, in accordance with Hirshman-Sigmar [45,46] and Hinton-Hazeltine [START_REF] Hinton | Theory of plasma transport in toroidal confinement systems[END_REF][START_REF] Hinton | Collisional Transport in Plasma[END_REF] papers. These definitions have been chosen such that the deflection and slowing-down frequencies ν d,ab and ν s,ab coincide with the definitions given by Hirshman and Sigmar [45]. The notation ν HS ab refers to the interspecies collision rate 1/τ ab defined in [40,47]. Hence ν HS ab is different from the inter-species momentum transfer rate ν ab . In particular ν HS aa differs from ν aa by a factor

Figure 1 . 6

 16 Figure 1. For a given radial profile of deuterium collisionality (blue dashed line), corresponding profiles of the collisionality of 3 trace impurities: helium (He), carbon (C) and tungsten (W, Z W = 40 here). The limits of the three collisionality regimes (banana, plateau, Pfirsch-Schlüter) are also displayed.

Figure 2 .

 2 Figure 2. Left: the 3 impurity density profiles considered for GYSELA simulations. Right: the time averaged neoclassical flux computed for each of these profiles is plotted versus the normalized density gradient in the case of 3 impurities (He, C and W) in a deuterium plasma. The slope gives the neoclassical diffusion coefficient.

Figure 3 .

 3 Figure 3. Diffusion coefficient of 3 impurities (He, C and W) in a deuterium plasma versus the impurity collisionality ν * z . The GYSELA results (o) are compared with the values predicted by Hirshman-Sigmar [45] (+) and the ones obtained with the NEO code (x).

Figure 4 .

 4 Figure 4. Detail of the method used to retrieve the pinch velocity of 3 impurities (He, C and W) in a deuterium plasma (no temperature gradient). The slope provides the diffusion coefficient, the intersection with the y-axis gives the pinch velocity. Intersection points with the x-axis correspond to the expected equilibrium gradient of the impurity density profiles as a result of the sole main ion density gradient (density peaking).

Figure 4

 4 Figure 4 shows the results of the 9 simulations. The extrapolated values on the y-axis provide the estimate of pinch velocities in GYSELA. These values are plotted on Fig. 5 as a function of collisionality. Theoretically, the normalized pinch velocity is expected to scale like RV pinch z

Figure 5 .

 5 Figure 5. Normalized pinch velocity of 3 impurities (He, C and W) in a deuterium plasma versus the impurity collisionality ν * z . The GYSELA results (o) are compared with the values predicted by Hirshman-Sigmar (+) and the ones obtained with the NEO code (x).

Figure 6 .

 6 Figure 6. Profile of the thermal screening factor H from a D-W simulation of GYSELA (dots) compared to the value predicted by Hirshman-Sigmar (line).

Figure 7 .

 7 Figure 7. Top row: poloidal fluctuations of the actual parallel velocities and flux obtained with GYSELA. Bottom row: same quantity computed from the radial gradients (cf. eqs.(36)(37)(38)). The agreement between theory and modeling looks satisfactory.

Figure 8 .

 8 Figure 8. Time and flux-surface average of the radial impurity flux in D-Helium, D-Neon and D-Tungsten simulations. The total flux (red) is decomposed into its electric (blue) and magnetic (black) drift components.

Figure 9 .

 9 Figure 9. Flux-surface average of the neon radial flux as a function of time and normalized radius. The total flux (c) is decomposed into its electric (b) and magnetic (a) drift components.

Figure 10 .

 10 Figure 10. Time-averaged radial tungsten flux in neoclassical (blue), turbulent (green) and full (red) simulations. The sum of the fluxes of the neoclassical and turbulent simulations is shown in black.

Figure 11 .

 11 Figure 11. Time-averaged poloidal cross-sections of the (m, n) = (1, 0) component of the electric potential for the 3 simulations in the case of a deuterium-tungsten plasma.

Figure 12 .

 12 Figure 12. Instantaneous poloidal cross-sections of tungsten density fluctuations δN W = N W -N W F S (with ... F S the flux-surface average) for the 2 reduced simulations in the case of a deuterium-tungsten plasma.

3 ) 6 )

 36 With these notations, the velocity, deflection and slowing-down collision frequencies are ν d,ab (x a ) = ν HS ab Here x b is a function of x a , i.e.

Ψ 2

 2 (Φ(x) -xΦ (x)) (A.12) Φ(x) = 2 √ π x 0 dy exp -y 2 (A.13)

√ 2 ,

 2 i.e. ν HS aa = √ 2ν aa . Useful asymptotic limits are slow particle (x → 0) : Ψ(x) → 1; Θ(x) → 1 (A.14) fast particle (x → ∞) : Ψ(x) → 3 √

Acknowledgments

The authors acknowledge the fruitful and lively discussions held during the 8th Festival de Théorie, held in Aix-en-Provence, July 2015. This work was granted access to the national HPC resources of CCRT, CINES, and IDRIS under GENCI (Grand Equipement National de Calcul Intensif) allocations and also with Helios/IFERC machine dedicated to Fusion at Rokkasho in Japan and with PRACE (Partnership for Advanced Computing in Europe) Research infrastructure. This project has received funding from the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 633053. The views and opinions expressed herein do not


Member, IEEE Joseph S Friedman
email: joseph.friedman@u-psud.fr

Member, IEEE Laurie E Calvet

Pierre Bessière

Jacques Droulez

Member, IEEE Damien Querlioz

I J S Belykh

L E Friedman

D Calvet

ON CIRCUITS AND SYSTEMS-I: REGULAR PAPERS 1 Bayesian Inference With Muller C-Elements

Keywords: Fault-tolerant circuit design, Muller C-element, nanotechnology, stochastic computing, variation-prone devices

Bayesian inference is a powerful approach for integrating independent conflicting information for decision-making. Though an important component of robotic, biological, and other sensory-motors systems, general-purpose computers perform Bayesian inference with limited efficiency. Here we show that Bayesian inference can be efficiently performed with stochastic signals, which are particularly adapted to novel low power nanodevices that exhibit faults and device variations. A simple Muller C-element directly implements Bayes' rule. Complex inferences are performed by C-element trees, which compute the probability of an event based on multiple independent sources of evidence. A naïve Bayesian spam filter circuit is demonstrated as a pedagogical application, and design techniques for improving circuit functionality are described. Limitations of this structure are discussed in terms of signal autocorrelation. The stochastic inference structure is exceptionally robust to faults, an essential feature of decision circuits, and can therefore leverage the increased efficiency of emerging nanodevices. This hardware implementation of Bayesian inference is extremely area and power efficient, with an area-energy-delay product several orders of magnitude less than the conventional floating point implementation. These results open a pathway for a direct stochastic hardware implementation of Bayesian inference, enabling a new class of embedded decision circuits for robotics and medical applications.

I. INTRODUCTION

B AYESIAN inference permits decision-making with max- imal information integration from distinct inputs [START_REF] Bessière | Probabilistic Reasoning and Decision Making in Sensory-Motor Systems[END_REF]. It computes the probability of an event through the incorporation of information from numerous sources. Bayesian inference has been suggested as a fundamental component of biological systems [START_REF] Ma | Bayesian inference with probabilistic population codes[END_REF]- [START_REF] Houillon | The probabilistic cell: Implementation of a probabilistic inference by the biochemical mechanisms of phototransduction[END_REF], and has been successfully applied to robotics and other sensory-motor systems [START_REF] Lebeltel | Bayesian robot programming[END_REF].

Conventionally, Bayesian inference is performed by generalpurpose processors in which the operations on probability values are executed by floating point arithmetic logic units.

These systems provide high precision, but consume significant energy and area. In addition to the substantial resources required by floating point operations, conventional approaches do not make efficient use of the high parallelism of Bayesian inference [START_REF] Bessière | Bayesian Programming[END_REF]. Furthermore, conventional systems require exact values throughout the computation, impeding the exploitation of variation and error-prone emerging nanodevices that can compute with reduced power consumption [START_REF] Bernstein | Device and architecture outlook for beyond CMOS switches[END_REF]. It is therefore highly desirable to devise a specialized hardware structure to directly compute Bayesian inference that is robust to the imperfections characteristic of novel nanotechnology.

Several distinct approaches have been proposed for a specialized Bayesian hardware, including fully digital systems with dedicated architectures [START_REF] Pournara | FPGAaccelerated Bayesian learning for reconstruction of gene regulatory networks[END_REF], [START_REF] Lin | High-throughput Bayesian computing machine with reconfigurable hardware[END_REF], and analog computation of probabilistic inference [START_REF] Vigoda | Analog logic: Continuous-time analog circuits for statistical signal processing[END_REF], [START_REF] Mroszczyk | The accuracy and scalability of continuoustime Bayesian inference in analogue CMOS circuits[END_REF]. Here we propose the stochastic computation of Bayesian inference, with probabilistic signals encoded as clocked bitstreams. Stochastic computing permits complex computations with minimal hardware, but its inherent imprecision and correlations have impeded its development.

We show that stochastic Bayesian inference can be performed directly with Muller C-elements, using surprisingly simple circuits. The increase in precision with run-time inherent to stochastic computation [START_REF] Neumann | Probabilistic logics and the synthesis of reliable organisms from unreliable components[END_REF], [START_REF] Gaines | Stochastic computing systems[END_REF] is particularly well-suited to decision circuits. The C-element circuit performs useful inferences despite the imprecision and correlations of stochastic computing and process variations of scaled nanodevices [START_REF] Gaudet | Iterative decoding using stochastic computation[END_REF]- [START_REF] Querlioz | Bioinspired programming of memory devices for implementing an inference engine[END_REF]. Furthermore, the fault tolerant nature of stochastic computing permits the exploitation of the extreme low energy, high speed, and small size of emerging nanotechnology that can serve as low-area, low-power computing [START_REF] Hamilton | Stochastic electronics: A neuro-inspired design paradigm for integrated circuits[END_REF] and random number generation units [START_REF] Fukushima | Spin dice: A scalable truly random number generator based on spintronics[END_REF], [START_REF] Choi | A magnetic tunnel junction based true random number generator with conditional perturb and real-time output probability tracking[END_REF].

In this work, we thus leverage stochastic computing to perform Bayesian inference with a simple circuit structure. We use the spam-detection problem as an instructive example to analyze circuit performance, as it is a common application of Bayesian inference in commercial products [START_REF] Bessière | Bayesian Programming[END_REF]. We explore limitations to the inference circuit structure, and discuss the mitigation of these concerns. Finally, we highlight that in comparison to the conventional computation of Bayesian inference, the proposed circuit structure can achieve sufficient precision with orders-of-magnitude increases in performance and energy efficiency.

II. STOCHASTIC BITSTREAMS

Stochastic computing uses traditional logic gates in a unique manner to perform operations on clocked binary stochastic single-bit streams rather than standard multi-bit binary signals.

. This allows for the compact implementation of many complex operations-for example, the product of probabilities is computed by a simple AND gate [START_REF] Gaines | Stochastic computing systems[END_REF]. Stochastic computing is therefore an attractive alternative to CMOS for inherently probabilistic computation, and its data encoding makes it particularly effective for exploiting error-prone nanodevices.

A. Encoded Probability

In a randomly generated stochastic bitstream, the dynamics of a 0-> 1 (1-> 0) switch can be defined by a state switching probability a (b). The probability R encoded by a stochastic bitstream is

R = a a + b (1)
which, for an uncorrelated bitstream (i.e., a + b = 1), is equivalent to

R = a. ( 2
)
Each bit in this bitstream thus has a probability R of being "1," and a probability 1 -R of being "0." Throughout this work, the percentage of "1"s in a bitstream encodes its probability.

B. Switching Rate

A particular probability can be represented by several distinct bitstreams that all encode the same value. For example, both bitstreams "0101010101" and "0000011111" encode the value 0.5, but the second bitstream contains longer "domains" of consecutive "0"s and "1"s. This first bitstream therefore allows for an evaluation of the value with significantly fewer bits of information. For this contrived example, the value encoded by the first N bits of the bitstream is closer to 0.5 for the first bitstream than the second bitstream, for any value of N . As will be seen in the following sections, this autocorrelation has profound implications for Bayesian inference.

It should be noted that the switching rate of a stochastic bitstream is significantly higher for signals with probabilities close to 0.5 than for extreme probabilities close to 0 or 1. For an uncorrelated stochastic bitstream with value R, the switching rate is given by

S = 2R(1 -R) = 2a(1 -a).
(3)

C. Precision

The probabilistic nature of stochastic bitstreams implies that there is no guarantee that the bitstream encodes the precisely correct value. In a properly operating circuit, however, the probability encoded by the bitstream converges toward the correct value as its length increases. The variance of an uncorrelated stochastic bitstream of length N is

Var = R(1 -R) N . ( 4
)
The expected variance of a bitstream thus decreases with 1/N , with smaller values for extreme probabilities. The strong effect of probability R on variance creates a bias that distorts the precision analysis, thereby making variance a poor metric of bitstream precision. A more appropriate metric for the error of a probabilistic signal is the Kullback-Leibler (KL) divergence

KL = (1 -R) ln 1 -R 1 -Q + R ln R Q (5)
which measures the information lost when a probability Q is used to approximate a probability R. The KL divergence is designed to measure error independently of R and, interestingly, the bitstream divergence approaches 1/2N for increasing N .

III. TWO-INPUT C-ELEMENT INFERENCE

To obtain the probability of a particular event V , Bayesian inference incorporates the probability of V given the prior P (V ) and evidence input E 1 as

P (V |E 1 ) = P (E 1 |V )P (V ) P (E 1 |V )P (V ) + P (E 1 |V )P (V ) (6)
which may be rewritten

P (V |E 1 )= P * (E 1 )P (V ) P * (E 1 )P (V )+(1-P * (E 1 )) (1-P (V )) (7)
with parameter P * (E 1 ) defined by

P * (E 1 ) ≡ P (E 1 |V ) P (E 1 |V ) + P (E 1 |V ) . ( 8
)
We now show that ( 7) is computed by a Muller C-element, enabling efficient inference circuits.

A. Muller C-Element

The Muller C-element, characterized by the truth table of Fig. 1(a), is a two-input memory element [START_REF] Muller | Theory of Asynchronous Circuits[END_REF]. The output Z maintains its state Z prev unless both inputs X and Y are opposite the current output state, in which case the output switches to the shared input value. Several circuits perform this function with eight to twelve transistors, as shown in Fig. 1(b)-(d), each with its own tradeoffs [START_REF] Shams | Modeling and comparing CMOS implementations of the C-element[END_REF]. Individual C-elements have recently found use in stochastic computing, particularly for error correction applications [START_REF] Gaudet | Iterative decoding using stochastic computation[END_REF], [START_REF] Winstead | A probabilistic LDPC-coded fault compensation technique for reliable nanoscale computing[END_REF], [START_REF] Winstead | C-element multiplexing for fault-tolerant logic circuits[END_REF].

B. Stochastic C-Element Inference

For C-element input signals i with no autocorrelation (a i + b i = 1 for switching probabilities a i and b i ), the output probability is

P (Z) = P (X)P (Y ) P (X)P (Y ) + (1 -P (X)) (1 -P (Y )) . ( 9
)
Substituting P * (E 1 ) for P (X), P (V ) for P (Y ), and P (V |E 1 ) for P (Z), ( 9) is equivalent to [START_REF] Bernstein | Device and architecture outlook for beyond CMOS switches[END_REF], demonstrating the Bayesian inference performed by C-elements. This equivalence motivates the implementation of Bayesian inference with stochastic computing, and shows that a simple C-element performs the complete inference of Bayes' rule.

C. C-Element Switching

Due to the inertia inherent to C-elements, the output switches equally often or or less frequently than the inputs. In contrast to a randomly generated bitstream, the switching rate of the C-element output bitstream is not uniquely determined by the probability it encodes, but is dependent on the combination of input probabilities. The upper bound of the output switching rate is the maximum switching rate of the two inputs. In particular, for uncorrelated bitstreams with opposite values (a i + b i = 1, R X + R Y = 1) and therefore equal switching rates, the expected switching rate of the output is one half of the input switching rate. Furthermore, the expected output switching rate when one input signal carries the value 0.5 is one half of the switching rate of the other input. The decreased switching causes a memory effect that results in autocorrelation-that is, the state of each bit is correlated to the states of other bits in the bitstream.

This autocorrelation causes a decrease in the amount of precision provided by a bitstream. The stochastic bitstream output of the C-element is also probabilistic and imprecise, but converges more slowly toward the precise Bayesian inference. In particular, the lower the switching rate of the output, the longer the "domains" of consecutive "0"s and "1"s, and the more imprecise the bitstream. Therefore, as can be seen in the simulation results shown in Fig. 2, the error is significantly increased for extreme inputs with low switching rates. The output signals polynomially converge to the exact Bayesian inference as N increases, but have consistently larger errors for opposite extreme input combinations. As these extreme opposite inputs can be considered as strong conflicting evidence, this error can be understood as inference uncertainty.

IV. MULTI-INPUT CASCADED C-ELEMENT INFERENCE

A. Cascaded C-Element Tree

The integration of additional inputs provides the capability for more complex inferences. Incorporating E 2 into the inference gives the probability of event V :

P (V |E 1, E 2 ) = P (V |E 1 )P (E 2 |V, E 1 ) P (V |E 1 )P (E 2 |V, E 1 ) + (1 -P (V |E 1 )) P (E 2 |V , E 1 ) . ( 10
)
Under the naïve assumption that evidence inputs E 1 and E 2 are conditionally independent, P (E 2 |V, E 1 ) = P (E 2 |V ), and (10) reduces to

P (V |E 1, E 2 ) = P (V |E 1 )P * (E 2 ) P (V |E 1 )P * (E 2 ) + (1 -P (V |E 1 )) (1 -P * (E 2 ))
. [START_REF] Mroszczyk | The accuracy and scalability of continuoustime Bayesian inference in analogue CMOS circuits[END_REF] This can in theory be computed using two cascaded C-elements, as illustrated in Fig. 3(a).

TABLE I SAMPLE MESSAGES FOR SPAM-DETECTION

This inference can be extended to multiple independent evidence inputs E 1 , E 2 , . . . , E i to give [START_REF] Neumann | Probabilistic logics and the synthesis of reliable organisms from unreliable components[END_REF], shown at the bottom of the page. This multi-input inference can be approximated with cascaded C-elements [Fig. 3(b)]. The final output is expected to be

P (V |{E 1 , . . . , E i }) = P (V ) i x=1 P * (E x ) P (V ) i x=1 P * (E x ) + (1 -P (V )) i x=1 (1 -P * (E x )) . ( 13
)
Straightforward arithmetic shows that the C-elements can equivalently be organized in trees as in Fig. 3(c).

B. Example: Bayesian Spam Filter Circuit

To illustrate the feasibility of multi-input Bayesian inference with cascaded C-elements, we evaluate its use for spam filtering [START_REF] Bessière | Bayesian Programming[END_REF], a canonical Bayesian inference. Final applications of our proposal will target situations in which decisions are automated, such as brain-machine interfaces and robotics, and for which an efficient direct hardware implementation would fuel tremendous advances.

We consider the messages in Table I: A and B are likely to be ham (i.e., not spam), C is almost definitely spam, D seems to be ham but contains many words that are associated with spam, E is probably spam, and F is entirely inconclusive. To infer the probability that a message is spam, a simple Bayesian spam filter considers the "spamicity" P * (W m x x ) of each word W x given its presence (m x = 1) or absence (m x = 0). Statistics are provided in Appendix A. The probability that a message is spam given the presence or absence of various keywords is computed in accordance with Section IV-A. Assuming no prior information (P (S) = 0.5) the posterior probability reads

P (S|{W m 1 1 , . . . , W m i i }) = i x=1 P * (W m x x ) i x=1 P * (W m x x ) + i x=1 (1 -P * (W m x x )) . ( 14
)
The three-stage C-element tree of Fig. 3 can be used to perform the spam filter operation on messages A-F for eight key words. The presence or absence of each word provides a Simulation results for several runs of the spam filter are shown as points in the plot of Fig. 4, with the exact naïve Bayesian inference shown as a solid line. This simple Bayesian system successfully filters the messages, with high outputs for message C and E, low outputs for messages A, B, and D, and an ambiguous output for F. These decisions can be reached reliably by observing bitstreams of 1000 bits. Despite the inertia and autocorrelation, these simulation results show that the inference performed by individual C-elements can be extended to multiinput Bayesian inference with cascaded C-elements.

As can be seen in the figure, there is significant variation of the output values for short bitstreams. As the bitstream length increases, the outputs converge toward values that are quite close to the exact Bayesian inference. However, the output is not perfectly exact, due to the inertia resulting from the autocorrelated domains. The largest discrepancy is for message E: the output converges to a value of 0.7, while an exact inferences gives a value of 0.78. It should further be noted that message D converges slower than other messages, which is a result of the extreme values at its input. We now discuss the origin of this non-ideal behavior.

C. Cascaded C-Element Switching & Inertia

Beyond the first stage of the tree, the inputs to C-elements are autocorrelated; that is, a i + b i = 1. More specifically, following the reasoning of Section III-C, a i + b i ≤ 2 -n , where n is the number of C-elements through which information has passed since the original stimulus. The C-elements can therefore no longer be considered as isolated black boxes, and the autocorrelation between consecutive bits B t and B t+1 of an output bitstream B can be described in terms of the mutual information I as

P (V |{E 1 , . . . , E i }) = P (V |{E 1 , . . . , E i-1 }) P * (E i ) P (V |{E 1 , . . . , E i-1 }) P * (E i ) + (1 -P (V |{E 1 , . . . , E i })) (1 -P * (E i )) (12)
I(B) = 1 + (1 -S) log 2 (1 -S) + S log 2 (S) (15)
for a bitstream encoding a value 0.5 with switching rate S defined as

S = P (B t = 1 ∧ B t+1 = 0) + P (B t = 0 ∧ B t+1 = 1). (16)
As signals propagate through successive stages, the inertial nature of the C-elements induces strings of consecutive "1"s and "0"s. As mentioned previously, these "domains" suggest an increase in autocorrelation, and are therefore detrimental to the operation of the inference circuit. This domain growth is shown in Fig. 5, which presents the distribution of domain sizes at each stage for a tree with input values of 0.5. Though the input bitstream domains are quite short, the domains grow as the signals propagate through the circuit. By the fifth stage, which is the result of 32 input signals, the domain distribution deviates from that expected for uncorrelated stochastic bitstreams, with the most common domain length being 17 bits.

As the output of a C-element only switches in the case of both inputs having values opposite that of the output, the C-elements have a tendency to propagate signals with values biased toward the input values. More specifically, "inertia" causes the C-elements to be biased toward strongly correlated inputs with large domains. These correlation issues have been investigated for stochastic decoding applications, and several techniques have been described to decorrelate stochastic bitstreams [START_REF] Sharifi Tehrani | Relaxation dynamics in stochastic iterative decoders[END_REF]. The proposed counter and moving average circuits can reshuffle the bitstream and decrease the correlation, but come at the cost of increased area, power dissipation, and delay. This overhead therefore reduces the efficiency of the cascaded C-element inference circuit. Inherently probabilistic nanodevices [START_REF] Knag | A native stochastic computing architecture enabled by memristors[END_REF], [START_REF] Fukushima | Spin dice: A scalable truly random number generator based on spintronics[END_REF], [START_REF] Choi | A magnetic tunnel junction based true random number generator with conditional perturb and real-time output probability tracking[END_REF], [START_REF] Saighi | Plasticity in memristive devices for spiking neural networks[END_REF], [START_REF] Vincent | Spintransfer torque magnetic memory as a stochastic memristive synapse for neuromorphic systems[END_REF] may be envisioned for the random number generators necessary in these decorrelation circuits, potentially bringing a substantial reduction in overhead. Decorrelation circuits can be periodically inserted between cascaded C-elements, with optimal design dependent on the tradeoff between efficiency and accuracy. Decorrelation circuits should therefore be considered for applications that require a level of accuracy unavailable with the simple cascaded C-element circuit.

D. Cascaded C-Element Precision

The importance of domain size is evident from Fig. 6, which shows the error for various bitstream lengths for several combinations of inputs. Similar to the case of a single C-element, the use of extreme probabilities leads to large domains and slow convergence. With cascaded gates, the intermediate signals can take on extreme values, resulting in massive errors for the green line with inputs of 0.01, 0.01, 0.99, and 0.99. The intermediate values are 10 -4 and 1-10 -4 , preventing the representation of meaningful information with fewer than 10 000 bits. There are thus large output errors despite the expected output value being 0.5.

The effect of cascading on precision is studied in more depth in Fig. 7. Here, there is a tree structure with n stages and 2 n input stochastic bitstreams. The input values for one half of the bitstreams are 0.5k, where k is a constant plotted on the x-axis. The second half of the bitstreams have input values of 0.5 + k. This is the least favorable way to organize inputs, as the inputs of the final C-element encode extreme probabilities (0.5k) n-1 and (0.5 + k) n-1 . Fig. 7(b) shows that the error and autocorrelation increases with increasing n and k, as the signals propagated by the second-to-last stage grow larger domains and become more extreme. The simulations thus show the importance of intermediate signals in determining output precision, and explores the autocorrelation limitations on the number of inputs and their values. Therefore, when possible, the inputs should be organized to minimize the potential for inertia and extreme values at intermediate nodes in the circuit.

Because of these considerations, the behavior of cascaded C-elements deviates from that of a single C-element, and the output of the C-element tree can diverge from the Bayesian inference of [START_REF] Gaines | Stochastic computing systems[END_REF]. Despite the challenges posed by autocorrelation, we have shown with the pedagogical example of spam filtering that cascaded C-elements can give useful results in practical situations.

V. COMPUTING IMPLICATIONS

The use of stochastic computation brings exceptional fault tolerance, and the application-specific hardware enables high efficiency. This direct hardware implementation of Bayesian inference provides significant computing advantages which are particularly pronounced with error-prone nanodevices.

A. Fault Tolerance

For use with emerging nanotechnology, circuits should be designed to tolerate faults. Stochastic computation was originally designed for imprecise components, and fault tolerance is one of its primary strengths. Each incorrect bit in a stochastic bitstream causes a deviation from the value by only 1/N , where N is the length of the bitstream. Additionally, if it is equally likely for a "0" to improperly represented by a "1" as for a "1" to be represented by a "0," then these faults may cancel out.

1) Not Switching: For the various C-element circuits shown in Fig. 1, the most likely fault is a missed switching event due to an insufficient current being provided to overpower the latch [START_REF] Brzozowski | Testing C-elements is not elementary[END_REF]. This type of fault will be referred to as a "not switching" fault, and can result from the use of a low supply voltage for circuits composed of variability-prone nanodevices.

The effect of this "not switching" fault can be seen in Fig. 8, which shows the impact of fault rates up to 99%. The simulation reflects the arbitrary case of the two inputs X and Y having probability values of 0.7 and 0.1; the expected output value from ( 9) is 0.206. One million bits are used in this simulation. As seen in the graph, even for extremely high error rates, the output is close to the expected value. The KL divergence remains below 0.01, which may be considered a "correct" value, as discussed in the next subsection.

The error resulting from the "not switching" fault is unbiased: the output deviates from the correct value, but is not consistently larger or smaller. This is due to the nature of the "not switching" fault: as this fault also causes incorrect values in proceeding bits and there are necessarily equal numbers of switching events from "0" to "1" as from "1" to "0," there is an increase in domain length without biasing the output. This "not switching" fault is therefore equivalent to decreasing the bitstream length by the fault rate. Replacing the x-axis of Fig. 8 by 1 000 000 × (1 -F notswitching ), these results can be compared to the C-element outputs of Fig. 2. In both cases, the KL divergence rises from slightly greater than 10 -6 to slightly greater than 10 -4 as the bitstream length decreases from 1 000 000 to 10 000.

The exceptional tolerance of the "not switching" fault permits the consideration of advanced low-power design techniques with emerging nanotechnology. Reduction of the supply voltage provides a drastic decrease in power dissipation, but increases the likelihood of faults. Additionally, the clock speed can be increased to improve performance in exchange for an increase in faults. This fault tolerance thus provides an opportunity to optimize the power dissipation and speed for various applications and required levels of precision.

2) False Switching: If the latch itself is faulty, there is a possibility for the C-element to switch in the absence of a proper switching stimulus. This "false switching" fault is less likely than the non-switching fault. For realistic fault rates, far less than 1% in conventional circuits, the C-element structure provides the approximately correct output (KL divergence less than 0.001 for 1% fault rate). When the fault rate is increased toward 100%, the output approaches a random bitstream with probability 0.5. This can be understood as the expected result of the output switching every clock cycle-both those in which it is expected to switch, and those in which no switching event is expected. While this is not a realistic concern, this bias toward 0.5 can be interpreted in the context of a decision-making process as uncertainty arising from unreliable nanoelectronic components.

B. Inference Efficiency

This C-element Bayesian inference structure uses a simple circuit to stochastically compute inferences. The precision increases as more bits are used for the inference, at the cost of increased energy and delay. In comparison to conventional implementations with floating point operations, the C-element inference structure provides significant improvements in area, performance, and energy consumption. However, when performing benchmark tests, the "correctness" of a Bayesian inference is application-specific. Here, a KL divergence of less than 0.01 is considered "correct" for an unspecified general circuit. Using [START_REF] Lebeltel | Bayesian robot programming[END_REF], this divergence corresponds to a value of 0.5 being approximated as 0.57, a value of 0.1 being approximated as 0.15, and a value of 0.01 being approximated as 0.03. Clearly, more precision is necessary for some applications; for others, less precision is necessary.

To perform a comparison, the circuit efficiency is evaluated under identical technology node and process. Standard cells from Synopsys' SAED_EDK90_CORE library are used [START_REF] Synopsys | Digital Standard Cell Library: SAED_EDK90[END_REF]; the efficiency of each component is considered at its ideal operating conditions. For the C-element, the circuit of Fig. 1(d) is chosen despite its non-ideality, as it can most easily be translated to standard cells. For the conventional floating point units, the parameters are optimized for the SAED_EDK90_CORE library by Tziantzioulis et al. [START_REF] Tziantzioulis | b-HiVE: A bit-level history-based error model with value correlation for voltage-scaled integer and floating point units[END_REF]. These parameters are described in Appendix B.

1) Area: For a two-input Bayesian inference, this stochastic circuit requires just one C-element. This is in contrast to the conventional floating point implementation, which requires addition, multiplication, and division units. With the standard cells described above, the floating point circuit consumes ∼16 000 times more area than a C-element. This estimate ignores the overhead associated with random number generation, which may be efficiently achieved with emerging nanodevices [START_REF] Gaba | Stochastic memristive devices for computing and neuromorphic applications[END_REF], [START_REF] Fukushima | Spin dice: A scalable truly random number generator based on spintronics[END_REF], [START_REF] Choi | A magnetic tunnel junction based true random number generator with conditional perturb and real-time output probability tracking[END_REF] and will therefore become essential elements of stochastic computers.

For multi-input Bayesian inference, a single floating point circuit can be used. The stochastic C-element inference structure, however, requires T -1 C-elements, where T is the number of inputs. Therefore, the C-element provides an area reduction of ∼5400x for a two-stage (four-input) inference and ∼2300x for a three-stage (eight-input) inference. Similar calculations can be performed for inferences with additional inputs.

2) Performance: Whereas the stochastic C-element inference structure provides increasing precision with increasing time, the floating point implementation provides a highly-precise inference after its characteristic latency period. Therefore, the performance benefit provided by the C-element structure is dependent on the necessary precision for the specific application.

For multi-input Bayesian inference, the multi-stage C-element circuit needs only one additional cycle per stage to compute a bitstream with a particular length. However, the autocorrelation and domain growth contribute to a loss of precision, requiring longer bitstreams and slightly more computation time to achieve a precise output. The floating point circuit, by contrast, performs multiple pipelined computations and computes a highly precise output after its characteristic delay time. These processes are compared in Fig. 9, which considers a tree with n stages and 2 n input bitstreams all encoding the probability 0.5. As can be seen in the figure, the times required to achieve a KL divergence of 0.01 are of the same order of magnitude.

3) Energy Consumption: The C-element inference circuits consume a static energy proportional to the bitstream length in addition to a dynamic energy related to the switching rate, with the dynamic energy consumption largely dominant under typical conditions. Therefore, the energy required for a C-element to output a given value with a particular level of precision is nearly independent of input probability values. For extreme inputs, a large number of bits is required, but the effect is counteracted by the low switching rate. In comparison to the floating point implementation, there is an enormous reduction in the energy required to produce a KL divergence of 0.01, of roughly three orders of magnitude.

The multi-input C-element inference provides similarly large reductions in energy consumption compared to the floating point implementation. As shown in Fig. 10, the energy consumption of both approaches increases with an increasing number of stages. These results show a ∼1000x reduction in energy consumption to produce an output with a KL divergence of 0.01.

VI. CONCLUSION

This work highlights the ability to perform Bayesian inference with simple circuits that can exploit emerging nanotechnology. In a stochastic approach, a circuit composed of Muller C-elements can integrate conflicting evidence from independent sources to provide an increasingly precise approximation of the probability of an event. Although extensions that provide a more formal analysis of the autocorrelation and explore inertia mitigation through signal rerandomization are envisioned, the C-element inference structure introduced here already achieves sufficient accuracy for some practical applications, and validates the basic concept of using stochastic computing to perform Bayesian inference.

The efficiency and fault tolerance of the C-element inference circuit provides significant advantages over conventional software implementations using floating point hardware. As shown in Fig. 11, the area-energy-delay product (AEDP) for producing a KL divergence of 0.01 is roughly one billionth of the floating point implementation. In addition, the exceptional fault tolerance provide an impetus for the use of stochastic C-elements in nanoelectronic decision circuits based on Bayesian inference.

For decision circuits, it is useful to have a constantly increasing precision with an initial approximation. Though floating point operations remain more effective for applications involving extreme values or requiring exact outputs, stochastic computing may provide important benefits for applications such as robotics. It also inspires a pathway toward compact embedded inference circuits with long battery life suitable for medical applications. The stochastic circuit structure proposed here has the potential to move Bayesian inference from software to nano-hardware, thereby enabling enormous improvements in decision-making efficiency.

APPENDIX A BAYESIAN SPAM FILTER

A naïve Bayesian spam filter searches a message for particular key words and computes the probability of the supposition the message is spam based on their presence or absence. As shown in Table II, eight key words are chosen and used for evaluating the messages. For each word, Laplace's rule of succession is applied to correctly handle incomplete observation:

P W 1 x |S = 1 + #SpamWithW x 2 + #Spams (17)
P W 1 x |H = 1 + #HamWithW x 2 + #Hams . ( 18
)
The spamicity is calculated in accordance with Bayes' rule:

P S|W 1 x = P (S)P W 1

x |S P (S)P (W 1

x |S) + P (H)P (W 1 x |H)

. [START_REF] Hamilton | Stochastic electronics: A neuro-inspired design paradigm for integrated circuits[END_REF] Assuming that P (S) = P (H) gives

P * W 1 x = P W 1 x |S P (W 1 x |S) + P (W 1 x |H) . ( 20
)
The probability that a message is spam given the absence of a particular key word is calculated similarly:

P * W 0 x = P W 0 x |S P (W 0 x |S) + P (W 0 x |H) . ( 21
)
The probability that a message is spam given the presence or absence of the various key words is calculated according to [START_REF] Gaudet | Iterative decoding using stochastic computation[END_REF].

For the case of message D, the equation is

P (S|Message D) = P S Nigeria 1 , Pizza 1 , Stochastic 0 , Check 1 , You 1 , Transfer 1 , Commute 1 , Fortune 0 . ( 22
)

APPENDIX B COMPARISON TO FLOATING POINT INFERENCE

A. Area Comparison

The areas of the various circuits are listed in Table III.

The circuit design has one hardware implementation for each function, and the results are pipelined when a particular function is used on multiple occasions. Whereas the area of the conventional inference circuit is independent of the number of inputs to the inference, the C-element approach requires 2 n-1 C-elements for an inference circuit with n stages.

B. Performance Comparison

The delays of the various circuits are listed in Table IV. For the conventional approach, ideal pipelining is assumed in which the operations are pipelined with a spacing equivalent to the number of issue cycles; there is therefore a decrease in delay per inference when more than one stage is combined for a cascaded circuit. These delay figures apply under the assumption that the separate hardware units of the floating point inference circuit each operate at their own ideal clock frequency. These assumptions are likely to underestimate the delay of the conventional inference.

C. Ernergy Consumption Comparison

The energy consumptions of the various circuits are listed in Table V. For the C-element, the expected energy consumed per bit is given for the worst case, with input values of 0.5. It is assumed that switching energy (but not power) is independent of frequency. For the conventional approach, it is assumed that the circuit dissipates static power throughout the complete inference calculation. In lieu of an accurate evaluation of the activity rate of the various circuit components, it is assumed that each cell of each arithmetic circuit switches once per arithmetic operation. These assumptions are likely to underestimate the energy consumption of the conventional inference.

Fig. 1 .

 1 Fig. 1. (a) Muller C-element truth table. (b) Compact C-element with weak latch. (c) Efficient one-stage C-element. (d) Standard cell C-element design.

Fig. 2 .

 2 Fig. 2. Mean simulated precision (KL divergence from exact value) of C-element inference for various input signal combinations as a function of observed bitstream length. Error bars depict 95% confidence interval.

Fig. 3 .

 3 Fig. 3. Cascaded C-element Bayesian inference circuits that integrate a prior and (a) two inputs; (b) four evidence inputs in a sequential structure; (c) seven inputs in a tree structure.

Fig. 4 .

 4 Fig. 4. Symbols: Spam probability value for messages A-E of Table I obtained by cascaded C-elements for several simulation runs as a function of observed bitstream length. Horizontal lines: exact inference. The output of the C-elments converges toward an approximated Bayesian inference as bitstream length increases.

Fig. 5 .

 5 Fig. 5. Normalized distribution of bitstream domain sizes for trees of various depths and input values of 0.5. The domain sizes of the output bitstreams increase with the number of C-element stages.

Fig. 6 .

 6 Fig. 6. Mean simulated precision of a two-stage C-element inference circuit, for various input signal combinations as a function of observed bitstream length. Error bars depict 95% confidence interval.

Fig. 7 .

 7 Fig. 7. (a) Cascaded C-element tree circuit with n = 3 stages to evaluate worst-case input probability combination. (b) Precision (solid lines) and autocorrelation (dashed lines) of tree as a function of input values for various depths with 500 000-bit bitstreams. Error bars depict 95% confidence interval.

Fig. 8 .

 8 Fig. 8. The simulated (black dots) and expected (black line) output value for "not switching" faults. The mean KL divergence (red) is below 10 -3 for a 99% fault rate. Error bars depict 95% confidence interval.

Fig. 9 .

 9 Fig. 9. Precision (KL divergence from exact inference) as a function of computing time for various inference circuits.

Fig. 10 .

 10 Fig. 10. Precision as a function of energy consumed by various inference circuits.

Fig. 11 .

 11 Fig. 11. AEDP as a function of precision for various inference circuits.

TABLE II WORD

 II COUNT STATISTICS FOR SPAMICITY CALCULATIONS

TABLE III

 III

	INFERENCE CIRCUIT AREA COMPARISON

ACKNOWLEDGMENT

The authors would like to thank N. Hardavellas, A.M. Gok, G. Tziantzioulis, D. Colliaux, E. Mazer, and J. Grollier for fruitful discussions.

This work was supported by the FP7 ICT BAMBI (FP7-ICT-2013-C) project and a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (Labex NanoSaclay, reference: ANR-10-LABX-0035). This paper was recommended