Shafranov Shift in Tokamak Plasmas

• Shafranov shift is the outward displacement of magnetic flux surfaces in toroidal plasmas caused by finite plasma pressure and modified force balance.
• It is derived from the Grad–Shafranov equation and quantified by parameters like Δ and A′, playing a key role in equilibrium reconstruction and heat-flux mapping.
• The shift impacts plasma stability, intrinsic rotation, and turbulent transport, with implications for both tokamak and stellarator designs.

The Shafranov shift is the outward displacement of magnetic flux surfaces, or equivalently of the magnetic axis in many tokamak formulations, that arises in toroidal magnetohydrodynamic equilibrium when finite plasma pressure and the associated poloidal magnetic field modify the force balance. In large-aspect-ratio tokamaks it is commonly written as a radial shift $\Delta(r)$ of each flux surface or as $\Delta = R_{\rm axis}-R_0$, while recent stellarator work introduces an averaged weighted shift $\delta S$ adapted to non-axisymmetric geometry (Mitteau et al., 2021, Ciro et al., 2014, Helander et al., 21 May 2026). Although often only of order centimeters in conventional devices, the shift enters equilibrium reconstruction, edge heat-flux mapping, ballooning stability, intrinsic rotation, Alfvén-eigenmode dynamics, and the design of plasma-facing components (Mitteau et al., 2021, Connor et al., 2016, Rofman et al., 12 Jan 2026).

1. Physical origin and definitions

In the idealized large-aspect-ratio tokamak with purely circular, concentric magnetic surfaces, each flux surface is centered on the major-radius axis $R=R_0$. Real equilibria depart from this limit because the plasma pressure gradient and the poloidal magnetic field produce a pressure-driven hoop force that pushes the plasma column outward, while the return poloidal field reacts to maintain force balance. The net result is an outward displacement of magnetic surfaces that increases with minor radius and with plasma pressure or poloidal beta (Mitteau et al., 2021).

In axisymmetric notation, the shift is often defined through the magnetic-axis position. One common normalized form is

$\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$

where $R_0$ is the geometric center of the poloidal cross-section and the magnetic axis lies at $R=R_0(1+\Delta)$ (Ciro et al., 2014). In equilibrium and stability studies with nearly circular surfaces, the flux-surface geometry is often written as

$$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$

with $\Delta(0)=0$ and $A'(r)\equiv d\Delta/dr$ used as the local Shafranov-shift parameter (Connor et al., 2016).

The terminology is not completely uniform. In tokamak heat-load calculations, $\Delta = R_{\rm axis}-R_0$0 denotes the outward displacement of each surface and directly enters the last closed flux surface (LCFS), scrape-off-layer (SOL) distance, and incidence angle of field lines on plasma-facing components (Mitteau et al., 2021). In some analytic equilibrium papers, the emphasis is instead on the outward displacement of the magnetic axis obtained from the condition $\Delta = R_{\rm axis}-R_0$1 (Palha et al., 2015, Li et al., 2019). In stellarators, where a single axis displacement is less informative, an “average” shift is defined through a volume functional of the perturbed flux (Helander et al., 21 May 2026). These usages describe related but not identical geometric diagnostics.

Several immediate consequences follow from the outward displacement. The center of the LCFS is shifted to $\Delta = R_{\rm axis}-R_0$2; the SOL width seen by a component at fixed $\Delta = R_{\rm axis}-R_0$3 is modified; and the local connection length and the e-folding length $\Delta = R_{\rm axis}-R_0$4 acquire a weak poloidal asymmetry that is shorter on the outboard side in the large-aspect-ratio picture (Mitteau et al., 2021).

2. Grad–Shafranov formulation and large-aspect-ratio expressions

The axisymmetric equilibrium is governed by the Grad–Shafranov equation for the poloidal flux $\Delta = R_{\rm axis}-R_0$5,

$\Delta = R_{\rm axis}-R_0$6

or in equivalent forms used across the equilibrium literature (Mitteau et al., 2021, Palha et al., 2015, Ciro et al., 2014). In this equation the right-hand side contains the pressure-gradient term and the toroidal-field-function term; the Shafranov shift is therefore a direct equilibrium consequence of finite pressure and current profiles.

