Y Camenen
email: yann.camenen@univ-amu.fr

F J Casson

P Manas

A G Peeters

Interplay between toroidal rotation and flow shear in turbulence stabilisation

Keywords:

The interplay between toroidal rotation u, parallel flow shear u ′ and perpendicular flow shear γE in the stabilisation of tokamak turbulence is investigated in non-linear flux-tube gyrokinetic simulations. The simulations are performed for a reference L-mode DIII-D plasma (the so-called shortfall case) at r/a = 0.8, varying the flow parameters around their nominal values. Depending on the respective signs of u, u ′ and γE, turbulence is found to be enhanced, reduced or unchanged. When the coupling is favorable, the overall effect on the non-linear heat fluxes can be very large, even at moderate flow values. The ion heat flux is for instance decreased by a factor of three when the direction of the parallel flow shear is reversed with respect to its nominal value. Even more surprising, keeping u ′ and γE at their nominal values, the ion heat flux decreases by more than 50% when the toroidal flow is reversed. The relevance of this mechanism in the experiments which depends on the ability to decouple u, u ′ and γE is discussed. The interplay between u and u ′ observed in the non-linear simulations qualitatively follows the linear stability results and is interpreted in the frame of a simple fluid model.

I. INTRODUCTION

Sheared flows have long been known to play an important role in the regulation of turbulence. Basically, two mechanisms are at play: on one side, the plasma advection by sheared flows tends to decorrelate the turbulent structures and decrease the turbulence saturation level, whereas on the other side sheared flows have a direct impact on linear modes that makes them more or less unstable depending on the cases [1]. For core tokamak plasma turbulence, the strong scale separation in the directions parallel and perpendicular to the magnetic field, typically k ⊥ /k ∼ qR/ρ i ≫ 1 (with q the safety factor, R the major radius of the torus and ρ i the ion Larmor radius), makes the shear in the perpendicular flow much more effective to decorrelate turbulent structures than the shear in the parallel flow. As a result, a radial shear in the perpendicular flow is usually stabilising, see e.g. [2], whereas the main effect of a radial shear in the parallel flow is to drive the parallel velocity gradient (PVG) instability [3,4] (also known as the Kelvin-Helmholtz instability). In the limit of pure toroidal rotation, a reasonable assumption in the plasma core in the absence of transport barriers, the perpendicular and parallel flows are geometrically linked by the ratio between the poloidal and toroidal magnetic field B p /B t (scaling as ε/q for circular flux surfaces at low inverse aspect ratio ε). When the toroidal flow shear increases, its perpendicular and parallel components increase in proportion and their impact on the turbulence is weighted according to the value of B p /B t : at large B p /B t , the perpendicular to parallel flow shear ratio is larger than at low B p /B t and stabili-sation is favored, the opposite being true at small B p /B t [5][6][7][8]. Recently, the impact of the toroidal flow itself on turbulence regulation has also been investigated. At low Mach number, the toroidal flow barely affects the turbulence saturation. However, for toroidal Mach numbers nearing or exceeding unity, ion temperature gradient and trapped electron mode turbulence is found to increase significantly with toroidal rotation, due to an increase of the linear trapped electron mode drive and to a reduced zonal flow residual induced by the centrifugal force [9]. The present work aims at investigating the interplay between toroidal rotation, parallel flow shear and perpendicular flow shear in turbulence stabilisation. The study is motivated by the unexpectedly strong impact of the toroidal flow on the non-linear heat fluxes observed at low Mach number with the gyrokinetic Eulerian code GKW [10] in the frame of a multi-code benchmark effort based on the DIII-D L-mode shortfall case [11]. In the following, the treatment of plasma flows in GKW is first presented, together with a brief reminder of the scaling of neoclassical flows in tokamaks. A simple fluid dispersion relation is then derived to highlight the coupling between toroidal rotation and parallel flow shear in the linear stability of the toroidal Ion Temperature Gradient mode (ITG). Guided by the fluid model results, the impact of the cross-term on the linear mode growth rate and its parametric dependencies are investigated in linear simulations. The impact of the flow parameters on the non-linear heat flux is then characterised and quantified in simulations based on the shortfall case. Finally, the relevance of this mechanism in the experiments is discussed and the results are summarized.

To proceed with the description of the flows implementation in GKW, we first define a right-handed toroidal flux-coordinate system (r, θ, ϕ). The minor radius r = 1 2 (R max -R min ) is a flux surface label, with R min and R max the minimum and maximum major radii of the flux surface, respectively. The poloidal angle θ is zero at the low field side midplane and increases towards the top of the torus. The toroidal angle ϕ is clockwise when viewed from above. As is customary, the axisymmetric equilibrium magnetic field is expressed as:

B = s b F ∇ϕ + s j ∇ϕ × ∇ψ (1)
where F = RB t > 0 and ψ is the poloidal magnetic flux. Note that the direction of the toroidal magnetic field and plasma current, s b and s j respectively (positive when clockwise from above), have been singled out in the expression above. As a consequence, all scalar quantities related to the magnetic field (the safety factor q for instance) are considered positive and the signs explicitly appear via s b and s j . Consistently with Eq. ( 1), the poloidal magnetic flux ψ is minimum on the magnetic axis and ∇ψ points in the same direction as ∇r.

B. Neoclassical flows

In tokamaks, under an expansion up to first order in ρ * = ρ i /R 0 , with R 0 the reference major radius, the flow of species s can be expressed as [12]:

V s = ûθ,s B + R 2 (ω p,s + ω Φ )∇ϕ (2)
where

ûθ,s = V s • ∇θ B • ∇θ (3)
ω p,s = - s j Z s en s ∂p s ∂ψ (4)
ω Φ = -s j ∂Φ ∂ψ (5)
with Z s , n s and p s the species charge number, density and pressure, respectively, and Φ the electrostatic potential. The poloidal component of the flow in Eq. ( 3) is strongly constrained by neoclassical physics and damped towards a species dependent offset ûneo θ,s on timescales of the order of the ion-ion collision frequency. The neoclassical poloidal flow is first order in ρ * and for main ions is given by

ûneo θ,i = -s b s j 1 2 ρ * v thi ∂r ∂ψ F B 0 < B 2 > k i R/L Ti (6)
where v thi = 2T i /m i is the ion thermal velocity, with T i the ion temperature and m i the ion mass, ρ * = m i v thi /(eB 0 R 0 ) is the normalised ion Larmor radius, with B 0 the reference magnetic field, R/L Ti = -R 0 ∂ log T i /∂r is the normalised ion temperature gradient and < . > denotes the flux surface average. The flow coefficient k i is dimensionless, typically of order unity, and depends on the plasma shape and collisionality (see [13] for typical dependencies on plasma parameters and comparisons between standard analytical limits and drift-kinetic simulations). In the banana regime (low collisionality), k i is positive which corresponds to a poloidal velocity in the ion diamagnetic drift direction. The diamagnetic contribution in Eq (4) is also first order in ρ * :

ω p,i = s j 1 2 ρ * v thi ∂r ∂ψ B 0 R/L pi (7)
In the limit of vanishing ρ * , the neoclassical flow is therefore purely toroidal [14] and given by:

V (0) i = ω Φ R 2 ∇ϕ (8)
with ω Φ constant on a flux surface. The flow physics implemented in δf flux-tube codes usually relies on a lowest order plasma flow described by Eq. ( 8), which is supported by the neoclassical theory.

C. Implementation in GKW

GKW is an Eulerian δf gyrokinetic code. The version used in the present study [10] is based on the gyrokinetic equations in the local approximation, assumes a Maxwellian background distribution function and a lowest order plasma flow given by Eq. (8). It describes the turbulence in a flux-tube wrapped around the flux surface labelled by r = r 0 . In this description, the flow physics enters in three distinct ways: via the flow itself (inertial effects), via the shear in the parallel flow and via the shear in the E × B flow. The implementation of these three components is briefly outlined below. A comprehensive description of the equations can be found in the GKW manual [15] and their derivation in the electrostatic limit is given in [16,17].

Inertial effects

GKW equations are formulated in a frame that rotates toroidally as a rigid body [16,17]. The angular frequency of the frame Ω is chosen to cancel the lowest order E × B drift at r = r 0 and is defined positive for a plasma flow in the direction of the toroidal magnetic field:

Ω = s b ω Φ (9)
The frame frequency is specified in the code via the normalised toroidal velocity u (ion Mach number) defined as:

u = R 0 v thi Ω (10)
In the rotating frame, the inertial effects associated to the background flow appear via the Coriolis and centrifugal forces. These two forces result in several modifications of the gyrokinetic equations presented in Appendix A. By definition of the frame frequency, all the inertial effects associated to the lowest order neoclassical flow given in Eq. ( 8) are taken into account. Inertial effects associated to departures from Eq. ( 8), i.e. finite poloidal rotation and/or species dependent toroidal rotation, are neglected in the present implementation, consistently with the local limit assumption, ρ * → 0. In most of the simulations presented in this paper, the Coriolis drift is the only inertial effect taken into account, i.e. centrifugal effects are neglected, which is appropriate at low Mach number.

