Proceedings of ITC/ISHW2007

Quantifying profile stiffness
B.Ph. VAN MILLIGEN1) , R. SÁNCHEZ2) , V. TRIBALDOS1) , and V.I. VARGAS1)
1)

2)

Asociación EURATOM-CIEMAT para Fusión, Avda. Complutense 22, 28040 Madrid, Spain
Fusion Energy Division, Oak Ridge National Laboratory, P.O. Box 2001, Oak Ridge, TN 37831-2001, USA
October 10, 2007

Profile stiffness is quantified using a simple technique. The approach is tested on a paradigmatic numerical
stiff transport model for one field (particles). The stiffness is found to exhibit radial structure and to depend on
collisionality, which might help explaining the observed lack of stiffness in stellarators, as compared to tokamaks.
The extension of the approach to heat transport requires some care. A proposal for a stiffness quantifier for heattransport is made, and it is tested on data from the TJ-II stellarator.
Keywords: profile stiffness, anomalous transport, critical gradient, tokamak-stellarator comparison
DOI: 10.1585/pfr.1.001
Profile stiffness (also known as profile consistency orresilience) is the striking phenomenon that temperature orpressure profiles tend to adopt the same shape, regardlessof the applied drive, at least in a certain parameter range.
The phenomenon is well-established for tokamaks, but elusive in stellarators [1]. This is slightly enigmatic, since (a)power degradation is a universal phenomenon in stellarators, with a similar power dependence as in tokamaks [2]and (b) it is believed that threshold-triggered instabilities(leading to enhanced transport) should be operative bothin tokamaks and stellarators. This leads to the expectationthat profile stiffness should also be present in stellarators,
if to a less degree (and less obviously) than in tokamaks.
The detection of profile stiffness based on the directcomparison of profiles does not allow a quantification of
the degree of stiffness, while the full 1-D modelling oftransport requires making assumptions not related to thestiffness issue. Therefore, a stiffness quantifier is needed toresolve this issue. In the present work we apply a standard
quantifier for profile stiffness to a paradigmatic stiff particle transport model. We then discuss the possible (nonstandard) extension of the method to heat transport andpresent first results for the TJ-II stellarator.
Profile stiffness can be understood as the sub-linear response of profile amplitudes to a (small) change in drive.
Pure diffusive transport models with fixed parameters produce a proportional response of profiles to changes in fuelling or heating, since the diffusion equation is linear in
the profile amplitude and the source strength. Thus, thesearch for profile stiffness is closely related to the study ofthe dependence of transport coefficients on (gradients of)
the transported quantity, since such a dependence would
break the linearity of the diffusion equation.
Such studies have been undertaken before [3], in theframework of the analysis of perturbative transport. In thecited reference, a distinction was made between the steadyauthor’s e-mail: boudewijn.vanmilligen@ciemat.es

state (or ”power balance”) transport coefficient Dpb =−Γ/∇n (here, Γ is the particle flux and ∇n the density gradient), and the perturbation response value Dinc = −∂Γ/∂∇n,
dubbed the ”incremental” transport coefficient. If Dinc >Dpb , the profiles will respond sub-linearly to changes inthe source term, thus producing stiffness. Accordingly, a”stiffness factor” can been defined (by analogy to [4]):Dinc ∇n ∂Γ
=
,
(1)
Γ ∂∇n
Dpb
so that C > 1 would indicate some degree of stiffness.
The evaluation of Dinc requires a (small) perturbationof the source term and profiles, either spontaneous or induced externally. However, the relevant variables of systems near a critical steady state tend to fluctuate spontaneously around a mean value. This property can be exploited to obtain another, equivalent estimate of the stiffness that does not require perturbing the system. Interpreting the mean amplitude of the fluctuations around thesteady state values (i.e. their standard deviation) as thesmall change symbolised by ∂ symbol in Eq. (1) [5]:C=