For large aspect ratio $\Delta = R_{\rm axis}-R_0$7, circular flux surfaces, and low $\Delta = R_{\rm axis}-R_0$8, an analytic expansion yields an outward shift of the flux surface labeled by minor radius $\Delta = R_{\rm axis}-R_0$9,

$\delta S$0

For a parabolic pressure profile $\delta S$1, this becomes

$\delta S$2

A frequently used edge estimate is

$\delta S$3

with $\delta S$4 the poloidal beta and $\delta S$5 the internal inductance (Mitteau et al., 2021).

A related representation emphasizes the derivative of the shift rather than the shift itself. For nearly circular equilibria,

$\delta S$6

which makes explicit that the local shift parameter is proportional to a pressure-weighted radial moment and scales as $\delta S$7 at the edge (Connor et al., 2016).

In semi-analytical linearized Grad–Shafranov models, the shift is tied to the eigenvalue of the homogeneous Helmholtz-type problem through the on-axis safety factor, toroidal-field fraction, and beta. In that framework, matching the safety factor at the magnetic axis establishes an algebraic relation between the eigenvalue $\delta S$8, the equilibrium parameters, and the prescribed shift $\delta S$9 (Ciro et al., 2014). This formulation is especially useful when scanning equilibria with specified elongation, triangularity, aspect ratio, and magnetic shear.

3. Extraction from analytic and numerical equilibria

In analytic equilibrium families, the Shafranov shift is obtained by locating the magnetic axis from $R=R_0$0. For the classical Soloviev model used in fixed-boundary studies, the constants in the polynomial flux solution are chosen so that $R=R_0$1 on a prescribed tokamak-like boundary, after which the axis position follows from the root of $R=R_0$2 (Palha et al., 2015). In the examples reported for that construction, an ITER-like case with $R=R_0$3, $R=R_0$4, and $R=R_0$5 gives $R=R_0$6, while an NSTX-like case with $R=R_0$7, $R=R_0$8, and $R=R_0$9 gives $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$0 (Palha et al., 2015).

A related semi-analytical Grad–Shafranov solver imposes $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$1 and $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$2, expands the solution as a particular piece plus a homogeneous Helmholtz-type contribution, and fits the coefficients to a prescribed D-shape with possible X-points. In the normal-shear example with $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$3, $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$4, $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$5, $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$6, $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$7, and $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$8, the fit converges at $\Delta \equiv \frac{R_{\rm axis}-R_0}{R_0},$9, $R_0$0, $R_0$1, and finds $R_0$2. In a reversed-shear case with $R_0$3, $R_0$4, $R_0$5, $R_0$6, $R_0$7, and $R_0$8, the fit yields $R_0$9, $R=R_0(1+\Delta)$0, $R=R_0(1+\Delta)$1, and $R=R_0(1+\Delta)$2 (Ciro et al., 2014).

High-order numerical solvers treat the shift as a derived quantity whose accuracy inherits the accuracy of the reconstructed flux. The mimetic spectral-element Grad–Shafranov solver reported in (Palha et al., 2015) enforces discrete topological laws exactly, achieves $R=R_0(1+\Delta)$3-convergence of order $R=R_0(1+\Delta)$4 and spectral convergence in polynomial degree $R=R_0(1+\Delta)$5, and recovers $R=R_0(1+\Delta)$6 from root-finding on $R=R_0(1+\Delta)$7. For moderate meshes and $R=R_0(1+\Delta)$8, the resulting shift is recovered to within $R=R_0(1+\Delta)$9–$$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$0. In high-shift cases with $$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$1–$$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$2 on a coarse $$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$3 mesh, the residual error peaks at order $$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$4 for $$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$5 and falls below $$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$6 for $$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$7 (Palha et al., 2015).