Flow shear

The parallel flow shear is assumed to arise entirely from the shear in a toroidal flow of constant angular frequency on a flux surface. It enters the Maxwellian background as detailed in Appendix B 1. It is specified for each species by the normalised toroidal rotation u ′ s defined as:

u ′ s = - R 2 0 v thi ∂ω f s ∂r (11)
where the species angular frequency in the rotating frame ω f s is assumed to be constant on a flux surface. In the present study, all species are considered to have the same parallel flow shear and this common flow shear is noted u ′ . The advection by a sheared E × B flow is taken into account by linearising the E × B flow around r = r 0 (see Appendix B 2) and specified via the normalised E × B shearing rate:

γ E = R 0 v thi 1 B 0 ∂ 2 Φ f ∂r 2 (12)
where Φ f is the lowest order electrostatic potential in the rotating frame.

Departure from purely toroidal flows

Strictly speaking, the flow treatment presented above is only valid in the limit of a purely toroidal rotation described by Eq. ( 8). The parallel and E × B flow shears are then coupled and given by:

u ′ s = -s b R 2 0 v thi ∂ω Φ ∂r (13
)

γ E = -s j R 0 B 0 v thi ∂ψ ∂r ∂ω Φ ∂r (14)
In reality, tokamak flows are not purely toroidal and differ for each species, i.e. ûθ,s and ω p,s are finite in Eq. ( 2). Considering that the poloidal flow is neoclassical and described by Eq. ( 6), the departure from a purely toroidal flow scales as

ûθ,i B ∼ 1 2 ρ * v thi B t B p R/L Ti , (15)
where ∂ψ/∂r = RB p |∇r| ∼ RB p and k i ∼ 1 have been assumed. Provided the collisionality is slowly varying with the minor radius, the corresponding contribution to the flow shear then scales as

R 0 B ∂ ûθ,i ∂r ∼ 1 2 ρ * v thi B t B p [R/L Ti ] 2 (16)
where

R 2 0 T i ∂ 2 T i ∂r 2 = [R/L Ti ] 2 + R 0 ∂ ∂r [R/L Ti ] ∼ [R/L Ti ] 2 (17)
has been assumed and where the derivatives related to the magnetic equilibrium have been considered of order unity, e.g. R 0 ∂ ln B/∂r ∼ 1.

Formally, the contribution to the flow shear arising from Eq. ( 16) is first order in ρ * and should be neglected in the local limit. However, due to the B t /B p factor and the presence of a second derivative this contribution is not particularly small for typical plasma parameters. For ρ * = 1/600, B t /B p ∼ q/ǫ = 10 and R/L Ti = 6, Eq. ( 16) yield B∂ ûθ,i /∂r ∼ 0.3 v thi /R 0 , which is to be compared to the value of the parallel flow shear, 1 -3 v thi /R 0 , typically measured in the tokamak core. The assumption in Eq. ( 17) is rather conservative and in practice the second derivative of the temperature can easily be several times larger, especially in the case of internal transport barriers or closer to the plasma edge where R/L Ti increases rapidly with the minor radius. For the modelling of experimental plasmas with a significant departure from purely toroidal sheared flows, the values of u ′ s and γ E can be specified as described in Appendix B 3. It is with these cases in mind that the values of u ′ s and γ E are scanned independently in the present study. The reader is cautioned, however, that this independent variation does not strictly comply with the local limit ordering ρ * → 0 and selects one first order contribution in ρ * believed to be of particular importance. Before moving to the simulation results, the interplay between toroidal rotation u and parallel flow shear u ′ in the linear stability of the toroidal Ion Temperature Gradient mode (ITG) is explored with a simple fluid model.

The purpose of the model is to help interpreting the gyrokinetic simulations by emphasizing the mechanisms at play in a simplified framework. A detailed derivation of the model can be found in [16]. For completeness, its main characteristics are reminded below. The model is built from the first three moments of the linearised gyrokinetic equation in the local δf approximation and evaluated at the low-field-side midplane (strong ballooning assumption). It assumes an adiabatic electron response, neglects finite Larmor radius effect and uses a Maxwellian closure. The magnetic field is considered to be purely toroidal with s b = -1. In the following, the inertial effects associated to the Coriolis force are retained but centrifugal effects are neglected (they are retained in [16]). Singly charged ions Z = 1 with an equal ion and electron background temperature T i = T e are considered. The complex mode frequency ω = ω R + iγ, the parallel wave vector k and the perturbed electrostatic potential δφ are normalised as follows:

ω = ω ω D , k = 1 2 k R k z ρ i , φ = eδφ T e ( 18
)
with k z the poloidal wavevector and ω D = -k z T i /(eBR) the drift frequency. With these normalisations, the perturbations are moving in the electron diamagnetic drift direction for ω R > 0. The evolution of the ion normalised perturbed density ñ = δn i /n i , parallel fluid velocity w = δv i /v thi and temperature T = δT i /T i is given by:

ωñ + 2[ñ + T ] + 4[u + k ] w = R L n -2 φ (19) ω w + 4 w + 2[u + k ][ñ + T ] = u ′ -2u -2 k φ (20) ω T + 4 3 ñ + 14 3 T + 8 3 [u + k ] w = R L Ti - 4 3 φ (21)
where the equation for the temperature has been obtained by subtracting Eq. ( 19) to the equation for the pressure. In the three equations above, the terms related to the perturbations are gathered on the left hand side, whereas the ones related to the background are on the right hand side. On the left hand side, one finds, in that order, the temporal derivative of the perturbation, the advection of the perturbation by 1) the curvature and ∇B drift, 2) the Coriolis drift (terms proportional to u) and 3) the parallel motion (terms proportional to k ). On the right hand side, one has the advection of the background by the E × B drift (terms involving the gradients of the background) and finally the advection in velocity space of the background (acceleration) due to the curvature and ∇B drift in the perturbed electric field for Eqs. ( 19) and (21) or to the Coriolis drift and parallel motion in the perturbed electric field for Eq. ( 20).

The Coriolis drift and the parallel dynamics couple parallel velocity fluctuations w to density ñ and temperature T fluctuations. When this coupling is neglected (u = 0 and k = 0), Eqs. ( 19) and ( 21) yield the "standard" picture of the toroidal ITG. This picture is reminded here assuming R/L n = 2 for simplicity. Starting from an initial temperature perturbation T , the advection by the curvature and ∇B drift generates an ion density perturbation (by compression):

ñ = - 2 ω + 2 T (22)
Quasineutrality (δn e = δn i ) and the adiabatic electron response imply that an electrostatic potential perturbation builds up in response to the density perturbation:

φ = ñ ( 23
)
which in turn generates a temperature perturbation by E×B advection in the background temperature gradient:

T = 1 ω + 14/3 R L Ti - 8 3 φ (24)
Depending on the phase shift between T and φ, which is governed by the curvature and ∇B drifts, the perturbation will grow or decay. This is best seen by computing the dispersion relation from Eqs (22)(23)(24):

ω2 + 20 3 ω + 2 R L Ti + 4 = 0 (25)
and solving for the mode growth rate:

γ = k z ρ i v thi R 0 2 R L Ti - 64 9 (26)
which shows that the mode requires a critical temperature gradient R/L Ti > 32/9 to be unstable. When the Coriolis drit and parallel dynamics are considered (finite u and k values), additional density and temperature perturbations are generated by the parallel velocity perturbations. Eqs ( 22) and ( 24) then become:

ñ = - 2 ω + 2 T -2 u + k w (27) T = 1 ω + 14/3 R L Ti - 8 3 φ - 8 3 u + k w (28)
Whether the density and temperature perturbations generated by the Coriolis drift and parallel compression add or subtract to the one generated by the curvature and ∇B compression depends on what mechanism drives the parallel velocity perturbation. This can be determined by computing w with Eq. ( 20). The corresponding dispersion relation has been derived in [18] in the limit of marginaly unstable modes (γ > 0 and γ ≪ |ω R |) and small symmetry breaking (|u + k | ≪ 1). It yields to the following mode growth rate:

γ = γ 2 0 + a 1 [u + k ]u ′ + a 2 [u + k ] 2 (29)
where the assumption R/L n = 2 has been relaxed, γ 0 is the growth rate in the absence of flows (i.e. u = 0 and u ′ = 0, which in this limit imply k = 0) and the coefficients a 1 and a 2 are given by:

a 1 = -16 R/L n -2 R/L n -2/3 , a 2 = - a 1 2 R L n + 26 3 (30)
To get a closed expression for the mode growth rate in the presence of flows, the value of k needs to be specified in Eq. ( 29). The kinetic electron response has been shown to be an essential ingredient in the parallel mode structure response to background flows [18] which makes the quantitative prediction of k well beyond the scope of the simple fluid model used in this section, and in fact beyond the scope of any fluid model assuming an adiabatic electron response. Although inaccurate, the selfconsistent parallel mode structure under the adiabatic electrons approximation is still of interest to investigate the dependencies of k on the flow parameters. A calculation in this spirit has been done in [19] and shows that k ∝ u ′ -2u. Interestingly, this linear dependence of k on u ′ and u with an opposite sign is recovered in gyrokinetic simulations, as shown in Fig. 3 of [20] and in Fig. 2 of the next section. Of course, the respective impact of u and u ′ is not always linked by a factor of 2 and departures from the linear dependence can also be observed. Nevertheless, assuming that k is proportional to u and u ′ with an opposite sign qualitatively captures the main impact of flows on the parallel mode structure. With this assumption, Eq. ( 29) can be written as

γ = γ 2 0 + b u ′ u ′2 + b uu ′ uu ′ + b u u 2 (31)
where b u ′ , b uu ′ and b u are parameter dependent coefficients left unspecified. Eq. ( 31) highlights the coupling of the flow parameters in their impact on the mode growth rate. Retaining centrifugal effects as in [16] would bring several additional terms proportional to u 2 in Eq. ( 29) and therefore modify the value of b u . The qualitative dependence of the mode growth rate on the flow parameters would not be modified however. In summary, the simple model tells us that toroidal rotation (u), parallel flow shear (u ′ ) and the parallel gradient of the perturbed pressure and electrostatic potential (k ) generate parallel velocity perturbations which modify the phase shift between temperature and density perturbation (and hence potential) via the coupling provided by the parallel compression and Coriolis drift. As a consequence, the ITG mode is made more or less stable. The overall effect depends on the respective signs of toroidal rotation u and parallel flow shear u ′ , as emphasised in Eq. (31). The existence of a cross-term proportional to uu ′ implies that the parallel flow shear can have a stabilising effect. In the present section, the impact of E × B shearing was not discussed, but previous analytical works [21,22] show the existence of a similar cross-term between u ′ and γ E in the stability of the linear ITG.

B. Gyrokinetic simulations

The effect of u and u ′ on linear stability is now explored in gyrokinetic simulations for two cases. First, for the Waltz standard case [5], to characterize the parametric dependencies and emphasize the robustness of the mechanisms at play. Second, for the shortfall case at r/a = 0.8 and with zero electron temperature gradient, R/L Te = 0, for which all the non-linear simulations presented in this paper have been performed. The main input parameters for these two cases are gathered in Table I. No E × B shear flow is applied in the linear simulations.

Waltz standard case

For the Waltz standard case, the magnetic equilibrium is described by circular concentric flux surfaces (with a consistent treatment of first order terms in ǫ, see [23]) and only electrostatic perturbations are retained. The impact of the toroidal flow u and parallel flow shear u ′ on the linear mode growth rate is investigated for k θ ρ i = 0.4 and k r ρ i = 0, near the peak of the linear mode growth rate. The variations of the growth rate as a function of u and u ′ are shown in Fig. 1. At zero toroidal rotation, the parallel flow shear is found to be destabilising, consistently with the conventional ITG-PVG coupling picture. However, as soon as the toroidal rotation is finite, the parallel flow shear can be stabilising or destabilising, depending on the sign of uu ′ . Conversely, at finite parallel flow shear, toroidal rotation can be stabilising or destabilising. The linear mode growth rate is maximum for positive uu ′ and minimum for negative uu ′ values, with a saddle point at u = 0 and u ′ = 0. Encouraged by the similarity of this behaviour with the prediction of the fluid model, a fit of the form given by Eq. ( 31) is attempted. The coefficients b u ′ , b uu ′ and b u are calculated from the growth rate values of the five cases given in Table II. The combination of proves to be of rather good quality: the relative difference between the fitted and actual growth rates is less than 5% over the considered (u, u ′ ) domain. The maximum difference is obtained in the top-left and bottom right corners (up to 5% for uu ′ < -0.6) and in the top-right and bottom-left corners (up to 3.5% for uu ′ > 0.6). In the remaining part (|uu ′ | < 0.6) the relative difference is less than 2%. The comparatively worst fit accuracy in the corners of the domain (|uu ′ | > 0.6) is due to the fact that the assumptions of the fluid model are not fully valid. In the top right and bottom left corners, |u + k | reaches up to 0.3, which is not particularly negligible compared to unity. In the top left and bottom right corners, a stronger dependence of the parallel vector on u and u ′ is obtained, Fig. 2, whereas a constant linear dependence on u and u ′ over the whole domain was assumed in Eq. (31). In Fig. 2, the effective parallel wavevector is computed as The relative impact of the toroidal rotation and parallel flow shear contributions (b u and b u ′ , respectively) and of the cross-term (b uu ′ ) on the mode growth rate are shown in Fig. 3 as a function of the local inverse aspect ratio ǫ, the safety factor q, the magnetic shear ŝ, the poloidal wavevector k θ ρ i , the electron to ion temperature ratio T e /T i and the normalised ion temperature gradient R/L Ti . The values of b u , b uu ′ and b u ′ are computed from the five (u,u ′ ) pairs given in Table II and normalised to the growth rate value without flow γ 0 = γ u=0,u ′ =0 . Additional 2D scans were performed to check that the fit accuracy remains acceptable over the range of scanned parameters. The relative difference between the fitted and actual growth rate values is less than 5%, except for q > 4 and k θ ρ i > 0.6 where it reaches about 10%. As for the reference case, the largest differences are obtained for uu ′ < -0.6 and are correlated to large values of k eff . Over the whole parameter range, with the exception of the R/L Ti scan that will be discussed separately, b u is systematically found to be negative (ITG stabilisation for increasing |u|), b u ′ positive (ITG destabilisation for increasing |u ′ |) and b uu ′ positive (ITG destabilisation for increasing uu ′ ). The impact of these three terms on the mode growth rate is of comparable magnitude (in normalised units u ′ is typically 10 times larger than u). The strongest stabilisation by toroidal rotation is obtained at high q, high k θ ρ i and low T e /T i . The destabilisation by parallel flow shear is particularly enhanced at high ŝ values and to a lesser extent at low ε and q values and at low and high k θ ρ i values. The magnitude of the crossterm mostly depends on ε, q an k θ ρ i and increases with these parameters. The R/L Ti scan has the specific feature that a transition from a Trapped Electron Mode (TEM) to an ITG is obtained when increasing R/L Ti . Interestingly the impact of the flows on the mode sta- 3: Normalised coefficients bu/γ0 (red), 10b uu ′ /γ0 (green) and 100b u ′ /γ0 (blue) as a function of the local inverse aspect ratio ǫ, the safety factor q, the magnetic shear ŝ, the poloidal wavevector k θ ρi, the electron to ion temperature ratio Te/Ti and the normalised ion temperature gradient R/LT i . The scans are performed around the reference Waltz standard case. Dashed and full lines indicate that the dominant instability is a TEM or an ITG, respectively.