Dfluct = RMS (Γ)/RMS (∇n)
(2)
D
E1/2
where RMS ( f ) = ( f − h f i)2, and the angular brackets refer to a time average. In steady state, and assumingthat the system response to perturbations is linear to first
approximation, we expect Dfluct ' Dinc .
In the following, we study the stiffness parameter C
using a simplified transport model, considered paradigmatic for transport controlled by a critical gradient. The
model is described in considerable detail elsewhere [6].
The simplified model is Markovian in nature and the time
evolution of the single field n(x, t), which may be interpreted as a (particle) density, obeys, in one dimension, aGeneralized Master Equation:∂n(x, t)= S (x, t) +∂t
Z 1n(x, t)1
dx0 p(x − x0 ; x0 , t)n(x0 , t) −
,
(3)
τD 0
τD

Proceedings of ITC/ISHW2007

When the local gradient is below the critical value (subcritical), transport is governed by the p0 step distribution,
and when it is above (super-critical), it is governed by the
p1 step distribution. Here, p0 and p1 are fixed and symmetric stable probability distributions (of the Lévy type,of which the Gaussian is a special case). Transport at any
given location x will therefore be sub- or super-critical asa function of the local value of the density gradient. This
introduces a mechanism for self-regulation into the model.
In this work, p0 is always chosen to be a Gaussian(with width σ0 ), while p1 can either be a Gaussian (withwidth σ1 ) or a Cauchy distribution p1 (x − x0 , x0 , t) =σ1 /π(σ21 + (x − x0 )2 ). While a Gaussian distribution models’normal’ diffusive transport, a Cauchy distribution (with a’long tail’) is used to model processes with long-range correlations, typically called ’avalanches’ or ’streamers’ in theplasma transport context, and representative of turbulent or’anomalous’ transport.
To compute Dinc , we will be comparing steady state
profiles at slightly different values of the amplitude of thesource S . For simplicity,
the steady state flux Γ is comRx
puted as Γ(x) = 0.5 dx0 S (x0 ). The lower limit of the integral corresponds to the system centre at x = 0.5. This
calculation is sufficiently accurate for the purpose of evaluating Dinc . However, for the calculation of Dfluct we will
use a different estimate of the flux that includes fluctuating
contributions (see below).
The numerical calculations are performed in the domain 0 ≤ x ≤ 1, on a grid with either N = 2000 (highresolution) or N = 200 (low resolution) grid points, using
standard integration techniques for stiff differential equations to advance Eq. (3) in time. In all cases, the source
S (x, t) = S 0 is taken constant. We first performed a scanof the source rate at high resolution (13 cases), and after<n>Gauss-Gauss

1000

500

12000
Gauss-Cauchy

1000
800
600
400
200
0