Spectral-element equilibrium solvers with toroidal rotation preserve the same logic: the magnetic axis is still obtained from $$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$8 at $$R(r,\theta)=R_0-r\cos\theta+\Delta(r)+O(\epsilon^2),$$9, but the equilibrium equation includes centrifugal terms. In the rotating Solov’ev family, this yields an explicit modification of the static shift and provides an analytic benchmark for extended equilibrium solvers (Li et al., 2019).

4. Incorporation into heat-flux mapping and plasma-facing-component design

In heat-flux deposition models, the Shafranov shift is not a minor geometric refinement; it enters the field-line geometry, the SOL mapping, and the local incidence angle on plasma-facing components. In the Tokaflu implementation within Castem 2000, the LCFS in the meridional plane is written

$\Delta(0)=0$0

while a general surface of minor radius $\Delta(0)=0$1 is represented by $\Delta(0)=0$2 (Mitteau et al., 2021).

The local magnetic field is then evaluated from large-aspect-ratio Shafranov expressions for the poloidal and toroidal components, combined into the unit vector $\Delta(0)=0$3 along the field. The SOL coordinate used in the heat-flux decay is

$\Delta(0)=0$4

so that the local heat-flux density is

$\Delta(0)=0$5

with $\Delta(0)=0$6 the surface normal. Because the shift changes both the position of the flux surfaces and the field-line direction, it changes the exponential decay factor and the cosine incidence factor simultaneously (Mitteau et al., 2021).

For typical Tore Supra advanced-scenario parameters $\Delta(0)=0$7, $\Delta(0)=0$8, $\Delta(0)=0$9, $A'(r)\equiv d\Delta/dr$0, $A'(r)\equiv d\Delta/dr$1–$A'(r)\equiv d\Delta/dr$2, and $A'(r)\equiv d\Delta/dr$3, the quoted edge shift is $A'(r)\equiv d\Delta/dr$4–$A'(r)\equiv d\Delta/dr$5, reaching up to $A'(r)\equiv d\Delta/dr$6 in the most extreme cases. At $A'(r)\equiv d\Delta/dr$7, $A'(r)\equiv d\Delta/dr$8–$A'(r)\equiv d\Delta/dr$9; at $\Delta = R_{\rm axis}-R_0$00, $\Delta = R_{\rm axis}-R_0$01–$\Delta = R_{\rm axis}-R_0$02 (Mitteau et al., 2021).

These values alter predicted heat-flux patterns. Outboard flux surfaces are shifted outward, which decreases the local connection length and slightly steepens the SOL decay, producing a narrower and more intense heat-flux footprint on outer limiters; the high-field side shows the opposite tendency. In Tokaflu calculations for the Tore Supra toroidal pumped limiter, including $\Delta = R_{\rm axis}-R_0$03 raised the peak local parallel heat flux by approximately $\Delta = R_{\rm axis}-R_0$04–$\Delta = R_{\rm axis}-R_0$05, shifted the hot-spot location by approximately $\Delta = R_{\rm axis}-R_0$06–$\Delta = R_{\rm axis}-R_0$07 toroidally and poloidally, and modified incidence angles by up to $\Delta = R_{\rm axis}-R_0$08 (Mitteau et al., 2021).

The design implications reported for the CIEL project are correspondingly direct. Shafranov-shift-aware heat-flux maps were used to select the optimal leading-edge radius on the toroidal pumped limiter so as to limit peak $\Delta = R_{\rm axis}-R_0$09 to approximately $\Delta = R_{\rm axis}-R_0$10, to define tile orientations such that the worst-case incidence angle never exceeded $\Delta = R_{\rm axis}-R_0$11, and to place cooling channels beneath predicted hot spots (Mitteau et al., 2021).

5. Consequences for stability, intrinsic rotation, and turbulent transport