u 0 0.3 0 0.3 -0.3 u ′ 0 0 3 3 3
k eff = k |φ k | 2 |φ k | 2 (32
bility is reversed at the ITG/TEM transition and the b coefficients all change signs. The magnitude of these coefficients gets particularly large close to the mode transition. Just around the transition, i.e. for R/L Ti ∼ 6 -7, changing the flow parameters can even trigger a change of the dominant instability (this is why the calculation of the b coefficients was not possible in this region). Deep in the TEM regime, when R/L Ti approaches zero, the impact of u and u ′ on the mode stability gets vanishingly small.

Shortfall case

The linear runs for the shortfall cases retain electrostatic and electromagnetic perturbations (φ and A ), collisions (pitch-angle scattering only) and flux surface shaping. The Miller parameterisation of the magnetic equilibrium at r/a = 0.8 is given in Table III. The updown asymmetry of the equilibrium governed by the parameter ∂Z mil /∂r was set to zero to focus on symmetry breaking driven by the flows (using the experimental value ∂Z mil /∂r = 0.0708 changes the mode growth rate by about 1-2%). The most unstable mode is an ITG over the whole spectrum and reaches its maximum growth rate around k θ ρ i = 0.4 and k r ρ i = 0. The impact of the toroidal rotation and parallel flow shear at about 3 times larger, which is mostly due to the higher q and ŝ values, and the change of slope of k eff as a function of u and u ′ occurs earlier when decreasing uu ′ . The largest k eff values are mirrored by a very strong distortion of the eigenfunctions and a shift of the maximum from the low field side midplane, see Fig. 6 , which makes the failure of a 1-point fluid model in the strong ballooning approximation understandable. As mentioned in the introduction, the importance of the synergy between the plasma flow components for turbulence stabilisation has been realised in the frame of a benchmark exercise based on the DIII-D shortfall case. The focus of the benchmark was on ion-scale turbulence which is why all the non-linear simulations presented in this paper are with R/L Te = 0. Non-linear simulations at finite R/L Te have not been performed yet, but some insight can be gained by looking at the impact of this parameter on the linear stability results. The variation of the linear growth rate with u and u ′ has therefore been explored at R/L Te = 13.6 (nominal value) and R/L Te = 18. At R/L Te = 13.6, the most unstable mode is an ITG over most of the domain and the variation of the flow parameters leads to about ±20% changes in the linear growth rate. Qualitatively, the growth rate changes with u and u ′ are between the results obtained for the Waltz standard case and the shortfall case at R/L Te = 0, i.e. more negative uu ′ values still tend to be stabilising but not as much as for the Waltz standard case. In the top-left and bottom-right corners, the stabilisation of the ITG leads to a transition to a TEM (TEM dominated points are indicated by black stars on the left plot of Fig. 7). At R/L Te = 18, the most unstable mode is a TEM over the whole domain and the effect of flows gets reversed as in the R/L Ti scan of Fig. 3. The variation of the growth rate with u and u ′ also gets smaller and is about ±10% with respect to the case without flows. When the linear spectrum is ITG dominated at low k θ ρ i and TEM dominated at high k θ ρ i , as for the shortfall case at nominal R/L Te , a partial compensation of the flow effect could therefore potentially occur in non-linear simulations. This point would deserve further investigation and is left for future studies.

Finally, the inclusion of centrifugal effects for these cases with moderate toroidal rotation (ion Mach number |u| < 0.3) do not strongly modify the general picture, as expected. For the highest uu ′ values, centrifugal effects are nevertheless not totally negligible and tend to increase by some 10 -15% the destabilising effect of parallel flow shear, Fig. 8.

IV. NON-LINEAR SIMULATIONS A. Simulations set-up

The non-linear simulations focus on the shorfall case at R/L Te = 0 and investigate the impact of plasma flows by varying u, u ′ and γ E around their nominal values. The simulations are performed for kinetic electron and deuterium species with their actual mass ratio m e /m i = 2.72 × 10 -4 . The 5D computational domain is discretised with N s = 32 points in the direction parallel to the magnetic field line (finite differences), N k θ × N kr = 21 × 339 wavevectors in the perpendicular plane (spectral decomposition) and N µ × N v = 16 × 64 points in velocity space. The 21 poloidal wavevectors range from k θ ρ i = 0 to 1.3626 (for the reference case, the linear growth rate is zero for k θ ρ i > 0.8) and the radial wavevectors extend up to k r ρ i = ±24.2457. For the runs presented in this section, the electrostatic potential perturbation is more than 600 times smaller at the maximum k r ρ i than at k r ρ i = 0, which guarantees that the numerical scheme used to describe E × B shearing (wavevector remapping) is applicable. To minimize the simulation duration, the E × B shearing is only applied after the non-linear overshoot. The simulations are run until a converged time average is obtained, starting after the overshoot or after the E × B shearing is turned on. Depending on the case, the required interval for the time average ranges between 400 and 1000R 0 /v thi . A typical case is illustrated in Fig. 9, showing the time evolution of the ion and electron heat fluxes, the convergence of the time averaged fluxes and the time average of < kr |φ(k r , k θ )| 2 > as a function of the poloidal wavevector. The electrostatic heat fluxes are defined as:

Q r s = 1 2 m s v 2 δv E • ∇r δf s dv ( 33
)
where δv E and δf s are the perturbed E × B velocity and distribution function, respectively, and normalised as:

Q N,s = Q r s n s T s ρ 2 * v thi (34)
The heat flux due to the magnetic flutter is less than 1% in all the cases considered here and therefore neglected. For comparison with other works based on the shortfall case, e.g. [11,24,25], the dimensional fluxes can be calculated from the following normalising factors: 526.43 W.m -2 for the ion heat flux and 448.42 W.m -2 for the electron heat flux. The radial derivative of the plasma volume, sometimes used to make the fluxes independent of the radial coordinate definition, is ∂V/∂r = 44.615 m 2 at r/a = 0.8.

B. Stabilisation by toroidal rotation

The impact of toroidal rotation is first explored at finite parallel flow shear, u ′ = -1.2584, and zero E × B shear flow. Increasing the toroidal rotation from negative to positive values has first almost no effect and then results into a nearly monotonic decrease of the non-linear ion heat flux, see Fig. 10. The heat flux is about 30% smaller at u = 0.3 compared to u = -0.15. Qualitatively similar results are obtained at finite E × B shear flow, γ E = -0.2106, but with an enhanced effect of toroidal rotation and a stronger stabilisation at u = -0.3. In the presence of E × B shearing, increasing the toroidal rotation reduces the heat flux by more than a factor of two. The behaviour of the electron heat flux closely follows the one of the ion heat flux (i.e. Q e /Q i remains almost constant). Surprisingly, at the nominal toroidal rotation u = -0.1549, turning on the E × B shearing has almost no impact on the ion heat flux, whereas a decrease of the transport level is usually expected and in fact observed at the other toroidal rotation values. Looking at the ion heat flux spectrum reveals that E × B shearing has a strong impact, as expected, but that the reduction of the heat flux at k θ ρ i < 0.25 is compensated by an increase of the heat flux at 0.25 < k θ ρ i < 0.4. For u = -0.1549, left plot of Fig. 11, this compensation is almost perfect and no effect is seen on the total heat flux. For the other cases, as for instance u = 0.1549 in the right plot of Fig. 11, the heat flux reduction at low k θ ρ i is larger than the increase at high k θ ρ i and overall the E × B shearing results in a smaller heat flux. The peak of the time averaged electrostatic potential spectrum, defined as

k α ρ peak i = kr,k θ k α |φ(k r , k θ )| 2 kr,k θ |φ(k r , k θ )| 2 (35)
with k α standing for k r or k θ , does not depend significantly on the toroidal rotation but is shifted towards higher k r ρ i and k θ ρ i when E × B shearing is applied, Fig. 12. The linear mode growth rate for the wavevectors around the peak location is shown as a function of the toroidal rotation in Fig. 13. It is interesting to note the similarity between the response of the linear growth rate at the peak of the non-linear spectrum (Fig. 13) and the non-linear fluxes (Fig. 10) to an increase of toroidal rotation. This similarity suggests that the non-linear behaviour finds its root in the linear stability response to toroidal rotation. Quantitatively, the impact on the nonlinear fluxes is however much larger than on the linear mode growth rate (more than 5 times). The reason of the stronger impact observed in non-linear simulations is not understood yet.

C. Impact of the parallel and perpendicular flow shears at finite toroidal rotation

We now turn our attention to the impact of the parallel and perpendicular flow shears in the presence of a finite toroidal rotation. The toroidal rotation is kept fixed at its nominal value and the parallel flow shear is scanned for three values of the E × B shear flow. As shown in Fig. 14, the heat flux response to an increase of u ′ strongly depends on the value of γ E . At the nominal value, γ E = -0.2106, the heat flux decreases linearly with increasing u ′ . This reduction is large and the heat flux is almost reduced by a factor of 3 by simply reversing the parallel flow shear. In contrast, at γ E = 0.2106, there is almost no dependence of the heat flux on the parallel flow shear. The response at γ E = 0 lies between this two extremes and a moderate stabilisation with increasing parallel flow shear is obtained. Focusing on the three cases consistent with the purely toroidal flow assumption (thick symbols in Fig. 14), the heat flux level is observed to be asymmetric with respect to a reversal of the parallel and perpendicular flow shears. This asymmetry is due to the finite toroidal rotation and results in a decrease of the ion heat flux by about 35% when the parallel and perpendicular flow shears are reversed. As in the previous section, the peak of the time averaged electrosatic potential is affected by E × B shearing. The peak position moves towards more negative k r ρ i values when γ E increases and to more positive k θ ρ i values at finite |γ E |. The parallel flow shear has a an effect on k r ρ peak i comparable in magnitude, but in the opposite direction, whereas it barely affects k θ ρ peak i . These dependencies are summarised in Fig. 15. The k r ρ i shift of the electrostatic potential spectrum caused by E × B shearing is known to play a key role in the non-linear turbulence saturation [26] and is associated to a strong modification of the parallel structure of the electrostatic potential [27]. In linear simulations, the parallel asymmetry of the electrostatic potential generated by a finite k r ρ i couples to u ′ and makes the impact of u ′ on the mode growth rate depend on the value of k r ρ i , see Fig. 16. The interplay between γ E and u ′ observed in non-linear simulations could therefore be due to the modification of the parallel mode structure caused by the k r ρ i shift. In fact, when the linear growth rate is computed at the non-linear electrostatic potential peak, hence with a finite k r ρ i , the heat flux variations with respect to γ E and u ′ obtained in non-linear runs qualitatively follow the linear growth rate variations, see Figs 14 and16. As for toroidal rotation, the variation of the non-linear heat fluxes is much larger than the one of the linear growth rate.

Finally, the ion heat flux is shown in Fig. 17 as a function of the time averaged electrostatic potential perturbation amplitude, |φ| noZF , for all the non-linear simulations. The potential amplitude is calculated without the radial mode contributions:

|φ| noZF = k θ =0 |φ(k r , k θ )| 2 dk r dk θ (36)
The linear increase of the ion heat flux with |φ| noZF shows that the main impact of flows is not a modification of the cross-phase between temperature and electrostatic potential fluctuations, but rather a change of the nonlinear saturation level.

V. CONCLUSIONS

The impact of toroidal rotation u and parallel flow shear u ′ on the linear stability of the toroidal ITG has been investigated in the frame of a simple fluid model. The key role played by the Coriolis drift and the parallel dynamics in coupling density and parallel velocity perturbations has been highlighted and the linear growth rate shown to have a dependence on |u|, |u ′ | and a cross-term proportional to uu ′ . As a consequence, the toroidal rotation and the parallel flow shear can have a stabilising or a destabilising effect, depending on their respective sign and of the magnitude of the cross-term. Guided by the fluid model results, linear gyrokinetic simulations have been performed to determine the scaling of the crossterm with respect to plasma parameters. The strongest growth rate dependence on the cross-term is obtained at high safety factor q, high local inverse aspect ratio ε, high poloidal wavevector k θ ρ i and close to the ITG/TEM transition. Interestingly, the growth rate dependence on u and u ′ is found to reverse sign at the ITG/TEM transition. The strong interplay between u and u ′ has been confirmed in non-linear simulations based on the DIII-D shortfall case at r/a = 0.8 and with R/L Te = 0. The flow effect on the non-linear fluxes qualitatively follows the linear growth rate variations but with a much stronger amplitude. Adding E × B shearing brings another parallel symmetry breaking term into the game, which can enhance or reduce the turbulence level depending on its sign with respect to u and u ′ . Qualitatively, the strongest impact is obtained when the different terms conjugate to generate a large parallel symmetry breaking (for instance u > 0, u ′ < 0 and γ E > 0). The most striking effect observed in the non-linear simulations is the factor of three reduction of the ion heat flux obtained when reversing the parallel flow shear and keeping u and γ E at their nominal values. The relevance of this mechanism in the experiments depends on the ability to decouple the flow parameters. In the limit of purely toroidal flows, u ′ and γ E are strictly proportional. This proportionality is broken when the poloidal rotation and the pressure gradient make a significant contribution to the radial electric field which is typically the case in transport barriers or close to the plasma edge. We caution again the reader that, formally, this effect is first order in ρ * and should be described by a theory derived up to the same order. Decoupling u from the flow shear is to some extent easier to achieve: increasing the ripple magnitude in the presence of central NBI heating [28] or moving the X-point location from the inboard to outboard side of the vacuum chamber [29] provides an effective way to change the value of the toroidal rotation at relatively constant parallel flow shear by acting on the boundary condition. There are therefore experimental conditions in which the sign dependent interaction of flows in turbulent transport stabilisation evidenced in this theoretical study could be significant. It may in particular be a relevant candidate to interpret the unexplained dependence of the E × B quenching strength on toroidal rotation in JET [30] or the dependence of ITBs strength on the direction of toroidal rotation in JT60U [31].

where v is the parallel velocity coordinate in the rotating frame and Ω ⊥ = Ω-(Ω•b)b with b the unit vector in the direction of the magnetic field B and Ω = RΩ∇R × ∇ϕ. There are two important differences between the Coriolis and centrifugal forces. First, the Coriolis force depends on the parallel velocity coordinate v , whereas the centrifugal force does not. Starting from a density perturbation extended along the magnetic field line and with zero average (fluid) velocity, the Coriolis drift will therefore split the particles according to their parallel velocity and generate parallel fluid velocity perturbations. This is the reason why the Coriolis drift plays a central role in the turbulent transport of parallel momentum [32]. Second, unlike the Coriolis force, the centrifugal force contributes to the parallel acceleration (because it is not perpendicular to the parallel motion). The corresponding modification of the parallel force balance gives rise to poloidal asymmetries in the equilibrium density and electrostatic potential [14]:

n s = n R0,s exp [-E s /T s ] (A3)
where n R0,s is the density of species s at the reference major radius R 0 and the energy E s is defined as

E s = Z s eΦ (1) - 1 2 m s Ω 2 (R 2 -R 2 0 ) (A4)
The poloidally varying electrostatic potential Φ (1) arises to balance the charge separation caused by the centrifugal force acting on plasma constituents of different mass and is solution of Σ s Z s n s = 0.

In GKW, the impact of the centrifugal force on the equilibrium is taken into account by specifying the Maxwellian background distribution function as:

F M,s = exp - Z s eΦ (1) T s - m s [Ω + ω f s ] 2 2T s (R 2 -R 2 0 ) × n R0,s π 3/2 v 3 ths exp - m s [v -(F/B)ω f s ] 2 2T s - µ s B T s (A5)
where µ s = m s v 2 ⊥ /B 0 and ω f s = ω s -Ω is the toroidal angular frequency of species s in the rotating frame. By assumption, ω f s (r 0 ) = 0 for all species and the Maxwellian background (entering in the Poisson equation [9]) is then given by:

F M,s = n R0,s π 3/2 v 3 ths exp - m s v 2 2T s - µ s B T s - E s T s (A6)
The radial gradient of the Maxwellian enters the gyrokinetic equation, but part of it cancels with other terms. The remaining contribution is noted ∇ p F M,s in GKW reference papers. Centrifugal effects require the following contribution to be added to ∇ p F M,s :

∇ p F M,s | cf = E s T 2 s ∂T s ∂r + m s Ω T s (R 2 -R 2 0 ) ∂ω f s ∂r F M,s (
A7) In addition to the modification of the equilibrium, the effect of the inertial forces on the equations of motion is taken into account by considering the following extra contributions to the gyro-centers velocity and parallel acceleration:

dX dt co = v co (A8) dX dt cf = v cf + b × ∇Φ (1) B (A9) m s v dv dt cf = - dX dt • ∇E s (A10)
The centrifugal contribution to the parallel acceleration, Eq. (A10), enhances the low field side trapping and is responsible for an enhanced destabilisation of the trapped electron mode at high toroidal rotation [9]. To summarize, GKW takes into account the inertial effects due to the lowest order neoclassical flows, i.e. a purely toroidal rotation of constant angular frequency on a flux surface that is independent of the species. Practically, this is achieved by formulating and solving the equations in the frame rotating at the plasma angular frequency. In this frame, the gyrokinetic equations are modified as follows. The Coriolis force results in the addition of the Coriolis drift to the gyro-center velocity.

The centrifugal force results in the modification of the background equilibrium (poloidal asymmetries), in the addition of the centrifugal drift and of the E × B drift due to the first order electrostatic potential to the gyrocenter velocity and to the corresponding contribution of these two terms to the parallel acceleration. The parallel flow shear is implemented as follows in the gradient of the Maxwellian background:

∇ p F M,s = F M,
+ m s v F BT s + m s Ω T s [R 2 -R 2 0 ] ∂ω f s ∂r (B1)
This implementation assumes that the parallel flow shear is due to the projection of the shear in the toroidal flow and that the species toroidal angular frequency in the rotating frame ω f s is constant on a flux surface. The possibility to have a species dependent toroidal angular frequency gradient ∂ω f s /∂r is retained.

E × B flow shear

In the rotating frame, the lowest order E × B flow, v f E , is zero at r = r 0 , by definition of the frame angular frequency. The possibilty to have a finite shear in the E × B velocity is retained and included in the code by expanding v f E around r = r 0 :

v f E (r) = (r -r 0 ) b × ∇r B ∂ 2 Φ f ∂r 2 r0 (B2)
with Φ f the lowest order electrostatic potentiel in the rotating frame. The advection of perturbed quantities by the background E × B sheared flow (i.e. terms appearing in the gyrokinetic equation as v f E •∇g, with g a perturbed quantity) is then treated in Fourier space by a periodic remapping of the radial wavectors [10].

Shear flow specification for the modeling of experimental plasmas

Starting from the parallel fluid velocity in the rotating frame: and by using Eq. ( 14) to specify γ E , the departure from a purely toroidal flow is approximately taken into account in the relationship between u ′ s and γ E . What is neglected is the variation of the parallel flow shear on a flux surface due to ûθ,s and ω p,s .

V f ,s = ûθ,s B + F B [s b ω p,s + s b ω Φ -Ω] , (

  III. IMPACT OF u AND u ′ ON THE LINEAR STABILITY OF THE TOROIDAL ITGA. Insights from a simplified fluid model

FIG. 1 :

 1 FIG.1: Linear mode growth rate at k θ ρi = 0.4 as a function of the toroidal rotation u and parallel flow shear u ′ for the Waltz standard case.

FIG. 2 :

 2 FIG.2: Effective parallel wave vector at k θ ρi = 0.4 for the Waltz standard case as a function of the toroidal rotation u for increasing parallel flow shear u ′ (from -3 in blue to +3 in red, with steps of 0.75). Note the change of slope for uu ′ < -0.6.

  FIG.3: Normalised coefficients bu/γ0 (red), 10b uu ′ /γ0 (green) and 100b u ′ /γ0 (blue) as a function of the local inverse aspect ratio ǫ, the safety factor q, the magnetic shear ŝ, the poloidal wavevector k θ ρi, the electron to ion temperature ratio Te/Ti and the normalised ion temperature gradient R/LT i . The scans are performed around the reference Waltz standard case. Dashed and full lines indicate that the dominant instability is a TEM or an ITG, respectively.

this wavector is shown in Fig. 4 .FIG. 4 :FIG. 5 :

 445 FIG.4: Linear mode growth rate at k θ ρi = 0.4 as a function of the toroidal rotation u and parallel flow shear u ′ for the shortfall case.

FIG. 6 :

 6 FIG.6: Parallel structure of the electrostatic potential (real part on the left and imaginary part on the right) for the shortfal case at k θ ρi = 0.4, shown for u = 0 and u ′ = 0 (blue curve), u = 0.3 and u ′ = 3 (green curve) and u = -0.3 and u ′ = 3 (red curve). The parallel coordinate s is zero at the low field side midplane and ±0.5 at the high field side midplane.

FIG. 7 :

 7 FIG. 7: Linear mode growth rate at k θ ρi = 0.4 as a function of the toroidal rotation u and parallel flow shear u ′ for the shortfall case at R/LT e = 13.6 (left plot) and R/LT e = 18 (right plot).

16 FIG. 8 :

 168 FIG.8: Left plot: linear mode growth rate at k θ ρi = 0.4 as a function of the toroidal rotation u and parallel flow shear u ′ for the shortfall case with centrifugal effects included. Right plot: relative difference between the linear mode growth rate obtained with (γ CF ) and without (γ) centrifugal effects.

FIG. 9 :

 9 FIG. 9: Shortfall case with γE = 0. Left: time evolution of the ion (blue) and electron (red) normalised heat flux. The vertical black dashed line indicates the starting point for the time averages. Top right: time averaged ion and electron heat fluxes for increasing time average interval. Bottom right: time averaged electrosatic potential spectrum (for the largest time interval).

FIG. 10 :

 10 FIG.10: Time averaged normalised ion heat flux as a function of toroidal rotation for non-linear simulations of the shortfall case with γE = 0 (blue stars) and γE = -0.2106 (red squares). The parallel flow shear is kept at its nominal value u ′ = -1.2584.

FIG. 11 :FIG. 12 :

 1112 FIG. 11: Time averaged ion heat flux spectra, integrated over the kr direction, as a function of k θ ρi for γE = 0 (blue) and γE = -0.2106 (red). Left plot: u = -0.1549. Right plot: u = 0.1549.

FIG. 13 :FIG. 14 :

 1314 FIG.13: Linear mode growth rate as a function of toroidal rotation u for the shortfall case at k θ ρi = 0.2044 (blue) and k θ ρi = 0.2725 (red). Dashed lines correspond to krρi = -0.05, full lines to krρi = 0 and dashed-dotted lines to krρi = 0.05.

FIG. 15 :

 15 FIG. 15: Position of the peak of the electrostatic potential as a function of the parallel flow shear u ′ in non-linear simulations of the shortfall case with γE = -0.2106 (red squares), γE = 0 (blue stars) and γE = 0.2106 (purple triangles).

FIG. 16 :

 16 FIG.16: Linear mode growth rate as a function of the radial wave vector krρi for the shortfall case at k θ ρi = 0.2044 (blue) and k θ ρi = 0.2725 (red) for u ′ = -1.2584 (full line) and u ′ = 1.2584 (dashed line).

FIG. 17 :

 17 FIG. 17: Time averaged normalised ion heat flux as a function of the total electrostatic potential perturbation amplitude without the contributions of radial modes, for all the non-linear simulations based on the shortfall reference case. Cases with and without E × B shearing are shown in red and blue, respectively.

Appendix B: Flow shear in GKW 1 .

 1 Parallel flow shear

TABLE II

 II

: Toroidal flow and parallel flow shear values used to compute the coefficients in Eq.

(31)

.

Acknowledgments

The authors would like to thank C. Bourdelle, J. Citrin, N. Dubuit, D. Escande, T. Görler, W. Guttenfelder, C. Holland, M. Honda, E. Narita and G. Staebler for many stimulating discussions from which this work has greatly benefited. The simulations were carried out using HPC resources from GENCI-IDRIS (Grant i2015057398), from the project Equip@Meso (ANR-10-EQPX-29-01) at Aix-Marseille and from the HELIOS supercomputer system at International Fusion Energy Research Centre, Aomori, Japan. This work was partially carried out with funding from the RCUK Energy Programme [grant number EP/I501045].

Appendix A: Inertial effects in GKW

In the rotating frame, the inertial effects appear through the Coriolis and centrifugal drifts:

and


Demian Wassermann

Dorian Mazauric

Guillermo Gallardo-Diez

Rachid Deriche

Extracting the Core Structural Connectivity Network: Guaranteeing Network Connectedness Through a Graph-Theoretical Approach

Keywords:

We present a graph-theoretical algorithm to extract the connected core structural connectivity network of a subject population. Extracting this core common network across subjects is a main problem in current neuroscience. Such network facilitates cognitive and clinical analyses by reducing the number of connections that need to be explored. Furthermore, insights into the human brain structure can be gained by comparing core networks of different populations. We show that our novel algorithm has theoretical and practical advantages. First, contrary to the current approach our algorithm guarantees that the extracted core subnetwork is connected agreeing with current evidence that the core structural network is tightly connected. Second, our algorithm shows enhanced performance when used as feature selection approach for connectivity analysis on populations.

Introduction 4 Conclusions

Our pure extrinsic based parcellation has good agreement with anatomical and functional parcellations from the literature. Particularly, the motor and sensory cortex appear to be found. Having the tractograms in a vectorial space allowed us to work e ciently with them and to create a populationrepresentative parcellation.

Conclusions

Our pure extrinsic based parcellation has good agreement with anatomical and functional parcellations from the literature. Particularly, the motor and sensory cortex appear to be found. Having the tractograms in a vectorial space allowed us to work e ciently with them and to create a populationrepresentative parcellation.

Conclusions

Our pure extrinsic based parcellation has good agreement with anatomical and functional parcellations from the literature. Particularly, the motor and sensory cortex appear to be found. Having the tractograms in a vectorial space allowed us to work e ciently with them and to create a populationrepresentative parcellation. Our pure extrinsic based parcellation has good agreement with anatomical and functional parcellations from the literature. Particularly, the motor and sensory cortex appear to be found. Having the tractograms in a vectorial space allowed us to work e ciently with them and to create a populationrepresentative parcellation. Our pure extrinsic based parcellation has good agreement with anatomical and functional parcellations from the literature. Particularly, the motor and sensory cortex appear to be found. Having the tractograms in a vectorial space allowed us to work e ciently with them and to create a populationrepresentative parcellation. Isolating the common core structural connectivity network (SCN) of a population is an important problem in current neuroscience [START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF][START_REF] Gong | Mapping Anatomical Connectivity Patterns of Human Cerebral Cortex Using In Vivo Diffusion Tensor Imaging Tractography[END_REF]. This procedure facilitates cognitive and clinical studies based on Diffusion MRI [e.g. 1, 5] by increasing their statistical power through a reduction of the number of analyzed structural connections. We illustrate this process in Fig. 1. Furthermore, recent evidence indicates a core common network exists in human and macaque brains and that it is tightly connected [START_REF] Bassett | Task-Based Core-Periphery Organization of Human Brain Dynamics[END_REF]. In this work we develop, for the first time, a group-wise core SCN extraction algorithm which guarantees a connected network output. Furthermore, we show the potential of such network to select gender-specific connections through an experiment on 300 human subjects.

The most used population-level core SCN extraction technique [START_REF] Gong | Mapping Anatomical Connectivity Patterns of Human Cerebral Cortex Using In Vivo Diffusion Tensor Imaging Tractography[END_REF] is based on an effective statistical procedure to extract a population SCN: 1) computing, for each subject, a connectivity matrix using a standardised parcellation; and 2) extracting a binary graph by analysing each connection separately and rejecting hypothesis is not in the population. The resulting graph can be a set of disconnected subgraphs. This is problematic, recent studies have shown the core network to be tightly connected [START_REF] Bassett | Task-Based Core-Periphery Organization of Human Brain Dynamics[END_REF]. However, extracting connected group-wise core SCN is far from simple: an algorithm to find the largest core network of a population cannot find an approximated solution in polynomial time.