Thus, the model is closely related to standard transportmodels in common use.
The step distribution p is chosen as follows to produce
the required critical gradient mechanism:(
p0 : |∇n| < [∇n]crit (sub − critical)p=(5)
p1 : |∇n| ≥ [∇n]crit (super − critical)<n> (S0=0.01)<n> (S0=0.02)<n> (S0=0.05)<n> (S0=0.1)<n> (S0=0.2)<n> (S0=0.5)

1500

<n>

The domain of the system is 0 ≤ x ≤ 1, implying a normalisation of the spatial scales of the system to the systemsize. τD is a waiting time and specifies the mean time aparticle remains at a given location before taking a step.
We set τD = 1, implying a normalisation of the time scalesof the system to the mean waiting time. S (x, t) is an external particle source, and compensates edge losses due tothe absorbing boundary conditions imposed at x = 0, 1.
The function p is a ”particle steph distribution”. Whenp isi√Gaussian, p(x − x0 , x0 , t) = exp −(x − x0 )2 /4σ2 /2σ π,
standard diffusion is recovered in the limit of small step
sizes σ (and assuming a smooth density profile [7]):
"
#
∂n
∂2 σ2
= 2
n + S.
(4)
∂t
∂x τD

0

0.5

x

1

Fig. 1 Profiles vs. S 0 : Gauss-Gauss (top) and Gauss-Cauchy(bottom).

checking that the results at lower resolution were equivalent, we performed a bi-dimensional parameter scan ofboth the source rate and the sub-critical diffusion coefficient, the latter being proportional to σ20 (130 cases).
Source scan: More details on the high-resolution calculations discussed here can be found in Ref. [8]. We set
σ0 = 0.002, while p1 is Gaussian with σ1 = 0.008 for thecases labelled Gauss-Gauss or ”GG” (both transport channels are Gaussian), or p1 is Cauchy with σ1 = 0.004 for thecases labelled Gauss-Cauchy or ”GC”. The critical gradient is chosen [∇n]crit = 2000.
For both series (GG and GC), a scan of
the source rate was performed, choosing S 0 ∈{0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 10}. The parameter Dinc
was computed by comparing the profiles corresponding to
two subsequent values of S 0 .
Steady state profiles are shown in Fig. 1. The GG and
GC profiles are very similar, except for the highest fuellingrate (S 0 = 0.5): whereas the GC profile remains criticalacross the system, the GG profile ”bulges”, i.e. becomessuper-critical. This difference can be ascribed to the larger
transport capacity of the super-critical transport channel in
the GC case.
In order to quantify the stiffness, we computed thestiffness factor C. Fig. 2 shows the radial profiles of thestiffness factor C for the GC case. The spikes in the figure
occur when the local values of the gradient ∇n(x) corresponding to the two subsequent values of S 0 being analysed are equal. Such points should be ignored, since theirstatistical error is large, and thus we will base the analysis
on the global traces while ignoring the spikes. The figure
shows that the critical region (where C  1) grows fromthe edge inward as the source is increased. Next, a supercritical region (with C ' 1) starts to grow from the edge inward. The super-critical state covers almost the whole sys-

Proceedings of ITC/ISHW2007
104

C (S0=0.01/0.02)

Gauss-Gauss

Gauss-Cauchy

C (S0=0.02/0.05)

106

<C>

C (S0=0.1/0.2)

105

C (S0=0.2/0.5)

4

10
C

C (S0=0.2)
C (S0=0.6)
C (S0=2.0)

1000

C (S0=0.05/0.1)

100
10

1000
100

1

10

0.1

1

Gauss-Cauchy
1000

0.75

x

1

<C>

0.1
0.5

100
10

Fig. 2 Case Gauss-Cauchy: profiles of C for different combinations of the fuelling S 0 .

1
0.1

Super-critical (GG)

104

10
1
0.1

S0 (lower)

0.1

1000
Gauss-Gauss
100

Gauss-Gauss
Gauss-Cauchy
0.01

σ0

Fig. 4 hCi vs. σ0 for various values of S 0 , cases Gauss-Gauss(top) and Gauss-Cauchy (bottom).

<C>

100

Sub-critical

<C>

1000

0.01

Cfluct (S0=0.2)
Cfluct (S0=0.6)
Cfluct (S0=2.0)

10

0.1
1
1000

Fig. 3 Source scan, hCi vs. S 0 .

Gauss-Cauchy

tem in the GG case at highest fuelling (not shown), while
it only affects a narrow boundary layer for the highest fuelling GC case.
Fig. 3 shows hCi, the radial average of C. The points
in the figure are labelled by the lower of the two S 0 values
used. In comparison to the GG cases, the GC cases do
not only yield (slightly) larger values of C, but the rangeof values of S 0 where C  1 exceeds the corresponding
range for the GG cases.
Bi-dimensional parameter scan (source rate anddiffusivity). Here we scan the source amplitude S 0 and
σ0 . The latter can be interpreted as a scan of the subcritical diffusivity (or ”collisionality”). As in Ref. [6], we
set τD = 1, S (x) = S 0 , and [∂n/∂x]crit = 50. To compute
C, the same runs were repeated with S (x) = 1.1 · S 0 .
Fig. 4 shows hCi for the GG cases with σ1 = 0.08,
and for GC with σ1 = 0.04. For GG, the largest possible
value of σ0 is σ1 . For GC, no such limit exists on σ0 . It isobserved that hCi is a sensitive diagnostic for criticality.
In all series studied, hCi is seen to increase graduallywith increasing σ0 , reach a maximum value and then drop
somewhat abruptly and make a sharp transition to its subcritical expectation value (C = 1) at a precise value of σ0 .
For the GC cases, the point where the system transits froma fully sub-critical state to a critical state has been computed in previous work [6]. This critical power threshold√is given by: S c = 2σ20 /τD [∇n]crit , or σ0 = S c /10 withour choice of parameters. This matches the results exactly.
Stiffness from fluctuations. In the case of our numer-

<C>

100
10
1
0.1

Fig. 5

0.01

σ0

0.1

D
E
C fluct vs. σ0 and S 0 , cases Gauss-Gauss (top) and
Gauss-Cauchy (bottom).

ical model, it is straight-forward to compute the fluctuating gradient, while the flux can be evaluated from particle
conservation (∂n/∂t = −∂Γ/∂x + S ) in combination with
Eq. (3). Thus:
Z 1n(x0 , t)Γ(x, t) =dx0 K(x − x0 , x0 , t),
(6)
τD (x0 )0where K(x − x0 , x0 , t) = Θ(x − x0 ) − P(x R− x0 , x0 , t), Θ(x) is∆the Heaviside function and P(∆, x0 , t) = −∞ dxp(x, x0 , t) isthe cumulative step probability distribution.
D
E
Fig. 5 shows the calculation of C fluct . Compare theseresults
D toEFig. 4. Although the maximum numerical valueof C fluct is somewhat lower than that of hCi, the global
trend is the same. The deviation between C and C fluct atpoints with large stiffness is to be expected, as the system
response will be strongly non-linear at such points. The
calculation of C fluct is not possible when the system is locally static, which explains why these figures have less data
points than Fig. 4.

Proceedings of ITC/ISHW2007

C LTe =

1 ∂(χ∇ ln T ).
χ ∂(∇ ln T )(7)This unorthodox proposal is designed to detect the dependence of the heat diffusivity, χ, on the expected criticalparameter for the ETG instability, ∇T/T . Note that many
alternative definitions are possible.
Application to the TJ-II Stellarator. Here we reporton the first attempt to estimate the stiffness of the temperature profile in the stellarator TJ-II. Profiles at TJ-II are
obtained using the single-pulse high-resolution Thomson
Scattering diagnostic [11], yielding around 200 data points
for the electron temperature T and density n along a chord
spanning most of the plasma cross section, with a spatial
resolution of 2.25 mm.
The goal of the present analysis is to determine theglobal transport response to a change in heating. Therefore, we fit the profiles to simple functional forms, thus
ignoring any detailed radial structure. This improves therobustness of the calculation of radial derivatives needed tocompute C LTe . The temperature profile is fit to the sum of
two Gaussians, while the density profile is fit to a Gaussian
multiplied by a second-order polynomial in ρ2 (for sym√
metry). ρ = ψ is a radial coordinate, where ψ is thenormalised poloidal magnetic flux, obtained from the theoretical calculation of the magnetic flux surfaces in vacuum. Finite pressure effects can safely be ignored. The
discharges studied here are those reported in Ref. [12].
The error in the profile reconstruction, evaluated using the Jacobian of the fit matrix, is of the order of 10%,lower in the centre but increasing towards the edge. While
the temperature profile reconstruction is reliable (i.e. withan error less than 10%) out to about ρ = 0.7, the density
profile reconstruction is reliable only out to about ρ = 0.4.
The calculation of C LTe is not very sensitive to the detailsof the density profile.
To compute χ and C LTe , an estimate of the heat flux q
is required. The heating source is assumed to have a Gaussian deposition profile, centred at ρ = 0, with a fixed widthof ∆ρ = 0.2. The heating efficiency is estimated to be

9 shot combinations 200/400 kW
1000
100
<C>

Heat transport. The preceding analysis was simplified by the fact that the thermodynamic force and the critical parameter were both equal to ∇n. For heat transport,these two quantities do not coincide and the definition ofa stiffness quantifier is less obvious. The traditional proposal, C = χinc /χpb = −∂q/∂(n∇T )/χ [4], where χ is theheat diffusivity, will respond only weakly when the system
criticality is not determined by n∇T ' const. In general, a
stiffness quantifier C that responds sharply to a given critical condition should be proportional to the inverse of thechange in that critical condition, as might indeed be deduced from the modelling efforts in, e.g., [9]. Since we expect Electron Temperature Gradient modes to play a role inthe stiffness (if any) of the temperature profile, we believethat the critical parameter must be ∇T/T [10], so that we
define the stiffness of the temperature profile by

10
1
0.1

-1

-0.5

0

ρ

0.5

1

Fig. 6 Stiffness estimate for TJ-II.

60%. Radiation and other losses are ignored. The heat flux
is obtained by integrating the net deposited power. In any
case, the stiffness factor C LTe is not very sensitive to thedetails of this calculation. Fig. 6 shows the mean stiffness
factor, averaged over 9 equivalent discharge combinations
with similar densities and different heating levels. Ignoring the Dspikes,
E one observes that a certain profile stiffness
exists ( C LTe > 1) in the radial range 0.15 ≤ ρ ≤ 0.55,
roughly coincident with the T gradient region.
Discussion. The quantification of profile stiffness isdirectly related to the detection of the dependence of thetransport coefficient on the profile gradient. In accordancewith this idea and with literature, we make use of a stiffness
quantifier C, and show that it provides a useful quantification of stiffness in a paradigmatic transport model. It appears that stiffness has a radial structure and a dependence
on system parameters (such as the source or drive, and thecollisionality), which could possibly shed some light on
the observed differences between tokamaks and stellarators.
The extension of these results to heat transport requires some care. We have suggested a definition and applied it to TJ-II data.
[1] U. Stroth, Plasma Phys. Control. Fusion 40, 9 (1998).
[2] B. Carreras, IEEE Trans. Plasma Sc. 25, 1281 (1997).
[3] N.J. Lopes Cardozo, Plasma Phys. Control. Fusion 37, 799(1995).
[4] F. Ryter, F. Leuterer, G. Pereverzev, et al., Phys. Rev. Lett.
86, 2325 (2001).
[5] B.A. Carreras, V.E. Lynch, P.H. Diamond and M. Medvedev,
Phys. Plasmas 5, 1206 (1998).
[6] B.Ph. van Milligen, R. Sánchez, and B.A. Carreras, Phys.
Plasmas 11, 2272 (2004) and ibid., 11, 3787 (2004).
[7] R. Sánchez, B.Ph. van Milligen, and B.A. Carreras, Phys.
Plasmas 12. 056105 (2005).
[8] B.Ph. van Milligen, B.A. Carreras, V.E. Lynch, and R.
Sánchez, Nucl. Fusion 47, 189 (2007).
[9] F. Imbeaux, F. Ryter, and X. Garbet, Plasma Phys. Control.
Fusion 43, 1503 (2001).
[10] X. Garbet, N. Dubuit, E. Asp, et al., Phys. Plasmas 12,082511 (2005).
[11] C.J. Barth, F.J. Pijper, H.J. van der Meiden, et el., Rev. Sci.
Instrum. 70, 763 (1999).
[12] V.I. Vargas, D. López-Bruna, J. Herranz, F. Castejón, and
TJ-II Team, Nucl. Fusion, 47, 1367 (2007).