In ideal high-$\Delta = R_{\rm axis}-R_0$12 ballooning theory, the Shafranov shift enters through the local shift parameter $\Delta = R_{\rm axis}-R_0$13 and modifies the coefficients of the low-shear ballooning equation. With finite aspect ratio and ellipticity included, increasing $\Delta = R_{\rm axis}-R_0$14 raises the critical pressure-gradient parameter $\Delta = R_{\rm axis}-R_0$15, so that larger shift stabilizes ballooning modes. In the reported self-consistent calculations, increasing $\Delta = R_{\rm axis}-R_0$16 from $\Delta = R_{\rm axis}-R_0$17 to $\Delta = R_{\rm axis}-R_0$18 shifts the marginal-stability curve to the right by approximately $\Delta = R_{\rm axis}-R_0$19–$\Delta = R_{\rm axis}-R_0$20 for the circular case, while the shaped case is even more stable (Connor et al., 2016). The stated interpretation is that outward displacement reduces the effective curvature drive seen by the mode.

In up-down asymmetric tokamaks, the shift also modifies intrinsic momentum transport. A large-aspect-ratio expansion for tilted elliptical boundaries shows that the Shafranov shift becomes up-down asymmetric and depends strongly on the tilt angle of the flux surfaces, yet is insensitive to the detailed shape of the current and pressure profiles when the geometry, total plasma current, and average pressure gradient are kept fixed (Ball et al., 2016). This establishes a useful distinction between global control parameters and detailed profile shaping.

Nonlinear electrostatic GS2 simulations of these equilibria show that the shift alone can enhance intrinsic momentum flux by approximately $\Delta = R_{\rm axis}-R_0$21, but self-consistent inclusion of $\Delta = R_{\rm axis}-R_0$22 strongly reduces the normalized momentum-to-heat-flux ratio $\Delta = R_{\rm axis}-R_0$23 and largely cancels the shift-driven enhancement. The combined effect broadens the rotation profile while leaving the on-axis rotation roughly unchanged; the on-axis Alfvén Mach number remains of order $\Delta = R_{\rm axis}-R_0$24 for ITER-like parameters (Ball et al., 2016). The same study also shows that a pure shift in an up-down symmetric circular cross-section drives very little $\Delta = R_{\rm axis}-R_0$25, so the shift must act together with flux-surface shaping to break symmetry.

Global gyrokinetic simulations with ORB5 extend this role of the Shafranov shift to Alfvén eigenmodes and microinstabilities in finite-$\Delta = R_{\rm axis}-R_0$26 plasmas with energetic particles. In one study, self-consistent finite-$\Delta = R_{\rm axis}-R_0$27 equilibria with EP pressure reduce the $\Delta = R_{\rm axis}-R_0$28 TAE growth rate from approximately $\Delta = R_{\rm axis}-R_0$29 to approximately $\Delta = R_{\rm axis}-R_0$30 in a standard-profile case at $\Delta = R_{\rm axis}-R_0$31, and from approximately $\Delta = R_{\rm axis}-R_0$32 to approximately $\Delta = R_{\rm axis}-R_0$33 at $\Delta = R_{\rm axis}-R_0$34; in the same framework, a peaked-profile KBM unstable at $\Delta = R_{\rm axis}-R_0$35 in $\Delta = R_{\rm axis}-R_0$36 geometry is fully suppressed by the consistent MHD shift, and an ITG mode is reduced from approximately $\Delta = R_{\rm axis}-R_0$37 to approximately $\Delta = R_{\rm axis}-R_0$38 (Rofman et al., 12 Jan 2026). A separate ORB5 study finds that consistent inclusion of EP pressure in the MHD equilibrium can reduce the $\Delta = R_{\rm axis}-R_0$39 TAE growth rate by approximately $\Delta = R_{\rm axis}-R_0$40–$\Delta = R_{\rm axis}-R_0$41, with the stabilization strongest at low toroidal mode number (Rofman et al., 16 May 2026). Together these results identify the shift as both an equilibrium effect and a transport-relevant ingredient in self-consistent finite-$\Delta = R_{\rm axis}-R_0$42 modeling.

6. Rotation, energetic-particle pressure, and stellarator generalizations

Rigid toroidal rotation modifies the Shafranov shift through centrifugal terms in the Grad–Shafranov equation. For the rotating Solov’ev equilibrium with constant $\Delta = R_{\rm axis}-R_0$43, constant $\Delta = R_{\rm axis}-R_0$44, and the same linear pressure and toroidal-field ansatz as in the static case, the shift is