In this work, we propose a graph-theoretical algorithm to obtain the connected core SCN of a subject sample. Our approach guarantees a connected core SCN, agreeing with novel evidences on structural connectivity network topology [e.g. 2]. We start by proving that we can formulate the problem such that core network extraction is NP-Complete in general but in our case we find an exact polynomial time algorithm to perform the extraction. Finally, we show that our algorithm outperforms that of Gong et al. [5] as a tool for selecting regressors in SC. For this, we use 300 subjects from the HCP database and comparing the performance of the networks obtained with both algorithms to predict connectivity values from gender in a subsection of the core network.

Definitions, Problems, and Contributions

A first approach to core sub-network identification can be derived from the binary connectivity model. In this model the cortical and sub-cortical regions are common across subjects and what varies is whether these regions are connected or not [START_REF] Gong | Mapping Anatomical Connectivity Patterns of Human Cerebral Cortex Using In Vivo Diffusion Tensor Imaging Tractography[END_REF]. Using this approach, a sample of human brain connectivity of a given population can be represented by k ≥ 1 graphs

G 1 = (V, E 1 ), . . . , G k = (V, E k ).
In this formalism each graph G i corresponds to a subject and, in accordance with Gong et al. [5], the vertices, V , stable across subjects are cortical and sub-cortical regions and the edges E i are white matter bundles connecting those regions. Note that all graphs have the same ordered set of nodes. A first approach to compute, or approximate, the core sub-network of the population sample consists in finding the core SCN graph G * = (V * , E * ) such that G * and every G i has some quantitative common properties, where V * ⊆ V . In this article, we model the difference between the core SCN, G * , and the subject ones, G i , by a function f λ . This function measures the difference between the sets of edges (and the sets of non-edges) of the core network and those of the subjects: difference between the set of edges of G and the set of edges of G i modulated by the parameter λ. In the following, we will refer to f λ (G * , G i ) as the difference threshold of a core sub-network G * wrt G i . Note that if λ = 1, we only consider edges excluded from the core network, |{e ∈ E, e / ∈ E(G i [V * ])}|, and if λ = 0, we only consider edges included in the core network, |{e / ∈ E, e ∈ E(G i [V * ])}|. In Definition 1, we formalize the problem of computing the core sub-network as a combinatorial optimization problem:

