Burning plasma achieved in inertial fusion – Nature.com

Terminology

Definitions of commonly used mathematical symbols are summarized in Extended Data Table 2.

Reproducibility

Experiments subsequent to those described in this Article have demonstrated the reproducibility of a burning-plasma state at NIF, with two additional experiments that have performance comparable to the highest-performing experiments in this Article. These newer experiments, N210307 and N210605, were conducted in the Hybrid E platform. N210307 repeated N210207, albeit using a capsule from a different fabrication batch and produced a yield of approximately 145 kJ with an experimentally inferred Q_α = 1.34 ± 0.07 from the Hurricane criterion. Experiment N210605 reduced the thickness of the ice layer relative to N210207 and resulted in a lower yield (135 kJ) but high ion temperature, with Q_α = 1.40 ± 0.10, again from the Hurricane criterion. These additional experiments confirm that the burning-plasma state is reproducible at NIF, and full details and analysis on them will be presented in future publications.

Inferred hot-spot conditions

Hot-spot conditions must be inferred from measured quantities using a model. The simplest hot-spot model is to assume an isobaric volume of uniform conditions, as used in a previous work¹⁴ between equations 2 and 3, in which case the hot-spot number density is given by

where Y is the fusion yield in J, ⟨σν⟩ is the fusion reactivity, which depends on the ion temperature (T_i), V_hs is the hot-spot volume in cm³, and τ is the burn duration in s, for equimolar DT mixtures. The remaining hot-spot quantities follow from the inferred density, including the pressure (P = (1 + Z)nk_BT_i, with k_B Boltzmann’s constant), hot-spot energy (E_hs = 1.5PV_hs), and areal density (({rho R}_{{rm{h}}{rm{s}}}=(2.5n/{N}_{{rm{a}}})sqrt[3]{3{V}_{{rm{h}}{rm{s}}}/4pi })).

A more detailed inference is to use a one-dimensional (1D) profile in radius for temperature and density, maintaining the isobaric assumption. A conduction-limited profile follows the expression⁴⁵:

where T_min is the temperature at the boundary, T₀ is the central temperature and R₀ is the hot-spot boundary. β is the thermal conductivity power law, 2.5 from classical Spitzer conductivity. Following a previous work²⁶ we use a lower value, β = 2/3, which accounts for additional physics, dynamical processes and reproduces radiation–hydrodynamics simulations. The density profile is then determined by the isobaric assumption through P ∝ nT being constant. T_min is taken as 1 keV leaving T₀, R₀ and P as free parameters in the model; the data are compared to synthetic data calculated from this 1D profile with the model parameters adjusted to minimize residuals. As in the zero-dimensional (0D) model, the hot-spot energy simply follows from pressure and volume, and the areal density is the mass density integrated over the inferred radial profile.

In either dimensionality the model’s radius is matched to the experimental measurements, which take a contour of emission level, by calculating synthetic emission images to calculate an equivalent contour radius. The measurements include 2D and three-dimensional (3D) asymmetries, so an equivalent spherical volume, and radius, are calculated using the modal decompositions, where the emission contour measured from the equator (R_eq) and pole (R_po) are

where P₀ and M₀ are the average measured sizes from each view, δP_ℓ (δM_m) is the relative modal amplitude, often referred to as P_ℓ/P₀ with the nomenclature above chosen for clarity. P_ℓ are the Legendre polynomials, and ϕ_m are the orientation of the azimuthal modes. R_hs is the hot-spot radius as a function of both θ and ϕ, which is integrated to obtain the volume. Here, the dominant parameters are P₀, M₀ and δP₂, which are given in ref. ⁷.

Implosion velocity (v_imp) is inferred using a rocket model of the implosion⁴⁶ constrained by both supporting experiments, especially in-flight radiography, and the measured time of peak nuclear production on each experiment. The inferred yield amplification given in Extended Data Table 1 is a function of the measured yield, shell compression and fuel mass (m_fuel); both the velocity and Y_amp inferences use the prescription given in a previous work²⁶. The fuel kinetic energy then follows from (frac{1}{2}{m}_{{rm{fuel}}}{v}_{{rm{imp}}}^{2}). Our techniques for inferring the PdV work done on the fuel are discussed in the following section.