$\Delta = R_{\rm axis}-R_0$45

with

$\Delta = R_{\rm axis}-R_0$46

Hence $\Delta = R_{\rm axis}-R_0$47 for any finite $\Delta = R_{\rm axis}-R_0$48, and in the weak-rotation limit $\Delta = R_{\rm axis}-R_0$49 (Li et al., 2019). This places toroidal flow in the same equilibrium category as pressure: both increase the outward displacement.

Energetic-particle pressure alters the shift by contributing directly to the equilibrium pressure profile. In CHEASE-generated finite-$\Delta = R_{\rm axis}-R_0$50 equilibria used by ORB5, the shift at $\Delta = R_{\rm axis}-R_0$51 in standard-profile cases increases from approximately $\Delta = R_{\rm axis}-R_0$52 at $\Delta = R_{\rm axis}-R_0$53 to approximately $\Delta = R_{\rm axis}-R_0$54 at $\Delta = R_{\rm axis}-R_0$55 and approximately $\Delta = R_{\rm axis}-R_0$56 at $\Delta = R_{\rm axis}-R_0$57; in peaked-profile cases the corresponding values are approximately $\Delta = R_{\rm axis}-R_0$58, $\Delta = R_{\rm axis}-R_0$59, and $\Delta = R_{\rm axis}-R_0$60 (Rofman et al., 12 Jan 2026). In a second ORB5 study with circular flux surfaces, the same radial location yields $\Delta = R_{\rm axis}-R_0$61, $\Delta = R_{\rm axis}-R_0$62, and $\Delta = R_{\rm axis}-R_0$63, emphasizing the rough proportionality $\Delta = R_{\rm axis}-R_0$64 (Rofman et al., 16 May 2026).

The concept also extends beyond axisymmetry. For stellarators, an average weighted shift is defined by

$\Delta = R_{\rm axis}-R_0$65

so that the first-order shift is

$\Delta = R_{\rm axis}-R_0$66

Using reduced MHD and an auxiliary Poisson problem for $\Delta = R_{\rm axis}-R_0$67, one obtains

$\Delta = R_{\rm axis}-R_0$68

which shows that, to lowest order, the Shafranov shift is carried by the parallel Pfirsch–Schlüter current (Helander et al., 21 May 2026). Dimensional estimates then give $\Delta = R_{\rm axis}-R_0$69 up to geometry factors.

For special optimized stellarator classes, more specific scalings are reported. In quasisymmetric $\Delta = R_{\rm axis}-R_0$70-periodic fields,

$\Delta = R_{\rm axis}-R_0$71

while in quasi-isodynamic fields

$\Delta = R_{\rm axis}-R_0$72

Fixed-boundary VMEC equilibria at $\Delta = R_{\rm axis}-R_0$73 show the largest weighted shift in the precise QA configuration, up to $\Delta = R_{\rm axis}-R_0$74; Wendelstein 7-X is somewhat smaller at $\Delta = R_{\rm axis}-R_0$75; the QI design SQuID is about $\Delta = R_{\rm axis}-R_0$76; and the precise QH configuration is smallest at $\Delta = R_{\rm axis}-R_0$77 (Helander et al., 21 May 2026). In that sense, quasi-helical and quasi-isodynamic stellarators with large field-period number are reported to be particularly robust to pressure-driven expansion.

A recurring theme across these extensions is that the Shafranov shift is not a single universal scalar but a geometry-dependent equilibrium response. In static tokamaks it may be characterized by $\Delta = R_{\rm axis}-R_0$78, $\Delta = R_{\rm axis}-R_0$79, or $\Delta = R_{\rm axis}-R_0$80; in rotating equilibria it acquires explicit Mach-number dependence; in energetic-particle-rich plasmas it scales with the total MHD pressure; and in stellarators it is naturally recast as a weighted average controlled by Pfirsch–Schlüter current. Across these formulations, the common content is the same: finite pressure, current, and flow reshape toroidal equilibria by displacing flux surfaces outward.