f λ (G * , G i ) = λ|{e ∈ E, e / ∈ E(G i [V * ])}| + (1 -λ)|{e / ∈ E, e ∈ E(G i [V * ])}|, where λ ∈ [0, 1] and G i [V * ] is the subgraph of G i induced by the set of nodes V * ⊆ V
Definition 1 (Core sub-network problem). Let G 1 = (V, E 1 ), . . . , G k = (V, E k ) be k ≥ 1 undirected graphs. Let λ be any real such that λ ∈ [0, 1].
Let n ≥ 0 be any integer. Then, the core sub-network problem consists in computing a connected graph G * = (V * , E * ) such that |V * | ≥ n and such that the sum of the difference thresholds

k i=1 f λ (G * , G i ) is minimum.
c) represent G 1 , G 2 , G 3 , respectively. Figure 2 (d) is a solution G * = (V * , E * ) when n = 5. Indeed, we have f λ (G * , G 3 ) = 1
2 because the difference between G * and G 3 is the edge connecting nodes 2 and 5 or the two element set {2, 5}; we have f λ (G * , G 2 ) = 2 because the difference between G * and G 2 is the four edges {1, 5}, {1, 4}, {3, 4}, and {4, 5}; and we have f λ (G * , G 1 ) = 1 because the difference between G * and G 1 is the two edges {3, 5} and {4, 5}. We