A comparison of inferred values using 0D and 1D models are shown in Extended Data Table 3. Inferred pressures are highly consistent between these calculations, whereas hot-spot energies and areal densities are higher in the 1D model owing to substantial mass near the 1 keV temperature cut-off.

Inferring G
_fuel

The total fusion yield produced by a mass (m) of DT, over a characteristic confinement time, τ, is Y ≈ 5mP_ατ—with P_α = 8.2 × 10²⁴ρ⟨σν⟩ in GJ g⁻¹ s⁻¹ the specific DT fusion power for a given mass density, ρ, of DT with reaction rate ⟨σν⟩—and the internal energy in that DT is E_hs = c_DTmT_i. Therefore, one can write (O.A.H. et al., manuscript in preparation)

with

where P is in Gbar, T_i in keV, and τ in s. In equation (7), the total energy delivered by PdV work, E_PdV,tot, is determined from the hot spot and compressed, but cold, DT fuel energy at stagnation, E_hs and E_fuel, respectively, at peak compression. The last term in the denominator represents a correction for additional energy retained by self-heating of the fuel from α-particle deposition but not then lost as bremsstrahlung. So, E_PdV,tot ≈ E_hs + E_fuel − qY/10, where q is a ‘quality’ factor, 0 ≤ q ≤ 1, measuring the ability of the implosion to retain self-heating energy (O.A.H. et al., manuscript in preparation). Here we use q ≈ 0.7, inferred from simulations, and the factor of 10 results from one-fifth of the fusion energy released as α-particles and half of those produced up until the time of peak fusion burn. Albeit generally arrived at in a different fashion than above, the product P(⟨σν⟩/T²)τ is Lawson’s⁹ parameter for ignition. Figure 2a uses the useful reaction-rate approximation (langle sigma vrangle approx 4.2times {10}^{-20}{T}_{{rm{i}}}^{3.6}) (in units of cm³ s⁻¹ for ion temperature range 3.5 < T_i < 6.5 keV) to simplify the abscissa.

An expression for the fuel gain is given in equation (7). The yield is measured and the hot-spot energy is inferred as described in the previous section. Precisely determining the cold-fuel energy from data is not straightforward. For the purposes of this analysis we actually require the total PdV work done on all the DT. This is at a minimum the fuel kinetic energy and internal energy at peak velocity, which are both inferred. This neglects any work done by the inflowing remaining ablator material on the fuel, which can occur in these implosions. In this case the hot-spot energy is more than half the previous estimate; in this scenario we assume equipartition between the hot spot and cold fuel to evaluate equation (7).

Inferred PdV work

The primary uncertainty in the Betti metric⁸ is in the inference of PdV work on the hot spot. Here we use three methodologies: two inferences using an analytic model, and a direct extraction of PdV work from simulations that match the experimental observables.

We use the hydrodynamic piston model of an implosion described previously³⁴. This analytic model abstracts the implosion process using opposed pistons to represent the imploding shell. In spherical geometry, the stagnation pressure from this mechanical work on the hot spot is given by (equation 24in ref. ³⁴):

where ρδR_ave is the average shell areal density, calculated from the measured neutron ‘down-scattered ratio’ (DSR) using the relation ρδR_ave ≈ 19.3DSR, v_imp is the implosion velocity and R_hs is the average hot-spot radius (which can be obtained from the volume, V_hs, given in Extended Data Table 1). The factor f ² represents the effect of mode-1 asymmetry and is a measure of the residual kinetic energy (kinetic energy that is never converted into internal energy) in the implosion.

From the piston pressure we obtain the hot-spot internal energy (E_hs) from

In the absence of α-heating (which adds energy to the hot-spot) and radiative X-ray losses, or when α-heating exactly balances X-ray losses, then E_hs = E_PdV,hs. For low yield amplification implosions (Y_amp < 1.5), X-ray losses dominate over α-heating energy gains, so E_hs < E_PdV,hs. For higher yield amplification implosions (Y_amp > 2), α-heating energy gains start to dominate over X-ray losses, so E_hs > E_PdV,hs. The estimated values for these four experiments are given in Extended Data Table 4 as the piston methodology.

We can also estimate the stagnated fuel mass in a similar fashion, using

