A 1-D Power Balance Model for an ECRH-Heated Fusion Reactor
A radially resolved (1-D) power balance for a steady-state D-T reactor heated only by
electron cyclotron resonance heating. Prescribed density and temperature profiles are integrated
over a toroidal volume to give fusion, bremsstrahlung and synchrotron power densities; an empirical
confinement scaling closes the plasma balance and fixes the ECRH power required to sustain the burn.
That result is then carried into a whole-plant balance in the style of the PROCESS systems
code [1]: blanket multiplication and thermal conversion on the output side, and on the
recirculating side the gyrotron wall plug, the 4.5 K cryoplant, coil and current-lead losses,
torus vacuum pumping, coolant pumping and balance of plant. The figure of merit is the engineering
gain $Q_\mathrm{eng}$, and the notebook is interactive so the recirculating budget can be traded
against machine size and field.
1. Why the plant balance, not the plasma balance, is the hard constraint
Lawson's original criterion asks only that the fusion power exceed what is needed to keep the plasma hot [2]. It is a necessary condition and a very weak one: it says nothing about the electricity a reactor spends on itself. A magnetic fusion plant runs a large parasitic load — heating systems, the cryoplant that holds tens of thousands of tonnes of structure near 4.5 K, coil supplies, torus and cryostat pumping, coolant circulation — and that load is charged against gross electrical output, not against fusion power. The useful figure of merit is therefore the engineering gain
$$ Q_\mathrm{eng} \;=\; \frac{P_\mathrm{gross}}{P_\mathrm{recirc}}, \qquad P_\mathrm{net} \;=\; P_\mathrm{gross} - P_\mathrm{recirc}, $$
and a plant is electrically self-sufficient only for $Q_\mathrm{eng} > 1$ [3].
Two choices set this notebook apart from a textbook 0-D estimate.
It is 1-D. Fusion power density goes roughly as $n^2 \langle\sigma v\rangle(T)$, and $\langle\sigma v\rangle$ is steeply nonlinear in $T$ over the reactor-relevant range, so a volume-averaged calculation systematically misestimates the fusion power for peaked profiles. Bremsstrahlung ($\propto n^2\sqrt{T}$) and synchrotron ($\propto n T^{2.5}B^{2.5}$, roughly) weight the profile differently again. Each term is therefore evaluated locally on a radial grid and integrated over flux-surface volume.
ECRH is the only heating system. That is the natural choice for a current-free configuration — a stellarator needs no current drive, so its recirculating power is dominated by cryogenics and balance of plant rather than by a neutral beam system [4]. It also has a hard cost: the whole auxiliary heating budget passes through gyrotrons at a wall-plug efficiency of order 50% [5], so every megawatt the plasma needs costs the plant about two.
The plant-side accounting — which loads are recirculating, what fraction of neutron power is
recovered as primary heat, and how thermal power converts to gross electric — follows the
structure of the PROCESS systems code [6][1].
using CairoMakie, Printf
use_slate_theme!()
# Physical constants (SI unless noted; temperatures in keV throughout)
const keV = 1.602_176_634e-16 # J per keV
const MW = 1e6
const E_DT = 17.59 # MeV released per D-T reaction
const f_α = 3.52 / 17.59 # fraction of E_DT carried by the alpha
const f_n = 1 - f_α # fraction carried by the 14.1 MeV neutron
const T_amb = 300.0 # K, cryoplant heat-rejection temperature
nothing2. D-T reactivity
The Maxwell-averaged D-T reactivity uses the Bosch-Hale parameterisation [7], accurate to better than 0.25% over 0.2-100 keV. With $\theta$ the Gamow-corrected temperature and $\xi = (B_G^2/4\theta)^{1/3}$,
$$ \langle\sigma v\rangle \;=\; C_1\,\theta\,\sqrt{\frac{\xi}{m_r c^2 T^3}}\;e^{-3\xi}. $$
This is the only steeply nonlinear ingredient in the model, and it is what makes the 1-D treatment worth the trouble.
"""
sigmav_DT(T) -> m³/s
Maxwell-averaged D-T fusion reactivity, Bosch-Hale parameterisation (Bosch & Hale 1992,
Table VII / Eqs. 12-13). `T` is the ion temperature in keV; valid for 0.2 ≤ T ≤ 100 keV.
"""
function sigmav_DT(T)
T <= 0 && return 0.0
B_G = 34.3827 # Gamow constant, √keV
mrc2 = 1.124656e6 # reduced mass energy, keV
C1, C2, C3 = 1.17302e-9, 1.51361e-2, 7.51886e-2
C4, C5 = 4.60643e-3, 1.35000e-2
C6, C7 = -1.06750e-4, 1.36600e-5
θ = T / (1 - T*(C2 + T*(C4 + T*C6)) / (1 + T*(C3 + T*(C5 + T*C7))))
ξ = (B_G^2 / (4θ))^(1/3)
sv = C1 * θ * sqrt(ξ / (mrc2 * T^3)) * exp(-3ξ) # cm³/s
return sv * 1e-6 # → m³/s
end
# Sanity check against the standard tabulated values
(; at_10keV = sigmav_DT(10.0), at_20keV = sigmav_DT(20.0), peak_at = argmax(sigmav_DT, 1:0.5:120))(at_10keV = 1.1361654705836232e-22, at_20keV = 4.330201895538215e-22, peak_at = 66.5)
3. Machine and plasma parameters
The reference point below is a HELIAS-class stellarator power plant: major radius ~18 m, on-axis field ~5 T, moderately peaked profiles [4]. Every downstream quantity is reactive, so the whole balance re-solves as these move.
4. Radial profiles and volume integration
On the normalised minor radius $\rho = r/a$ the profiles are the usual power-law forms
$$ n_e(\rho) = n_{e0}\,(1-\rho^2)^{\alpha_n}, \qquad T(\rho) = T_0\,(1-\rho^2)^{\alpha_T}, $$
with $T_i = T_e$ assumed. For an elongated cross-section of area $\pi a^2 \kappa$ the volume element on a flux surface is $dV/d\rho = 4\pi^2 R_0 a^2 \kappa\,\rho$, so $V = 2\pi^2 R_0 a^2 \kappa$. Any power density $p(\rho)$ integrates as $P = \int_0^1 p(\rho)\,(dV/d\rho)\,d\rho$.
The fuel is diluted by helium ash and the seeded impurity: with $f_j = n_j/n_e$, $n_\mathrm{DT} = n_e\,(1 - 2 f_\mathrm{He} - Z\,f_\mathrm{imp})$ and $n_D = n_T = n_\mathrm{DT}/2$.
const NR = 401
ρ = range(0, 1; length=NR)
dρ = step(ρ)
Zimp = Z_imp.value # unwrap the dropdown Choice
ne = @. ne0 * 1e20 * (1 - ρ^2)^αn # m⁻³
Tp = @. T0 * (1 - ρ^2)^αT # keV (T_i = T_e)
# Fuel dilution and Z_eff
dilution = 1 - 2f_He - Zimp * f_imp
n_DT = ne .* dilution
n_D = n_DT ./ 2
n_T = n_DT ./ 2
Z_eff = dilution + 4f_He + Zimp^2 * f_imp
# Toroidal volume element for an elliptical cross-section
V_plasma = 2π^2 * R0 * aa^2 * κ
dVdρ = @. 4π^2 * R0 * aa^2 * κ * ρ
"Integrate a power density [W/m³] over the plasma volume → W."
volint(p) = sum(p .* dVdρ) * dρ
ne_avg = volint(ne) / V_plasma
T_avg = volint(Tp) / V_plasma
(; V_plasma, ne_avg, ne_bar = ne_avg/1e20, T_avg, Z_eff, dilution)(V_plasma = 1421.2230337568676, ne_avg = 1.691646272214501e20, ne_bar = 1.691646272214501, T_avg = 8.181790884091384, Z_eff = 2.6100000000000003, dilution = 0.8300000000000001)
let
fig = Figure(size=(860, 330))
ax1 = Axis(fig[1,1]; xlabel="ρ = r/a", ylabel="nₑ [10²⁰ m⁻³]", title="Density")
lines!(ax1, ρ, ne ./ 1e20; linewidth=3)
lines!(ax1, ρ, n_DT ./ 1e20; linewidth=2, linestyle=:dash, label="D+T (diluted)")
axislegend(ax1; position=:lb, framevisible=false)
ax2 = Axis(fig[1,2]; xlabel="ρ = r/a", ylabel="T [keV]", title="Temperature")
lines!(ax2, ρ, Tp; linewidth=3, color=Cycled(2))
ax3 = Axis(fig[1,3]; xlabel="ρ = r/a", ylabel="⟨σv⟩ [10⁻²² m³/s]",
title="D-T reactivity on the profile")
lines!(ax3, ρ, sigmav_DT.(Tp) ./ 1e-22; linewidth=3, color=Cycled(3))
ylims!(ax1, 0, nothing); ylims!(ax2, 0, nothing); ylims!(ax3, 0, nothing)
fig
end
5. Local power densities
Fusion. $p_\mathrm{fus} = n_D n_T \langle\sigma v\rangle(T_i)\,E_\mathrm{DT}$, split 20/80 between the alpha (which stays in the plasma) and the 14.1 MeV neutron (which leaves it).
Bremsstrahlung. $p_\mathrm{br} = 5.35\times10^{-37}\,Z_\mathrm{eff}\,n_e^2\sqrt{T_e}$ W m⁻³ with a first-order relativistic correction $(1 + 2T_e/m_ec^2)$ [8]. Line radiation from the seeded impurity is not modelled separately — it is folded into $Z_\mathrm{eff}$, which is the usual systems-code simplification and is why the impurity choice matters here.
Synchrotron. Cyclotron emission is reabsorbed and wall-reflected, so it depends on the reflectivity $r$ of the first wall. The Trubnikov form [9], as commonly implemented in systems codes and superseded in accuracy by [10], is
$$ p_\mathrm{syn} = 6.25\,B_0^{2.5} T_e^{2.5} n_{20}^{1/2} a^{-1/2}(1-r)^{1/2} \Big[1 + 0.12\,T_e\,(n_{20}(1-r))^{-1/2}\Big]^{-3/2} \;\;\mathrm{W\,m^{-3}}. $$
Its $B^{2.5}T^{2.5}$ dependence is why synchrotron, negligible in today's devices, becomes a real term in a high-field, high-temperature reactor.
E_DT_J = E_DT * 1e6 * 1.602_176_634e-19 # J per reaction
# --- fusion -----------------------------------------------------------------
p_fus = @. n_D * n_T * sigmav_DT(Tp) * E_DT_J # W/m³
p_α = f_α .* p_fus
p_n = f_n .* p_fus
# --- bremsstrahlung ---------------------------------------------------------
p_br = @. 5.35e-37 * Z_eff * ne^2 * sqrt(Tp) * (1 + 2Tp/511)
# --- synchrotron (Trubnikov, wall-reflected) ---------------------------------
function p_sync_local(Te, n20, B, a, r)
(Te <= 0 || n20 <= 0) && return 0.0
num = 6.25 * B^2.5 * Te^2.5 * sqrt(n20) * (1-r)^0.5 / sqrt(a)
den = (1 + 0.12 * Te / sqrt(n20 * (1-r)))^1.5
return num / den
end
p_sy = @. p_sync_local(Tp, ne/1e20, B0, aa, r_wall)
p_rad = p_br .+ p_sy
# --- volume integrals -------------------------------------------------------
P_fus = volint(p_fus)
P_alpha = volint(p_α)
P_neut = volint(p_n)
P_brems = volint(p_br)
P_sync = volint(p_sy)
P_rad = P_brems + P_sync
(; P_fus_MW = P_fus/MW,
P_alpha_MW = P_alpha/MW,
P_neut_MW = P_neut/MW,
P_brems_MW = P_brems/MW,
P_sync_MW = P_sync/MW,
rad_frac = P_rad/P_fus,
peak_MW_m3 = maximum(p_fus)/MW,
wall_load_MW_m2 = P_neut/MW / (4π^2 * R0 * aa * sqrt((1+κ^2)/2)))(P_fus_MW = 3125.940578346412, P_alpha_MW = 625.5435381341314, P_neut_MW = 2500.3970402122814, P_brems_MW = 193.67497118980282, P_sync_MW = 17.96786142977199, rad_frac = 0.06770532814527508, peak_MW_m3 = 8.736628192444886, wall_load_MW_m2 = 1.7593276922924055)
6. Closing the plasma balance: how much ECRH does the burn need?
Prescribed profiles fix the stored energy $W = \tfrac32\int (n_eT_e + n_iT_i)\,dV$, the alpha heating and the radiated power. What they do not fix is whether the machine can actually hold that plasma — that is the job of a confinement scaling. For a current-free configuration the relevant one is ISS04, regressed over the International Stellarator Database [11]:
$$ \tau_E^{\mathrm{ISS04}} = 0.134\;a^{2.28} R^{0.64} P^{-0.61} \bar n_{19}^{0.54} B^{0.84} \iota_{2/3}^{0.41}, $$
renormalised by a configuration factor $H$ (about 0.9 for W7-AS, ~1.2 assumed for an optimised HELIAS). Steady state then requires
$$ \underbrace{P_\alpha + P_\mathrm{ECRH}}_{\text{heating}} \;=\; \underbrace{P_\mathrm{rad}}_{\text{brems + sync}} + \underbrace{W/\tau_E}_{\text{transport}}, $$
and since $\tau_E$ itself degrades with heating power, this is implicit in $P_\mathrm{ECRH}$ and is solved by bisection. A negative root means the configuration is over-ignited at this operating point: alpha heating alone exceeds the losses, and the machine would run at a small control-level ECRH power instead.
# Stored thermal energy: 3/2 ∫ (nₑTₑ + n_iT_i) dV
n_i = ne .* (dilution + f_He + f_imp)
W_th = 1.5 * volint((ne .+ n_i) .* Tp .* keV) # J
"ISS04 stellarator energy confinement time [s]; `P_MW` is the net heating power."
tauE_ISS04(P_MW) = H_fac * 0.134 * aa^2.28 * R0^0.64 * max(P_MW, 1e-3)^(-0.61) *
(ne_avg/1e19)^0.54 * B0^0.84 * iota^0.41
# Steady state: P_α + P_ECRH = P_rad + W/τ_E(P_α + P_ECRH)
residual(P_ec) = (W_th / tauE_ISS04((P_alpha + P_ec*MW)/MW)) + P_rad - P_alpha - P_ec*MW
P_ec_req = let lo = -0.999 * P_alpha/MW, hi = 5000.0
if residual(lo) * residual(hi) > 0
NaN
else
for _ in 1:200
mid = (lo + hi)/2
residual(lo)*residual(mid) <= 0 ? (hi = mid) : (lo = mid)
end
(lo + hi)/2
end
end
ignited = P_ec_req < P_ec_min
P_ECRH = max(P_ec_req, P_ec_min) # MW delivered to the plasma
P_heat = P_alpha + P_ECRH*MW
tau_E = tauE_ISS04(P_heat/MW)
P_trans = W_th / tau_E
Q_plasma = P_fus / (P_ECRH*MW)
(; W_th_MJ = W_th/1e6,
tau_E,
P_ec_required_MW = P_ec_req,
P_ECRH_MW = P_ECRH,
ignited,
P_transport_MW = P_trans/MW,
Q_plasma,
triple_product = ne_avg * T_avg * tau_E / 1e21)(W_th_MJ = 1014.4231687026526, tau_E = 1.6335980912623023, P_ec_required_MW = 207.0740459426973, P_ECRH_MW = 207.0740459426973, ignited = false, P_transport_MW = 620.9747514572539, Q_plasma = 15.095762311088665, triple_product = 2.261013464757069)
let
fig = Figure(size=(860, 380))
ax = Axis(fig[1,1]; xlabel="ρ = r/a", ylabel="power density [MW/m³]",
yscale=log10, title="Local power densities")
flr(v) = max.(v ./ MW, 1e-6)
lines!(ax, ρ, flr(p_fus); linewidth=3, label="fusion (total)")
lines!(ax, ρ, flr(p_α); linewidth=2, linestyle=:dash, label="alpha heating")
lines!(ax, ρ, flr(p_br); linewidth=2.5, label="bremsstrahlung")
lines!(ax, ρ, flr(p_sy); linewidth=2.5, label="synchrotron")
ylims!(ax, 1e-5, 20)
axislegend(ax; position=:lb, framevisible=false)
ax2 = Axis(fig[1,2]; xlabel="ρ = r/a", ylabel="dP/dρ [MW]",
title="Volume-weighted contribution")
lines!(ax2, ρ, p_fus .* dVdρ ./ MW; linewidth=3, label="fusion")
lines!(ax2, ρ, p_rad .* dVdρ ./ MW; linewidth=2.5, color=Cycled(3), label="radiation")
axislegend(ax2; position=:lt, framevisible=false)
ylims!(ax2, 0, nothing)
fig
end
7. The plant balance
7.1 Output side
Neutrons deposit in the breeding blanket with an energy multiplication $M_n \approx 1.1-1.3$ from exothermic $^6$Li(n,α)T and Be(n,2n) reactions. Everything the plasma itself loses — alpha heating plus absorbed ECRH — eventually lands on the first wall and divertor, and with a helium- or water-cooled blanket it is recovered at primary-coolant temperature. Coolant pumping work is dissipated into the coolant and so is recovered too. The thermal power and gross electric are
$$ P_\mathrm{th} = M_n P_n + P_\alpha + P_\mathrm{ECRH} + P_\mathrm{pump}, \qquad P_\mathrm{gross} = \eta_\mathrm{th}\,P_\mathrm{th}, $$
with $\eta_\mathrm{th}$ set by the blanket coolant — about 31% for water and 41% for helium at a 500 °C top temperature [1].
7.2 Recirculating side
ECRH. The dominant controllable load. Wall plug to plasma is the product of gyrotron efficiency (~50% with a depressed collector), transmission-line efficiency (~90%) and single-pass absorption [5]: $P_\mathrm{ECRH}^\mathrm{wall} = P_\mathrm{ECRH}/(\eta_g\eta_t\eta_a)$.
Cryogenics. A 4.5 K refrigerator costs about $(T_\mathrm{amb}/T_c - 1)/\eta_2 \approx 260$ W of compressor power per watt lifted, taking a second-law efficiency $\eta_2 \approx 25\%$ typical of large helium plants [12][13]. The cold load has three parts: nuclear heating in the coils (what survives the shield), static heat leak through the cryostat and supports, and the conduction/Joule load of the current leads. The nuclear term is the one that scales with the reactor, which is why shield thickness is a cryoplant design variable as much as a neutronics one.
Coils. A superconducting winding dissipates nothing in steady state, so what is left is power supply and rectifier loss plus any normal-conducting trim and control coils — small, but not zero, and a fully resistive machine would swamp everything else here.
Vacuum, tritium plant and balance of plant. Torus exhaust and cryostat pumping, the tritium loop, and site services. These are weak functions of the physics and are entered directly.
# ---- output side (MW) -------------------------------------------------------
P_neut_MW = P_neut / MW
P_alpha_MW = P_alpha / MW
P_th_bare = M_n * P_neut_MW + P_alpha_MW + P_ECRH # before pumping work
P_pump = f_pump * P_th_bare / (1 - f_pump) # pumping work is recovered in the coolant
P_th = P_th_bare + P_pump
P_gross = η_th * P_th
# ---- recirculating side (MW) -----------------------------------------------
# 1. ECRH wall plug
P_ecrh_wall = P_ECRH / (η_gyro * η_trans * η_abs)
# 2. Cryoplant. Cold load [kW] → compressor power via a fraction of Carnot.
A_cryostat = 4π^2 * (R0 + 1.0) * (aa + 2.5) # ≈ torus surface at the cryostat
q_nuclear = f_nuc * 1e-5 * P_neut_MW * 1e3 # kW at T_cold
q_static = h_stat * A_cryostat / 1e3 # kW
q_leads = w_lead * I_leads / 1e3 # kW
Q_cold = q_nuclear + q_static + q_leads # kW
cop_inv = (T_amb / T_cold - 1) / η_carnot # W of compressor per W lifted
P_cryo = Q_cold * cop_inv / 1e3 # MW
# 3-5. coils, vacuum/tritium, balance of plant
P_bop = f_bop * P_gross
P_recirc = P_ecrh_wall + P_cryo + P_coil + P_vac + P_pump + P_bop
P_net = P_gross - P_recirc
Q_eng = P_gross / P_recirc
f_recirc = P_recirc / P_gross
(; P_th, P_gross, P_recirc, P_net, Q_eng, f_recirc,
Q_cold_kW = Q_cold, cryo_W_per_W = cop_inv,
loads = (; ECRH=P_ecrh_wall, cryo=P_cryo, coils=P_coil,
vacuum=P_vac, pumping=P_pump, bop=P_bop))(P_th = 4134.520318609094, P_gross = 1695.1533306297283, P_recirc = 948.1924336328518, P_net = 746.9608969968765, Q_eng = 1.7877735262398289, f_recirc = 0.5593549660080661, Q_cold_kW = 101.14850302004302, cryo_W_per_W = 262.6666666666667, loads = (ECRH = 484.38373319929195, cryo = 26.56834012659797, coils = 6.0, vacuum = 12.0, pumping = 351.4342270817729, bop = 67.80613322518913))
7.3 Where the power actually goes
The Sankey below is the same accounting as the PROCESS power-flow chart [1, Fig. 9].
Read it left to right: fusion power splits 20/80 into alphas and neutrons, the blanket adds
multiplication energy, everything recoverable becomes primary heat, and the thermal cycle then has
to pay the recirculating loads out of gross electric before anything reaches the grid.
let
P_mult = (M_n - 1) * P_neut_MW
P_reject = P_th - P_gross
r(x) = round(x; digits=1)
links = [
("D-T fusion", "Neutrons", r(P_neut_MW)),
("D-T fusion", "Alphas", r(P_alpha_MW)),
("Neutrons", "Primary heat", r(P_neut_MW + P_mult)),
("Blanket Mₙ", "Primary heat", r(P_mult)),
("Alphas", "Primary heat", r(P_alpha_MW)),
("ECRH → plasma", "Primary heat", r(P_ECRH)),
("Pumping work", "Primary heat", r(P_pump)),
("Primary heat", "Gross electric", r(P_gross)),
("Primary heat", "Rejected heat", r(P_reject)),
("Gross electric", "Net to grid", r(max(P_net, 0.0))),
("Gross electric", "Recirculating", r(P_recirc)),
("Recirculating", "Gyrotron plant", r(P_ecrh_wall)),
("Recirculating", "Coolant pumps", r(P_pump)),
("Recirculating", "Balance of plant", r(P_bop)),
("Recirculating", "Cryoplant", r(P_cryo)),
("Recirculating", "Vacuum + tritium", r(P_vac)),
("Recirculating", "Coil supplies", r(P_coil)),
]
echart(:sankey, links;
title = @sprintf("Plant power flow [MW] — Q_eng = %.2f, P_net = %.0f MW", Q_eng, P_net),
height = 520,
label = (color="#c8d0e0", fontSize=12),
emphasis = (focus="adjacency",))
end8. Sweeping the design space
Everything above is packaged into one function so the balance can be re-solved off the reactive grid. The sweeps that follow answer the three questions this kind of model exists to answer: how large does the machine have to be, how good does the gyrotron plant have to be, and how much does the cryoplant actually cost you.
"""
solve_plant(; kwargs...) -> NamedTuple
Full 1-D plasma + plant power balance. Every keyword defaults to the value of the corresponding
notebook control, so `solve_plant()` reproduces the reactive result and `solve_plant(R0=12.0)`
re-solves with only that changed.
"""
function solve_plant(; R0=R0, aa=aa, κ=κ, B0=B0, ne0=ne0, T0=T0, αn=αn, αT=αT,
f_imp=f_imp, Zimp=Zimp, f_He=f_He, H_fac=H_fac, iota=iota,
r_wall=r_wall, P_ec_min=P_ec_min,
M_n=M_n, η_th=η_th, f_pump=f_pump,
η_gyro=η_gyro, η_trans=η_trans, η_abs=η_abs,
f_nuc=f_nuc, h_stat=h_stat, I_leads=I_leads, w_lead=w_lead,
η_carnot=η_carnot, T_cold=T_cold,
P_coil=P_coil, P_vac=P_vac, f_bop=f_bop)
# --- profiles and geometry ---
ne = @. ne0 * 1e20 * (1 - ρ^2)^αn
Tp = @. T0 * (1 - ρ^2)^αT
dil = 1 - 2f_He - Zimp * f_imp
Zef = dil + 4f_He + Zimp^2 * f_imp
nDT = ne .* dil
n_i = ne .* (dil + f_He + f_imp)
V = 2π^2 * R0 * aa^2 * κ
dV = @. 4π^2 * R0 * aa^2 * κ * ρ
vi(p) = sum(p .* dV) * dρ
navg = vi(ne) / V
# --- local power densities ---
P_fus = vi(@. (nDT/2)^2 * sigmav_DT(Tp) * E_DT_J)
P_a = f_α * P_fus
P_n = f_n * P_fus
P_rad = vi(@. 5.35e-37 * Zef * ne^2 * sqrt(Tp) * (1 + 2Tp/511)) +
vi(@. p_sync_local(Tp, ne/1e20, B0, aa, r_wall))
# --- plasma balance → required ECRH ---
W = 1.5 * vi((ne .+ n_i) .* Tp .* keV)
τ(P) = H_fac * 0.134 * aa^2.28 * R0^0.64 * max(P, 1e-3)^(-0.61) *
(navg/1e19)^0.54 * B0^0.84 * iota^0.41
res(Pe) = W / τ((P_a + Pe*MW)/MW) + P_rad - P_a - Pe*MW
lo, hi = -0.999 * P_a/MW, 5000.0
Pec_req = if res(lo) * res(hi) > 0
NaN
else
for _ in 1:200
m = (lo + hi)/2
res(lo)*res(m) <= 0 ? (hi = m) : (lo = m)
end
(lo + hi)/2
end
P_ec = max(Pec_req, P_ec_min)
# --- plant ---
Pn_MW, Pa_MW = P_n/MW, P_a/MW
th_bare = M_n*Pn_MW + Pa_MW + P_ec
P_pmp = f_pump * th_bare / (1 - f_pump)
P_th = th_bare + P_pmp
P_gr = η_th * P_th
P_ew = P_ec / (η_gyro * η_trans * η_abs)
Acryo = 4π^2 * (R0 + 1.0) * (aa + 2.5)
Qc = f_nuc*1e-5*Pn_MW*1e3 + h_stat*Acryo/1e3 + w_lead*I_leads/1e3 # kW
P_cry = Qc * (T_amb/T_cold - 1)/η_carnot / 1e3
P_bp = f_bop * P_gr
P_rec = P_ew + P_cry + P_coil + P_vac + P_pmp + P_bp
(; P_fus = P_fus/MW, P_ECRH = P_ec, P_ec_req = Pec_req,
Q_plasma = P_fus/MW/P_ec, tau_E = τ((P_a + P_ec*MW)/MW),
P_th, P_gross = P_gr, P_recirc = P_rec,
P_net = P_gr - P_rec, Q_eng = P_gr/P_rec,
Q_cold_kW = Qc, V,
ECRH = P_ew, cryo = P_cry, coils = P_coil,
vacuum = P_vac, pumping = P_pmp, bop = P_bp)
end
# must reproduce the reactive cells
let s = solve_plant()
(; s.P_fus, s.P_ECRH, s.Q_eng, s.P_net, matches = isapprox(s.Q_eng, Q_eng; rtol=1e-9))
end(P_fus = 3125.940578346412, P_ECRH = 207.0740459426973, Q_eng = 1.7877735262398284, P_net = 746.9608969968763, matches = true)
let
Rs = range(10.0, 22.0; length=60)
S = [solve_plant(R0=R) for R in Rs]
fig = Figure(size=(880, 640))
ax1 = Axis(fig[1,1]; xlabel="R₀ [m]", ylabel="Q_eng", title="Engineering gain vs machine size")
lines!(ax1, Rs, getfield.(S, :Q_eng); linewidth=3)
hlines!(ax1, [1.0]; linestyle=:dash, color=:gray, linewidth=2)
text!(ax1, 10.3, 1.05; text="electrical break-even", color=:gray, fontsize=11)
ylims!(ax1, 0, nothing)
ax2 = Axis(fig[1,2]; xlabel="R₀ [m]", ylabel="MW", title="Gross, net and recirculating")
lines!(ax2, Rs, getfield.(S, :P_gross); linewidth=3, label="gross electric")
lines!(ax2, Rs, getfield.(S, :P_recirc); linewidth=3, label="recirculating")
lines!(ax2, Rs, getfield.(S, :P_net); linewidth=3, label="net to grid")
hlines!(ax2, [0.0]; color=:gray, linestyle=:dot)
axislegend(ax2; position=:lt, framevisible=false)
ax3 = Axis(fig[2,1:2]; xlabel="R₀ [m]", ylabel="recirculating power [MW]",
title="Recirculating budget breakdown")
keys_ = (:ECRH, :cryo, :pumping, :bop, :vacuum, :coils)
labs = ("gyrotron plant", "cryoplant", "coolant pumps", "balance of plant",
"vacuum + tritium", "coil supplies")
lower = zeros(length(Rs))
for (k, lab) in zip(keys_, labs)
upper = lower .+ getfield.(S, k)
band!(ax3, Rs, lower, upper; label=lab)
lower = upper
end
lines!(ax3, Rs, getfield.(S, :P_gross); color=:white, linewidth=2.5,
linestyle=:dash, label="gross electric")
axislegend(ax3; position=:lt, framevisible=false, nbanks=2)
fig
end
let
fig = Figure(size=(880, 640))
# (a) gyrotron wall-plug efficiency, for a few confinement qualities
ax1 = Axis(fig[1,1]; xlabel="η gyrotron wall-plug", ylabel="Q_eng",
title="ECRH plant efficiency")
ηs = range(0.30, 0.70; length=50)
for H in (1.1, 1.25, 1.4, 1.6)
lines!(ax1, ηs, [solve_plant(η_gyro=e, H_fac=H).Q_eng for e in ηs];
linewidth=2.5, label="H = $H")
end
axislegend(ax1; position=:lt, framevisible=false)
# (b) coil temperature — the HTS argument, priced in recirculating power
ax2 = Axis(fig[1,2]; xlabel="coil operating temperature [K]", ylabel="cryoplant [MW]",
title="Cryoplant load vs coil temperature")
Ts = range(4.5, 20.0; length=60)
for fn in (2.0, 4.0, 10.0)
lines!(ax2, Ts, [solve_plant(T_cold=t, f_nuc=fn).cryo for t in Ts];
linewidth=2.5, label="nuclear heating $(fn)e-5·Pₙ")
end
axislegend(ax2; position=:rt, framevisible=false)
# (c) coolant pumping fraction: helium vs water
ax3 = Axis(fig[2,1]; xlabel="coolant pumping fraction of P_th", ylabel="Q_eng",
title="Coolant choice")
fs = range(0.0, 0.12; length=60)
lines!(ax3, fs, [solve_plant(f_pump=f).Q_eng for f in fs]; linewidth=3)
vlines!(ax3, [0.0005]; linestyle=:dash, color=:gray)
vlines!(ax3, [0.085]; linestyle=:dash, color=:gray)
text!(ax3, 0.004, 1.35; text="water", color=:gray, fontsize=11)
text!(ax3, 0.074, 1.35; text="helium", color=:gray, fontsize=11)
# (d) confinement quality → required ECRH → everything
ax4 = Axis(fig[2,2]; xlabel="H confinement factor", ylabel="Q_eng",
title="Confinement quality")
Hs = range(0.9, 1.9; length=60)
lines!(ax4, Hs, [solve_plant(H_fac=H).Q_eng for H in Hs]; linewidth=3, color=Cycled(4))
for ax in (ax1, ax3, ax4)
hlines!(ax, [1.0]; linestyle=:dot, color=:gray, linewidth=2)
end
fig
end
let
nn = range(1.0, 3.6; length=70) # nₑ(0) in 10²⁰ m⁻³
tt = range(8.0, 30.0; length=70) # T(0) in keV
S = [solve_plant(ne0=n, T0=t) for n in nn, t in tt]
Qe = getfield.(S, :Q_eng)
Pnet = getfield.(S, :P_net)
Pfus = getfield.(S, :P_fus)
ign = getfield.(S, :P_ec_req)
fig = Figure(size=(880, 400))
ax = Axis(fig[1,1]; xlabel="nₑ(0) [10²⁰ m⁻³]", ylabel="T(0) [keV]",
title="Engineering gain over the operating space")
hm = heatmap!(ax, nn, tt, Qe; colormap=:viridis, colorrange=(0, 4))
contour!(ax, nn, tt, Qe; levels=[1.0], color=:white, linewidth=3)
contour!(ax, nn, tt, Qe; levels=[1.5, 2.0, 2.5, 3.0], color=(:white, 0.4), linewidth=1)
contour!(ax, nn, tt, ign; levels=[P_ec_min], color=:orangered, linewidth=2.5, linestyle=:dash)
scatter!(ax, [ne0], [T0]; color=:red, markersize=14, strokecolor=:white, strokewidth=1.5)
Colorbar(fig[1,2], hm; label="Q_eng")
ax2 = Axis(fig[1,3]; xlabel="nₑ(0) [10²⁰ m⁻³]", ylabel="T(0) [keV]",
title="Net electric [MW] and fusion power")
hm2 = heatmap!(ax2, nn, tt, Pnet; colormap=:magma)
contour!(ax2, nn, tt, Pfus; levels=[1000, 2000, 3000, 4000, 5000],
color=(:cyan, 0.7), linewidth=1.5)
scatter!(ax2, [ne0], [T0]; color=:red, markersize=14, strokecolor=:white, strokewidth=1.5)
Colorbar(fig[1,4], hm2; label="P_net [MW]")
fig
end
9. The operating point
| Plasma volume | 1421.0 m³ (R₀ = 18.0 m, a = 2.0 m, B₀ = 5.0 T) |
| Volume-average density / temperature | 1.69 ×10²⁰ m⁻³ / 8.2 keV |
| Energy confinement time | 1.63 s (ISS04 × H = 1.4) |
| Fusion power | 3126 MW |
| Neutron wall load | 1.76 MW/m² |
| Radiated fraction | 6.8% of P_fus |
| ECRH to plasma | 207.1 MW → Q_plasma = 15.1 |
| Cryogenic heat load at 4.5 K | 101.1 kW, costing 263.0 W/W |
| Gross electric | 1695 MW |
| Recirculating | 948 MW (55.9% of gross) |
| Net to grid | 747 MW |
| Engineering gain | Q_eng = 1.79 |
10. What the model says
The plasma is not the binding constraint; the recirculating budget is. At this operating point the plasma runs at $Q_\mathrm{plasma} = 15.1$ — comfortably burning — yet $Q_\mathrm{eng}$ is only 1.79, because 56.0% of gross electric never leaves the site. Scanning $R_0$ barely moves $Q_\mathrm{eng}$: fusion power and every major parasitic load scale together, so buying a bigger machine buys net megawatts, not efficiency. This is the same conclusion systems-code studies reach for both tokamaks and stellarators [1][3][14].
With ECRH as the sole heating system, gyrotron efficiency is a plant-level parameter. The gyrotron plant draws 484 MW here, the largest single recirculating item. Because the chain efficiency is only $\eta_g\eta_t\eta_a \approx$ 0.43, each megawatt the plasma needs costs about 2.3 MW at the wall plug. Moving the gyrotron from 40% to 60% — the gap multi-stage depressed collectors are meant to close [5] — is worth more $Q_\mathrm{eng}$ than several metres of major radius.
Confinement quality is the strongest single lever, and it saturates. $Q_\mathrm{eng}$ climbs steeply with $H$ until the configuration ignites; past that point ECRH drops to its control level and the curve flattens. Everything below the ignition knee is being paid for twice — once in gyrotron wall plug, once in the gross electric that never reaches the grid.
Cryogenics is real but second-order, and it is where higher-temperature superconductors pay. The cold load is only 101.0 kW, but at 263.0 W of compressor power per watt lifted it becomes 26.6 MW of electricity. Raising the coil operating temperature from 4.5 K to 20 K cuts that by roughly a factor of four, since the Carnot penalty falls as $T_\mathrm{amb}/T_c - 1$ — a systems-level argument for HTS that is independent of the field-strength argument usually made for it.
Coolant choice is not a detail. Helium pumping at ~8.5% of thermal power [1] costs 351 MW here — comparable to the entire gyrotron plant. Its higher cycle efficiency does not fully pay for it.
Limitations
Profiles are prescribed rather than solved from a transport equation, so the 1-D treatment captures profile weighting of the source terms but not profile consistency; $\tau_E$ still comes from a global empirical scaling. Impurity line radiation is folded into $Z_\mathrm{eff}$ rather than computed from a coronal-equilibrium cooling curve, which understates radiation for high-Z seeding. Divertor heat handling, the tritium breeding ratio, pulsed operation and cost are all out of scope; these are exactly the things a real systems code adds [6][1].
References
- Kovari, M. and Fox, F. and Harrington, C. and Kembleton, R. and Knight, P. and Lux, H. and Morris, J.. PROCESS: A systems code for fusion power plants --- Part 2: Engineering. 2016.
- Lawson, J. D.. Some Criteria for a Power Producing Thermonuclear Reactor. 1957.
- Menard, J. E. and Brown, T. and El-Guebaly, L. and others. Fusion nuclear science facilities and pilot plants based on the spherical tokamak. 2016.
- Warmer, F. and Beidler, C. D. and Dinklage, A. and Wolf, R.. From W7-X to a HELIAS fusion power plant: motivation and options for an intermediate-step burning-plasma stellarator. 2016.
- Franke, T. and Jelonnek, J. and others. Design considerations for future DEMO gyrotrons: A review on related gyrotron activities within EUROfusion. 2017.
- Kovari, M. and Kemp, R. and Lux, H. and Knight, P. and Morris, J. and Ward, D. J.. PROCESS: A systems code for fusion power plants --- Part 1: Physics. 2014.
- Bosch, H.-S. and Hale, G. M.. Improved formulas for fusion cross-sections and thermal reactivities. 1992.
- Wesson, J.. Tokamaks. 2011.
- Trubnikov, B. A.. Universal coefficients for synchrotron emission from plasma configurations. 1979.
- Albajar, F. and Johner, J. and Granata, G.. Improved calculation of synchrotron radiation losses in realistic tokamak plasmas. 2001.
- Yamada, H. and Harris, J. H. and Dinklage, A. and others. Characterization of energy confinement in net-current free plasmas using the extended International Stellarator Database. 2005.
- Strobridge, T. R.. Cryogenic Refrigerators --- An Updated Survey. 1974.
- Radebaugh, R.. Cryocoolers: the state of the art and recent developments. 2009.
- Costley, A. E. and Hugill, J. and Buxton, P. F.. On the power and size of tokamak fusion pilot plants and reactors. 2015.