get f λ (G * , G 1 ) + f λ (G * , G 2 ) + f λ (G * , G 3 ) = 7
2 . In the rest of this section we state our main contribution, an optimal polynomial time exact algorithm for the core sub-network problem if the number of nodes is sufficiently large (optimal means here that there is no exact algorithm with better complexity). Solving the problem in Def. 1, is hard: it can be proved that given an integer n ≥ 0 and a real number δ ≥ 0, then the decision version of the SCN problem is NP-complete even if k = 2. However, focusing on the problem of minimizing f λ we obtain a polynomial time algorithm for SCN extraction.

The main point of this work is to present an algorithm for the core graph extraction and assess its potential for clinical and cognitive studies. Even if the problem is very difficult to solve in general, we design our polynomial time core subnetwork extraction algorithm and show that it is optimal, when we focus on the problem of minimizing the difference threshold and when the number of nodes of the core sub-network is large. In the following, we aim at proving Theorem 1. Consider k ≥ 1 undirected graphs G 1 = (V, E 1 ), . . . , G k = (V, E k ) and consider any real number λ such that λ ∈ [0, 1]. Let us define some notations and auxiliary graphs. Let G = (V, V × V ) be the completely connected cortical network graph with V × V the set all pairs from the elements from V . We define two edge-weighting functions w 0 and w 1 :

Theorem 1. Consider k ≥ 1 undirected graphs G 1 = (V, E 1 ), . . . , G k = (V, E k )
w 0 (e) = (1 -λ)|{i, e ∈ E i , 1 ≤ i ≤ k}|, w 1 (e) = λ|{i, e / ∈ E i , 1 ≤ i ≤ k}|. (1)
Intuitively, w 0 (e) represents the cost of not adding the edge e in the solution and w 1 (e) represents the cost of adding the edge e in the solution. From this, we define the graph induced by the set of edges to keep in the core subnetwork

G 1 = (V, E 1 ), E 1 = {e | w 1 (e) ≤ w 0 (e)} ⊆ V × V. (2)
If G 1 is a connected graph, then it is an optimal solution. Otherwise, we have to add edges in order to obtain a connected graph while minimizing the cost of such adding. To add such edges we define a graph representing the fully connected graph where each node represents a maximal connected component:

G cc = (V cc , E cc ) with V cc = {u 1 , . . . , u t } and E cc = V cc × V cc , (3)
where cc(G 1 ) = (cc 1 (G 1 ), . . . , cc t (G 1 )) is the t maximal connected components of G 1 . Then, to select which maximal connected components to include in our core subnetwork graph, we define a weight function w cc :

w cc ({u i , u j }) = min v∈V (cci(G1)),v ∈V (ccj (G1))
w 1 (e)-w 0 (e), where 1 ≤ i < j ≤ t. (4)

We formally prove in Lemma 1 that the problem of obtaining a minimum connected graph from G 1 , that is solving the core sub-network problem when n = |V |, consists in computing a minimum spanning tree of G cc .

Lemma 1. The core sub-network problem when n = |V | is equivalent to compute a minimum spanning tree of G cc = (V cc , E cc ) with weight function w cc .

Proof. The core sub-network problem when

n = |V | consists in computing a graph G * = (V * , E * ) such that V * = V and δ * = k i=1 f λ (G * , G i ) is mini- mum. Consider the graph G 1 = (V, E 1 ) previously defined. Observe that: δ * ≥ e∈V ×V min(w 0 (e), w 1 (e)) = e∈E 1 w 1 (e) + e∈V ×V |e / ∈E 1 w 0 (e). Indeed, for every pair of nodes v, v of V , either we set {v, v } ∈ E * if w 1 ({u, v}) ≥ w 0 ({v, v }) or we set {v, v } / ∈ E * . Thus, if G 1 is a connected graph, then G * = G 1 is
an optimal solution such that δ * = e∈V ×V min(w 0 (e), w 1 (e)). Otherwise, we have to add edges in E 1 in order to get a connected graph (that is a spanning graph) and the "cost" of this addition has to be minimized.

Thus, suppose that the graph G 1 contains at least two maximal connected components. Let cc(G 1 ) = (cc 1 (G 1 ), . . . , cc t (G 1 )) be the t ≥ 1 maximal connected components of G 1 . We have to connect these different components minimizing the increasing of the difference threshold. Let E 0 be the set of candidate edges constructed as follows. For every i, j, 1 ≤ i < j ≤ t, let {v i , v j } be an edge such that for every v ∈ V (cc i (G 1 )) and for every v ∈ V (cc j (G 1 )), we have w 1 ({v i , v j }) -w 0 ({v i , v j }) ≤ w 1 ({v, v }) -w 0 ({v, v }). In other words, {v i , v j } is an edge that minimizes the marginal cost for connecting cc i (G 1 ) and cc j (G 1 ). We add {v i , v j } in E 0 . Thus, we have to add exactly t -1 edges of E 0 in order to get a connected graph and we aim at minimizing the cost of this addition. More precisely, we get our optimal core network by finding, for every i, j, 1 ≤ i < j ≤ t, an edge e such that w cc ({u i , u j }) = w 1 (e) -w 0 (e), that is an edge e of minimum marginal cost between the maximal connected component cc i (G 1 ) and the maximal connected component cc j (G 1 ). Let E 0 * be such a subset of t -1 edges. Thus, we get that

δ * = e∈E 1 w 1 (e) + e∈V ×V |e / ∈E 1 ,e∈E 0 * w 1 (e) + e∈V ×V |e / ∈E 1 ,e / ∈E 0 * w 0 (e). Ob- serve that e∈V ×V |e / ∈E 1 ,e∈E 0 * w 1 (e)
is exactly the cost of a mnimum spanning tree of the graph G cc defined before.