which allows us to then estimate the total mass that stagnates from m_shell + m_hs, with m_hs from the hot-spot inferences described earlier. We then estimate the PdV work from

where m_fuel is the initial fuel mass. The factor of 0.73 is derived from 1D simulations in which the imploding mass stagnates efficiently, and we drop the residual kinetic energy factor f ² because the inferred shell mass does not include non-stagnated material. This estimate leads to smaller estimates of E_PdV,hs than the first empirical estimate, and are given in Extended Data Table 4 as the stagnated mass estimate.

For analysis of previously published campaigns we use the simple relation E_PdV,hs ≈ (0.5–0.7)KE_fuel(1 − f ²), this is easy to evaluate with the available data and the factor 0.5–0.7 accounts for a wide range of 1D to 2D/3D behaviour observed on past experiments. For comparison, the proportionality constant inferred from the first methodology (equation (10)) is between 0.60 and 0.73 for our four experiments.

We also use radiation–hydrodynamics simulations to estimate the PdV work done on these implosions. The first simulation-based methodology is to use 2D simulations with degradation mechanisms that match the observed performance, and interrogate the work done upon the mass elements that form the hot spot to infer E_PdV,hs. The simulation methodology is described in ref. ⁶, and the values of E_PdV,hs for this method are given in Extended Data Table 4. The same fusion performance can be generated with varying application of degradation mechanisms that either degrade E_PdV,hs or do not; an estimate of the 2D simulation uncertainty of ±0.5 kJ is estimated by studying multiple simulations.

A similar energy-balance analysis can be done with 1D simulations, in which the work done upon the hot spot is well defined with a Lagrangian mesh. The 1D simulations are tuned to match the measured yields, but are expected to underestimate E_PdV,hs since they cannot properly incorporate residual kinetic energy. This estimate is given in Extended Data Table 4 as an upper bound.

We have thus develop four methodologies for estimating E_PdV,hs. In the main analysis we use a combination of the empirical piston model estimate as the more pessimistic data-based inference, and use the 2D simulated E_PdV,hs as the most robust computational description of the experiments.

Modified Hurricane metric

At peak burn, the time rate of change of hot-spot volume, dV/dt, is nearly zero, and therefore so is the heating rate, so time integration is needed. Mathematically, a statement of a burning plasma appropriate for ICF is

where t_pf is the time of peak fusion rate, and t_minV is the time of minimum hot-spot volume.

The integrals in equation (13) are easily approximated² without knowing the details of the actual implosion using the mathematical method of steepest descent; assuming that the thermodynamic quantities of interest, such as T, P, ρ, and so on, are impulsive, being highly peaked around the time of stagnation. Ultimately, the solution to equation (13), in terms of only burn-average hot-spot areal density, ρR_hs, T_i and v_imp is equation (1) after a correction to the original derivation.

A recent note from our colleagues at Los Alamos⁴⁷ discovered an arithmetic error in the derivation of the criteria as published in ref. ³. The error is in going from equation 8 to equation 9 in ref. ³, in which the conversion to peak temperature (T₀) to burn-averaged temperature (T_hs) should be, for n ≈ 4,

Additionally, we now believe that the inclusion of the fraction of α-particles stopping in the hot spot (f_α) in ref. ³ was inappropriate. When considering the temperature evolution of a defined mass—for example, the self-heating criterion in equation (17)—this is necessary because f_α is fundamentally the fraction of α-particle energy deposited into that mass. On the other hand, the burning-plasma criteria is one on the energy of the hot spot,

and α-particles that escape the hot spot still contribute to its energy via generation of additional hot-spot mass, as seen by examining the time derivative of the above:

Therefore, the inclusion of f_α in a burning-plasma criterion is inappropriate. We note that not including an f_α factor is consistent with other criteria, for example, ref. ⁸. With these two modifications to the criterion published in ref. ³ we use a new criterion (equation (1)). This modified criterion is slightly more restrictive for the burning-plasma threshold in the regime relevant to these experiments.

Model uncertainties for Hurricane’s metric

The Hurricane metric³ depends on more quantities than the Betti metric, although these quantities are more straightforward to infer than E_PdV,hs. The metric reduces to equation (1) where ρR_hs and v_imp are inferred as described previously, and T_i is measured. ⟨σν⟩ contains some systematic uncertainty from the evaluation used. Data uncertainties are well defined for T_i and in the inference of ρR_hs and v_imp, and are propagated as described in the next section; the inferred ρR_hs can also vary between models, which will be discussed.