We are now able to prove Theorem 1.

Proof.[of

* = E 1 ∪ E 0 * .
Having developed the core subnetwork extraction guaranteeing a connected core network (Algorithm 1). We proceed to assess its performance. Unstable conns. 2.1% Fig. 3: Consistency analysis for extracted core graphs. We performed a Leave-N-Out procedure to quantify the consistency across methods at 4 different parameter levels. The results show similar graphs for both methods. However our method, in blue, has a smaller number of connections that are not present or absent across all experiments, i.e. unstable connections (marked in red).

To assess the performance of our method, we compared our novel approach with the currently used [START_REF] Gong | Mapping Anatomical Connectivity Patterns of Human Cerebral Cortex Using In Vivo Diffusion Tensor Imaging Tractography[END_REF]: first, we compared the stability of the obtained binary graph across randomly chosen subpopulations; second, we compared connectivity prediction performance.

For this comparisons, we used an homogeneous set from the HCP500 datatest [START_REF] Sotiropoulos | Advances in diffusion MRI acquisition and processing in the Human Connectome Project[END_REF]:all subjects aged 21-40 with complete dMRI protocol, which resulted in 309 subjects (112 male). We obtained the weighted connectivity matrices between the cortical regions defined by the Desikan atlas [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] as done by Bassett et al. [START_REF] Bassett | Conserved and variable architecture of human white matter connectivity[END_REF]. To verify the untresholded graph construction, we computed the average degree, number connections over number of possible connections, on each subject. Bassett et al. [START_REF] Bassett | Conserved and variable architecture of human white matter connectivity[END_REF] reported an average degree of 0.20 and we obtained 0.20 ± 0.01 (min: 0.17, max: 0.25) showing our preprocessing in agreement.

Consistency of the Extracted Graph

To quantify the consistency of the core graph extraction procedure we performed 500 Leave-N-Out experiments. At each experiment randomly sampled 100 subjects from the total and computed the core graphs with both techniques. We performed the extraction at 4 different levels of the parameter for each technique, choosing the parameters such that the density of the resulting graph connections is stable across methods. Also, we reported the number of unstable connections, selected as the connections that were not present or absent in all experiments. We show the results of this experiment in Fig. 3. In this figure we can observe that the resulting graphs are similar, while the number of unstable connections is larger for Gong et al. [START_REF] Gong | Mapping Anatomical Connectivity Patterns of Human Cerebral Cortex Using In Vivo Diffusion Tensor Imaging Tractography[END_REF] by an order of magnitude.

Predicting Gender-Specific Connectivity

Fitting Fig. 4: Performance of core network as feature selection for a linear model for gender specific connectivity. We evaluate model fit (left) and prediction (right), Gong et al. [5] in green, and ours, in blue. We show the histograms of both values from our nested Leave-1 3 -Out experiment. In both measures, our approach has more frequent lower values, showing a better performance.

To assess model fit and prediction we implemented a nested Leave- 1 3 -Out procedure. The outer loop performs model selection on 1 3 of the subjects randomly selected. First, it computes the core graph of a population, with our approach and that of Gong et al. [START_REF] Gong | Mapping Anatomical Connectivity Patterns of Human Cerebral Cortex Using In Vivo Diffusion Tensor Imaging Tractography[END_REF]. Then, it selects the features F that are more determinant of gender classification using the f-test feature selection procedure. The features are taken from the core graph adding the connectivity weights to each subject. The inner loop performs model fitting and prediction using the selected features F . First, we randomly take 1 3 of the remaining subjects and fits a linear model on F for predicting gender. Second, we predict the values of the features F from the gender column. The outer loop is performed 100 times and the inner loop 500 times per outer loop. This totals 50,000 experiments. Finally, for each experiment, we quantify the prediction performance of the linear model at each inner loop with the mean squared error (MSE) of the prediction and Akaike Information Criterion (AIC) for model fitting.

We show the experiment's results in Fig. 4. In these results we can see that our approach, in blue, performed better than Gong et al. [START_REF] Gong | Mapping Anatomical Connectivity Patterns of Human Cerebral Cortex Using In Vivo Diffusion Tensor Imaging Tractography[END_REF], in green as the number of cases with lower AIC and MSE is larger in our case.

Discussion and Conclusion

We present, for the first time, an algorithm to extract the core structural connectivity network of a subject population while guaranteeing connectedness. We start by formalizing the problem and showing that, although the problem is very hard (it is NP-complete), we produce a polynomial time exact algorithm to extract such network when its number of nodes is large. Finally, we show an example in which that our network constitutes a better feature selection step for statistical analyses of structural connectivity. For this, we performed a nested leave- 1 3 -out experiment on 300 hundred subjects. The results show that performing feature selection with our technique outperforms the most commonly used approach.

Figure 1 :

 1 Figure 1: Clustering Process: the tractograms are transformed into a vectorial space where we can manipulate and average them.

Figure 1 :

 1 Figure 1: Clustering Process: the tractograms are transformed into a vectorial space where we can manipulate and average them.

Figure 1 :

 1 Figure 1: Clustering Process: the tractograms are transformed into a vectorial space where we can manipulate and average them.

Figure 1 :

 1 Figure 1: Clustering Process: the tractograms are transformed into a vectorial space where we can manipulate and average them.

Figure 1 :

 1 Figure 1: Clustering Process: the tractograms are transformed into a vectorial space where we can manipulate and average them.

Fig. 1 :

 1 Fig. 1: Scheme of analyses involving the core structural connectivity matrix.

Fig. 2 :

 2 Fig. 2: Instance of the common sub-network problem. (a,b,c) Brain connectivity of different subjects, namely G 1 , G 2 and G 3 . (d) Extracted common sub-network G * that is optimal for n = 5 with λ = 1 2 : the difference threshold is 7 2 .

  Small example: Consider the instance depicted in Figure2with λ = 1 2 . Figure2(a), (b), and (

  and consider any real number λ ∈ [0, 1]. Then, Core-Sum-Alg (Algorthim 1) is an O(max(k, log(|V |)).|V | 2 )-time complexity exact algorithm for the core subnetwork problem when n = |V |. Algorithm 1 Core-Sum-Alg: Exact polynomial time complexity algorithm for the core sub-network problem when n = |V |. Require: SC graphs for each subject G1 = (V, E1), . . . , G k = (V, E k ), and λ ∈ [0, 1] Start Computing a Core Graph that can have Multiple Connected Components: 1: Construct G = (V, V × V ), the completely connected graph G 2: Compute w0(•), w1(•) across all subject graphs as in Eq. 1. 3: Compute the set of edges to be added to the core graph E 1 = {e|w1(e) ≤ w0(e)} and construct G1 = (V, E 1 ). Compute the Connected Components of G1 and Connect Them 4: Compute the set of maximal connected components cc(G1) = (cc1(G1) . . . cct(G1)) 5: Construct Gcc (Vcc, Vcc × Vcc) with Vcc = {u1, . . . , ut} 6: Compute wcc as in Eq. 4. 7: Compute the set of edges E 0 that correspond to argument minimum of Eq. 4. In other words, for every {ui, uj}, select the edge e connecting the maximal connected components cci(G1) and ccj(G1) such that w1(e) -w0(e) = wcc({ui, uj}). 8: Compute a minimum spanning tree Tcc of Gcc 9: Compute the set of edges E 0 * ⊆ E 0 that corresponds to the set of edges of the previous minimum spanning tree. 10: Construct G * = (V * , E * ) with V * := V and E * := E 1 ∪ E 0 * 11: return G * the connected Core Structural Connectivity Network

  Theorem1] Core-Sum-Alg (Algorithm 1) follows the proof of Lemma 1 and so solves the core sub-network problem when n = |V |. The construction of G (line 1) can be done in linear time in the number of edges, that is in O(|V | 2 )time. The time complexity of line 2 is O(k|V | 2 ). The construction of G 1 (line 3) can be done in O(|V | + |E 1 |)-time, O(|V | 2 )-time in the worst case. The computation of the maximal connected components of G 1 (line 4) can be done in linear time in the size of G 1 , that is in O(|V | + |E 1 |), O(|V | 2 )-time in the worst case. The construction of G cc (line 5) can be done in linear time in the size of G cc , that is, in the worst case, in O(|V | 2 )-time. The time complexity of line 6 and 7 is O(|V | 2 ). There is an O(m log(n))-time complexity algorihm for the problem of computing a minimum spanning tree of a graph composed of n nodes amd m edges (line 8). Thus, in our case, we get an O(log(|V |)|V | 2 )-time algorithm. The time complexity of line 9 is O(|V |). Finally, the construction of G * = (V * , E * ) (line 10) can be done in constant time because V * = V and E

Acknowledgments

This work has received funding from the European Research Council (ERC Advanced Grant agreement No 694665).