Equation (1) depends on the fusion reactivity; in this work we use the evaluation from Bosch and Hale³². Recent publications have presented alternative evaluations⁴⁸ which differ by about 2%. We note that the inferred (rho Rpropto 1/sqrt{langle sigma vrangle }) from equation (2), so the condition in equation (1) depends on the reactivity as (1/sqrt{langle sigma vrangle }). f_α is also weakly increasing with ρR, leading to the condition being slightly less than square-root dependent on ⟨σν⟩, so this criterion has <1% uncertainty from the choice of ⟨σν⟩ evaluation.

The Hurricane criterion is sensitive to the inferred hot-spot ρR, which can vary between models depending on the spatial dependence of ρ. As shown in Extended Data Table 3, the 0D and 1D hot-spot models agree quite well. We also check these values using a 3D reconstruction of the hot-spot density and temperature profiles (a yet unpublished method of L. Divol, but briefly described in ref. ³⁵): for N201101 this gives a value of ρR_hs ≈ 0.36–0.38 g cm⁻² to the 1-keV contour for N201101 and ρR_hs ≈ 0.35–0.36 g cm⁻² for N201122. These values are consistent with the simple models described earlier.

Self-heating regime

The hot-spot per unit mass power balance is:

which describes the temporal evolution of the temperature (T) in terms of the balance of self heating (P_α) versus bremsstrahlung (P_b) and electron conduction (P_e) losses plus PdV work. Here electron conduction losses are calculated relative to a hot-spot boundary that is defined relative to a fraction of the peak burn rate or a specified ion temperature. Thermal conduction cools the hot spot while increasing the mass of the hot spot. Because the fusion burn rate is more strongly dependent on the temperature of the spot than its mass in the temperature range achieved by compression alone, α-heating must provide sufficient heating for the hot-spot temperature to increase in the presence of this conduction into an increasing mass. Hot-spot volume change, dV/dt, is negative on implosion, increasing T. During expansion the PdV term becomes an energy loss term. The bremsstrahlung loss can be enhanced beyond the emission of clean DT by the presence of high-Z contamination of the DT (that is, mix), by a fraction f_b. In equation (17), f_α is the fraction of α-particles stopped in the hot spot, evaluated using fits with modern stopping-power theory³⁶.

Uncertainty analysis

We perform uncertainty analysis for all hot-spot quantities by propagating the normally distributed uncertainties in measured quantities through the 0D and 1D models described earlier. The model input parameters are those that fully describe the system, and are constrained by the measured yield, ion temperature, burn widths (from both X-rays and γ-rays), and volume from the 17% contour of neutron emissivity. Distributions of model parameters are generated using Markov chain Monte Carlo (MCMC), calculated with the tensorflow⁴⁹ probability package. The log-likelihood function for MCMC is defined by the measurements and calculated with the log-likelihood function

which is summed over all observables (i) where m_i is the model value, y_i is the measured value and δy_i is the uncertainty in the measurement. This methodology produces full distributions of the model parameters including any correlations, from the model parameter distributions we generate full distributions of all hot-spot parameters, some of which exhibit correlation, such as in the temperature and areal density required to evaluate the Hurricane metric, which are partially anti-correlated (evident in Fig. 3a). Other inferences, such as the implosion velocity or kinetic energy, are treated with normally distributed uncertainties that are uncorrelated with the hot-spot inferences.

Power-balance relations

In evaluating the power-balance relations relevant to equation (17) we use the following expressions for the individual terms:

In these expressions the specific powers are given in units of GJ g⁻¹ s⁻¹ and thus are multiplied by the inferred hot-spot mass to obtain power. ρ is the hot-spot mass in g cm⁻³, ⟨σν⟩ is the fusion reactivity evaluated as a function of temperature in cm³ s⁻¹, T is the temperature in keV, and ρR is the hot-spot areal density in g cm⁻². The self-heating power P_α is multiplied by the fraction of α-particle energy deposited in the hot spot (f_α) using the evaluation published in ref. ³⁶; for all four experiments, f_α ≈ 0.77–0.80.