Testing nonlocal models of electron thermal conduction for magnetic and inertial confinement fusion applications

Performing full integrated simulations of large-scale fusion devices, such as the National Ignition Facility (NIF) or the ITER tokamak, is very challenging due to the wide range of scales over which the important physical processes occur. Codes, such as HYDRA¹ and BOUT++,² used to perform integrated simulations of inertial and magnetic confinement fusion (ICF/MCF), respectively, often include reduced models to capture the complex aspects of the physics. The accuracy of the models used naturally affects the predictive capability of these codes. In this paper, we compare three different models for kinetic (i.e., nonlocal) effects on electron thermal conduction against Vlasov-Fokker-Planck (VFP) simulations: (i) the eigenvector integral closure (EIC)^3–5 and (ii) the non-Fourier Landau-fluid (NFLF)^6,7 models, which have recently been suggested for application in the tokamak edge and scrape-off layer (SOL); and (iii) the Schurtz, Nicolaï, and Busquet’s multigroup diffusion (SNB) model,^8–12 which is currently the most widely used in inertial fusion and laser-plasma applications.

In collisional plasmas, where the mean free path (mfp) is much less than the temperature scalelength, the electron heat flow parallel to any macroscopic magnetic field in the plasma has been shown by Braginskii¹³ to obey Fourier’s law

where $κ(B)$ is the dimensionless thermal conductivity, n_e the electron density, and $vT=kBTe/me$ is the thermal velocity

is an averaged electron-ion mean free path this can be ommitted in Gaussian units (which shall be used throughout this paper), k_B is Boltzmann’s constant, T_e is the electron temperature, Z is the average ionisation, e is the magnitude of the electron charge, and $log Λei$ is the Coulomb logarithm for electron-ion scattering which is typically between 2 and 10 in cases of interest here. Here and for the entirety of this paper, we assume there to be no drift velocity or current (hence, the ion and electron rest frames are equivalent). Note that Epperlein and Haines¹⁴ have calculated $κ(B)$ to an increased accuracy and Epperlein and Short¹⁵ later suggested that this can be approximated well by $κ(B)≈ξ(Z)128/3π$⁠, where $ξ(Z)=(Z+0.24)/(Z+4.2)$⁠.

Equation (1) breaks down if the collisionality of the electrons becomes low. This is due to the inadequacy of the assumption that the electron distribution function is close to Maxwellian; electrons with mfp’s larger than the temperature scalelength can in fact escape gradients before being scattered and depositing their energy into the plasma, leading to a distortion of the distribution function away from Maxwellian.

The largest contribution to the heat flow comes from suprathermal electrons with a velocity of approximately $4vT$⁠. Due to the v⁴ scaling of the appropriate mfp’s, these suprathermals can travel over a hundred times further than thermal electrons enabling excess heat to be deposited beyond the steepest part of the temperature profile (often referred to as “preheat” in the literature¹⁵) A reduced population of suprathermals is left behind in the region of steep temperature gradient, contributing to a reduction in the heat flux. These “nonlocal” effects can become important even for temperature scalelengths as long as $∼200$ thermal mfp’s.⁸ 

Such situations occur frequently in important regions of both MCF and ICF experiments: In tokamaks, nonlocal thermal transport is thought to play a role in heat flow from the core plasma to the “divertor,”¹⁶ a region of the tokamak edge specifically designed to absorb and exhaust excess heat from the plasma. Thermal electrons entering the SOL at the separatrix have mfp’s ranging from 1% (C-Mod) to 20% [DIII-D/Tokamak de Varennes (TdeV)] of the distance to the divertor target (connection length). For ITER, this is estimated to be 8%. In fact, the ratio of $λei(B)$ to the local temperature scalelength $LT=1/∇∥ log Te$ tends to vary along the SOL from approximately 1 (TdeV) or 0.1 (DIII-D) near the separatrix, to 0.01 near the colder divertor.¹⁷ These ratios are all two orders of magnitude higher for suprathermal electrons, rendering the heat transport inherently nonlocal. Furthermore, transient events—Edge Localised Modes (ELMs), disruptions and filaments—can create even higher temperatures and steeper gradients, with which the associated suprathermals would be almost collisionless.⁴ 

Current state of the art codes, such as SOLPS^18,19 and UEDGE,²⁰ have been shown to significantly underestimate the outer divertor target electron temperature and overestimate its density compared to the experiment in the existing tokamaks, which in turn affects other plasma parameters. Chankin and Coster²¹ have suggested that nonlocal effects in addition to inadequate treatment of neutrals (which has largely been ruled out by a sensitivity analysis) and inappropriate use of time-averaging could explain this discrepancy. The plausibility of this hypothesis is supported by recent gyrokinetic simulations performed by Churchill et al.²² Another important factor is the effect of the enhanced suprathermal population on Langmuir probe characteristics:^23–26 Ďuran et al.²⁷ have shown that a more sophisticated interpretation of probe results can reduce but not eliminate the discrepancy. Resolution of this discrepancy is critical as excessive heat loads could erode and severely limit the lifetimes of the divertor target plates.²⁸ 

For the case of indirect-drive inertial fusion, steep temperature gradients of approximately $100 μm$ are set up near the inner surface of the gold hohlraum that contains both the helium gas-fill and the fuel capsule. This is induced by the high-energy lasers which ionise and ablate the hohlraum wall. Ratios of $λei(B)/LT$ exceeding 10%–20% in this region have been reported.⁸ Significant nonlocal effects on the thermal conduction are consequently observed, particularly in the neighbouring low-density gas-fill where the mfp’s of heat-carrying electrons can be very long. Failure to simulate this nonlocality accurately can have implications for hohlraum temperatures and implosion symmetry.¹ 

A common approach to incorporate the flux reduction aspect of nonlocal transport is to restrict the local heat flow to some fraction $fL$ of the free-streaming limit $Qfs=nekBTevT$⁠. However, at best the flux-limiter $fL$ must be tuned against the previous experiments, limiting predictive capability—values ranging from 0.03 to 0.15 have been suggested for NIF design codes^1,29 and up to 3 for SOL modeling³⁰—and preheat effects cannot be predicted.

A more complete way to take nonlocal effects into account, however, is with a fully kinetic approach. By solving the Vlasov-Fokker-Planck (VFP) equation for the electron distribution function directly (along with self-consistent electric and magnetic fields), we need not assume it is close to Maxwellian; nonlocal effects are calculated ab-initio. Such an approach typically assumes binary collisions and small-angle scattering limiting its applicability in regions where the plasma is strongly coupled (⁠$log Λ$ approaching unity) such as in ICF fuel capsules or the colder part of the partially ionised hohlraum wall. While it is possible for VFP codes to treat collisions between multiple ion species³¹ and even neutrals³² (though the latter might require coupling to a neutral Monte Carlo code such as EIRENE^33–35 due to the importance of large-angle collisions), here we only consider the collisions of electrons with a single stationary ion species. Quantum-mechanical effects such as Fermi degeneracy, which could be of importance in solid density material, are also typically negelcted;³⁶ nevertheless, methods to incorporate these have been suggested.³⁷ 

However, due to the extra dimensionality associated with solving in velocity-space, VFP codes are computationally intensive and difficult to incorporate into the existing integrated modeling codes. Such demands of resolving the distribution function and collision time are especially restrictive in cold, dense regions such as deep in the hohlraum wall where a fluid treatment might even be sufficient. Therefore, alternative models that have more predictive capability than flux-limiters, and reduced computational requirements compared to a full kinetic simulation, are desirable. A dedicated experiment to measure nonlocal effects performed by Gregori et al.³⁸ has shown that a model of this kind can approximate measured temperature profiles better than flux-limiters.

A large number of reduced models for nonlocal electron thermal transport have been suggested for applications in inertial fusion and laser-plasmas^{8–12,15,39–42} and to the SOL.^{3–7,26,43,44} This paper focuses on three of these models, here referred to as the SNB,^8–12 EIC,^3–5 and NFLF^6,7 models (described in Sec. III), and compares them for accuracy against VFP simulations. While the SNB model has recently been compared to VFP results by Marocchino et al.,⁴⁵ this has shown that the two approaches agree well for Z = 1 but begin to deviate from each other as the ionisation increases. This is surprising as the SNB model was originally derived in the high-Z (Lorentz) limit, and we observe here that such findings are sensitive to precise implementation details of the model. Additionally, while the EIC and NFLF models have been shown to predict similar heat-fluxes,⁷ they have not yet been validated against a full time-dependent VFP code.

The first problem we investigate here is the damping of a small-amplitude sinusoidal temperature profile of various wavelengths in Sec. IV. This test will be used to justify a tuning parameter which can be applied to the SNB model to improve agreement with VFP simulations. We will additionally suggest a new analytic fit for the conductivity reduction and use this to obtain improved coefficients for the NFLF model.

In Sec. V A, we will then consider the case, more relevant to both hohlraums and the SOL, of a plasma with a large temperature variation. We will show that the SNB model agrees very well with VFP simulations using the same tuning factor as in the linearised problem described above and that the EIC and NFLF models overpredict the peak heat flux. While this suggests that the SNB model may also be useful in SOL simulations, we also consider potential improvements to the other models to improve their performance.

Finally, we will show in Sec. V B that the SNB model can significantly overpredict the heat flow relative to VFP in the low-density helium gas-fill using a problem particularly relevant to the ablated hohlraum wall. The importance of gradients in both average ionisation Z and electron density n_e here could mean the findings could also be important for the detached divertor scenario.

The evolution of the electron distribution function f_e, assuming small-angle scattering from binary collisions, can be described by the Vlasov-Fokker-Planck equation⁴⁶ 

where $v$ is the electron velocity, E and B are the electric and magnetic fields respectively, and m_e is the electron mass. Two of the three VFP codes used here, IMPACT and KIPP, employ a zero-current constraint, $∫fv d3v=0$ to calculate the electric field, which ensures quasineutrality. The third VFP code, SPRING, uses a more sophisticated approach which solves the Poisson and ion continuity equations with an implicit-moment method.^47,48

In this work, we consider only collisions of electrons with themselves and a single ion species using the standard Trubnikov-Rosenbluth^49,50 form of the Fokker-Planck collision operator $C=Cee+Cei$⁠, where

(applying standard Einstein summation over repeated indices). Here, we have defined

where $Zi=Z$ is the average ionisation and $Ze=−1$⁠, along with m_i the ion mass. The ion distribution function f_i is assumed by KIPP to be Maxwellian; here, we enforce the ion temperature to be equal to the electron temperature but this is not necessary,⁵¹ and all other codes and models assume a cold ion population and neglect terms of order $me/mi$⁠, simplifying the electron-ion collision operator to

where $δ(v)$ is the Dirac delta function and δ_ij is the Kronecker delta tensor.

For the case where variations only occur along magnetic field lines, symmetry in the perpendicular direction allows for elimination of the magnetic field by “gyro-averaging” (integrating azimuthally around the $v∥$ axis, this process is still valid even in the absence of magnetic fields); this yields the 1d2v (one-dimensional in space, two-dimensional in velocity) equation

where $〈·〉$ represents a gyro-average (an explicit representation of $〈C〉$ can be found in previous work by Xiong et al.⁵² and Chankin et al.³³) and $∥$ denotes components of vectors parallel to the magnetic field.

The KIPP code³³ is designed to solve this equation using Cartesian spatial and velocity grids. The code uses an operator splitting method suggested by Shoucri and Gagne^53,54 that treats the spatial derivative using a second-order explicit scheme followed by the electric field and collision operator terms using a first-order (in time, second-order in velocity) implicit scheme. The velocity grid covers the domain $v∥∈[−vmax,vmax],v⊥∈[0,vmax]$ where v_max is a user-defined parameter. The distribution function is simply taken to be zero at the exterior of the grid and reflected along the interior $v⊥=0$ axis.

A simplified approach is the diffusion approximation, which consists of expanding the distribution function in Cartesian tensors and truncating all but the first two terms (⁠$fe=f0(v)+v·f1(v)/v$⁠). This reduces the number of velocity-space dimensions to one thereby increasing efficiency and has been observed to correctly predict heat flows to within 5% for temperature scalelengths $LT≲10λei(B)$⁠.⁵⁵ The IMPACT code⁴⁶ (two-dimensional in space) employs this approach and makes a further simplification of ignoring angular scattering due to electron-electron collisions, valid in the Lorentz limit. In order to recover the correct local thermal conductivity for lower-Z plasmas, the factor $ξ(Z)$ is used in the electron-ion collision frequency. Our comparisons between IMPACT and KIPP suggest that these approximations do not greatly affect the results for the problems studied in Sec. V A. The equations solved by IMPACT, along with Ampere and Faraday’s Law, are thus

where

is the isotropic contribution of the electron-electron collision operator and

is the velocity-dependent electron-ion collision frequency.

IMPACT is fully implicit and first order in time, and uses fixed-point/Picard iterations to handle nonlinear terms. Note that the implicit treatment of the electric field enforces charge conservation at every iteration. The magnetic field and electron inertia (⁠$∂f1/∂t$⁠) terms have not been included in the simulations appearing in this paper. The main reason for using IMPACT in Sec. V B is that KIPP has not yet been extended to spatially varying ionisation along $s∥$⁠.

Finally, we also include the results previously obtained with the SPRING⁴⁷ VFP code, which takes a Cartesian tensor expansion to arbitrary order and does not neglect/approximate anisotropic electron-electron collisions. This code uses a linearised approach, i.e., the spatial gradient operator $∇$ is replaced by $ik$⁠, making it advantageous for analysing the small-amplitude sinusoidal temperature perturbations featured in Sec. IV, but not problems with large temperature perturbations.

The first model investigated here was originally proposed by Ji, Held, and Sovinec³ and is directly derived from simplifications of the VFP equation. Necessarily, the time-derivative term is neglected to allow the heat flow to be calculated based on input density and temperature profiles only and not the history of the distribution function; this assumption should be reasonable over “mean” SOL profiles (i.e., averaged over time to eliminate fine-scalelength fluctuations), but could break down for transient events with faster timescales such as filaments and ELMs.

The distribution function is expressed as $fe=f(0)+δf$⁠, where δf is a perturbation to an initial, Maxwellian, guess $f(0)$⁠. Assuming the perturbation δf is small, the nonlinear collision and electric field terms in the gyro-averaged VFP equation are linearised, which yields

where

is the linearised collision operator.

The next step is to attempt a separation of variables into $s∥$ and $v/vT$⁠, where $v$ is made up of parallel and perpendicular components $v∥$ and $v⊥$⁠, by expressing

where ψ_n are eigenfunctions of the operator $〈CL〉/v∥$⁠, which depends only on $v/vT$⁠, and λ_n the inverse of their eigenvalue with dimensions of length. Substituting (14) into (12) and assuming that the dependence of ψ_n on space through v_T is negligible (only valid when relative temperature perturbations are small globally) yields

By similarly decomposing the right-hand side into (orthogonal) eigenfunction contributions, a set of independent first-order ordinary differential equation’s (ODE’s) for A_n is formed that can be solved efficiently. Consequently, δf can be reconstructed and the nonlocal correction to the heat flux computed through an integral in $v∥$ (hence the nomenclature Eigenvector Integral Closure or EIC).

The advantage of this approach is that it circumvents the need to solve in velocity-space at every timestep (as a VFP code must). The main challenge is identifying a discrete description of the eigenfunctions ψ_n that converges rapidly for use in a numerical scheme. In practice, this is done by using an orthonormal polynomial moment-basis to express ψ_n as a vector and the operator $CL/v∥$ as a matrix. Ji et al.³ proposed a Legendre-Laguerre (LL) basis in pitch angle and total speed. This converges rapidly in the hydrodynamic limit but slowly in the collisionless limit. As an alternative, Omotani et al.⁵ proposed a Hermite-Laguerre (HL) basis, decoupling parallel and perpendicular velocity components, which allows for easier implementation of sheath boundary conditions.

While there are a lot of computational benefits to the EIC model over a full VFP code, a large number of eigenfunctions (at least 120 according to Omotani et al.⁵) are needed for convergence. The NFLF model^6,7 provides a cheaper approach, potentially solving as few as three second-order ODE’s, but without a direct link to the VFP equation.

The popular Landau-fluid approach initially proposed by Hammett and Perkins^43,56,57 provides a closure for the nonlocal heat flux $Q̃$ in Fourier space. This recovers the correct damping rate of a sinusoidal temperature perturbation in both the collisional and collisionless limits (where the wavelength is much longer/shorter than the thermal mfp). However, the Fourier transforms involved are inconvenient for complex SOL geometries and large temperature and density variations.

The innovation by Dimits, Joseph, and Umansky⁶ was to enable direct calculation of the nonlocal parallel heat flux in configuration space by approximating the closure as a sum of Lorentzians

where $Q̃(B)$ is the (parallel) Braginskii heat flow in reciprocal space, a parametrises the behaviour in the collisionless limit and is determined analytically, k is the wavenumber of the Fourier mode, N is the number of Lorentzians chosen by the user for the fit, and $αj,βj$ are fit parameters.

Equating the contribution from each Lorentzian to a dummy contribution $qj$⁠, rearranging and taking the inverse Fourier transform gives a set of N second-order ODE’s for each spatial direction of interest that can be used to recover the nonlocal heat flow

This approach also conveniently avoids the issue of defining the mean free path in reciprocal space.

The final model being investigated is the multigroup diffusion or “SNB” model named after the original authors Schurtz, Nicolaï, and Busquet.⁸ It is widely used in inertial fusion codes such as Lawrence Livermore National Laboratory’s HYDRA,¹ CELIA laboratory’s CHIC,⁵⁸ and the University of Rochester Laboratory for Laser Energetics’ (LLE) DRACO;¹² and it is applicable in multiple spatial dimensions.

The SNB model is best explained starting from the diffusion approximation of the kinetic equations [see Eqs. (8) and (9) above], along with neglecting time-derivatives for similar reasons as the EIC model. The isotropic part of the distribution function f₀ is then split into a Maxwell-Boltzmann distribution $f0(mb)=ne exp (−v2/2vT2)/(2πvT)3/2$ and a deviation $δf0=f0−f0(mb)$⁠. The anisotropic part $f1$ is similarly split, but the “Maxwellian” part $f1(mb)$⁠, obtained from substituting $f0(mb)$ into Eq. (9), is replaced by an alternative, $g1(mb)$

This modification achieves positive-definiteness without affecting the integral used to calculate the heat flow, and is argued to be compensated by other approximations of the model.⁸ Here, we have defined $λei*=ξλei=ξv/νei$⁠. Note that the factor of 3 difference from the original paper in $f1(mb)$ is simply due to the use of spherical harmonics by Schurtz et al. while we use a Cartesian tensor expansion.

Electric field terms in Eq. (8) are neglected and instead incorporated phenomenologically by defining a limited mfp

where the local form for $E=kBTe(∇ log ne+γ∇ log Te)$ is used, with $γ=1+3(Z+0.477)/2(Z+2.15)$⁠. Substituting $f1=g1(mb)+λei(E)∇δf0$ into equation, the stationary form of (8) obtains

This can be solved to compute the deviation from the local heat flow

The main computational advantage of this approach is through the use of efficient model collision operators that are local in velocity-space. This allows for a more effective discretisation into velocity/energy groups and removes the need to store the entire distribution function in memory. The original authors suggested using a velocity-dependent Krook (BGK) operator due to its simplicity, but Del Sorbo et al.¹⁰ have also employed a more realistic operator suggested by Albritton, Williams, Bernstein, and Swartz (AWBS)⁵⁹ 

where we have introduced the dimensionless number r to account for variation in SNB model implementations/description across publications: The original authors⁸ halved the geometrically averaged mfp $λe=Zλei*λei$ [see Eq. (23) of Schurtz et al.⁸ and also Sec. III C of Del Sorbo et al.¹⁰) which is equivalent to setting r = 4 (except for the treatment of electric field). However, in a later section of the original paper⁸ (III F) as well as Sec. II of a later publication,⁹ this technicality is not restated when demonstrating a link to the kinetic equations, from which a value of r = 1 could be interpreted.

Furthermore, our attempts to replicate previous comparisons between SNB and VFP⁴⁵ suggest that Marocchino et al. used r = 16. Using this value for r in the SNB model along with neglecting corrections to angular scattering from electron-electron collisions (i.e., ξ is set to one) happens to give a good agreement with VFP codes for Z = 1. However, this agreement is observed to get progressively worse as Z increases. In this paper, we show that using the BGK collision operator with a different value (r = 2) and $ξ=(Z+0.24)/(Z+4.2)$ gives a very good agreement of SNB with VFP across a wide range of problems (and values of Z).

Note that, despite the differential form of the AWBS operator, its use does not actually require a significant increase in computational time unless an attempt to parallelise over energy groups is being made. This is because the velocity-space first-order differential equation is simply closed from above with the boundary condition $δf0(v=∞)=0$⁠. In a finite difference scheme, this could simply be implemented by enforcing the highest energy group to be zero and solving for the next highest group first. (Bear in mind that discretisation in velocity-space need not be uniform.) However, we identify other issues with the AWBS operator in Sec. IV A that limit its usefulness and the SNB model using this operator is therefore not explored further.

First we consider the damping of a small-amplitude temperature perturbation $Te=T0+δT cos (kx)$ (often referred to as the Epperlein-Short¹⁵ test) with a constant uniform background density and ionisation. Due to nonlocal effects as the wavelength is reduced, the dimensionless thermal conductivity κ decreases from that predicted in the local limit, $κ(B)$⁠. In this section, we first detail the methodology used in setting up simulations of this problem and assessing the respective conductivity reductions before briefly commenting on the agreement between the EIC model and VFP results. Analysis of the long-wavelength limit will then be presented in Sec. IV A so as to motivate a suitable choice for the SNB model parameter r. Finally, a new fit function for the conductivity reduction as a function of $kλei(B)$ is derived in Sec. IV B by connecting the collisional and collisionless regimes, and is used to calculate fit coefficients for the NFLF model.

A sinusoidal perturbation with a relative amplitude of $δT/T0=10−3$ was initialised for the KIPP simulations. We used a uniform spatial grid of 127 cells over a half-wavelength with a non-uniform Cartesian velocity grid extending to $vmax=7vT$ (with parameters $mmax=256,EPS=1.01$ as defined by Chankin et al.³³). The two methods employed by Marocchino et al.⁴⁵ were used to calculate the conductivity reduction $κ/κ(B)$⁠: (1) directly from the peak heat flow divided by the predicted local heat flow (⁠$κ/κ(B)=Q̃/Q̃(B)$⁠) and (2) inferred from the exponential decay rate ρ of the temperature perturbation (⁠$κ/κ(B)=ρ/ρ(B)$⁠, where $ρ(B)=2κ(B)k2vTλei(B)/3$⁠).

The thermal conductivity obtained by both these methods fluctuated in time initially (due to initialising as a Maxwellian) and was tracked until both methods approached constant values. Due to incomplete convergence in timestep, these values were slightly different and an average was then taken between the two final conductivity reductions. In order to improve accuracy without using unnecessary amounts of computational time due to a tiny timestep (KIPP is only first-order accurate in time), extrapolation to zero timestep was performed by fitting a straight line of conductivity reduction against timestep. Such complications were unnecessary when using the inherently stationary models: Instead of evolving in time, it was possible to calculate heat flow (and thus conductivity reduction) from a single profile with a lower relative amplitude of $10−5$ for each wavelength.

Results obtained for thermal conductivity reduction $κ/κ(B)$ as a function of nonlocality parameter $kλei(B)$ are shown in Figs. 1 and 2 for an ionisation of Z = 1. The choice of two separate figures for the case of Z = 1 is to allow for clear identification of features at both high and low collisionality and to avoid an excessive number of comparisons on a single figure. Kinetic results from the linearised VFP code SPRING calculated by Epperlein⁴⁷ and provided numerically by Bychenkov et al.⁶⁰ are also shown in Fig. 2 for comparison against the nonlocal models.

Both the LL³ and the HL⁵ bases for the EIC model were investigated using 40,40 moments to achieve convergence to within 1% for $kλei(B)<1$⁠. Figure 1 shows that both bases agree well with KIPP (within 5% and 10%, respectively) in this limit. For higher $kλei(B)$ (see Fig. 2), the SPRING VFP results begin to deviate from both the EIC and KIPP results for a number of reasons: First, the onset of electron inertia effects at high $kλei(B)$⁠, ignored by the EIC model, introduces a phase shift between the heat flow and temperature perturbation in frequency space (i.e., the perturbation goes from being critically damped to possessing an oscillatory component) making evaluation of the decay rate difficult with KIPP (the linearised formulation of the SPRING code makes this easier, and likely more accurate; note that it is the modulus of the complex thermal conductivity that has been provided in this case).

Additionally, while the HL basis only requires two Laguerre modes in the collisionless limit due to the parallel-perpendicular decoupling, we found that even 160 HL modes were insufficent to achieve convergence to within 10% for $kλei(B)>2$⁠. The LL basis, however, manages to converge to within 1% for $kλei(B)<50$ using 20, 20 modes. The collisionless limit predicted by Chang and Callen⁶¹ is approached as the total number of LL modes is increased (see below and also Figs. 2 and 3 of Ji et al.³ whose results we have successfully replicated); this is about a factor of 1.8 less than the true collisionless heat flow predicted analytically and by the SPRING code (see Sec. IV B).

Results for an ionisation of Z = 8 are shown in Fig. 3. Here, 50, 50 moments in the LL basis were required to achieve convergence at high $kλei(B)$ with the EIC. The diffusion approximation made by IMPACT is shown to break down around $kλei(B)≈0.3$⁠. Note that the thermal conductivity reduction at which the SNB begins to deviate from kinetic results is about the same (⁠$κ/κ(B)≈0.2$⁠) for both Z = 1 and 8; the lower nonlocality parameter $kλei(B)$ at which this occurs is due to the reduced importance of electron-electron collisions at higher ionisations.

Bychenkov et al.⁶⁰ have shown that for long wavelength perturbations (i.e., in the hydrodynamic limit)

to second-order in $kλei(B)$⁠, where $b≈264$ in the Lorentz limit (⁠$Z=∞$⁠). As the assumptions of the EIC model are valid in this linear and collisional limit (except for neglection of electron inertia which would only introduce a time-dependent component into the heat flow), and convergence of the LL basis is relatively rapid (only 2 Legendre modes are theoretically needed), we have used it to calculate the value of b for various Z (while the KIPP prediction for Z = 1 was within 4% of the EIC, this was considered less accurate due to insufficient extension/resolution of the velocity grid). This was done by fitting a straight line on a graph of heat flow against $Zk2λei(B)2$ for $kλei(B)<10−3/Z$ (the lower range of $kλei(B)$ extended below $2×10−4/Z$ and there were typically at least six points on the graph).

Results using the EIC model are summarised in Table I and Fig. 4, which also includes numerical results¹⁴ and rational approximations^15,62 for the low-Z conductivity correction $κ(B)(Z)/κ(B)(∞)$⁠. We find that the approximation $b(Z)/b(∞)=Z/(Z+11/2)$ fits numerical results to within 7%, whereas simply using ξ overestimates b by as much as 43% at Z = 1. However, the implications of this for the validity of using ξ in IMPACT and the SNB model are not as serious as they seem because b only quantifies the initial deviation from the local limit, and the total heat flux is not very sensitive to marginal errors in b in the hydrodynamic limit.

Table II outlines the values of b predicted by the SNB model depending on the model collision operator and choice of source term. This has been derived in the low-amplitude and continuum-velocity limit as detailed in  Appendix A. Using the AWBS operator and the kinetic source term $∇·f1(mb)$ on the right-hand side of Eq. (16) gives a priori the closest value of $b=316.9ξ$ (top right) to within 20% of that predicted analytically in the Lorentz limit (Table I).

The ability of the AWBS collision operator to predict the deviation in the hydrodynamic limit fairly accurately might suggest that it provides an improvement to the original SNB model; however, we find that coupling it with the original source term leads to negative values of the thermal conductivities at $kλei(B)>0.124/ξZ$ due to it not being positive-definite (see  Appendix B). This should never occur in the linearised problem considered here (i.e., decay of a small-amplitude temperature perturbation) as it would result in instabilities at these wavelengths. However, this issue does not necessarily imply that the AWBS operator is an inappropriate choice for other nonlocal models. For example, the M1 model presented by Del Sorbo et al.^10,11 does not appear to exhibit this issue of positive-definitiveness; we leave a detailed analysis of this model for future work.

Setting r = 2 exactly in the original implementation of the SNB model (BGK collision operator with the source term $∇·g1(mb)$⁠) remarkably gives the same value of $b=316.9ξ$ as with the AWBS operator and the source term $∇·f1(mb)$ (compare bottom left and top right entries of Table II) and in fact the entire distribution function in this limit (see  Appendix A). However, to match the kinetic results for b, a value of r = 2.4 is required in the Lorentz limit and r = 3 for Z = 1. We suggest that matching coefficients to such accuracy is not necessary and that using r = 2 achieves much better agreement for problems involving large temperature variations (see below). Results using both r = 2 and r = 3 for Z = 1 have been provided in figures to enable the reader to compare.

Faithfulness to kinetic results for b can be guaranteed with the NFLF model by modifying the analytical Fourier closure and constraining the fit coefficients appropriately as described in Sec. IV B.

With decreasing wavelength, the heat flow is predicted to slowly approach a constant value. By fitting to the results of both the EIC and SPRING models (we were unable to obtain meaningful KIPP results at low enough collisionalities due to issues mentioned above), we find that

where $Q̃fs$ is antiparallel to the wave-vector and χ₁, $c∞$⁠, and ϵ are dimensionless fit parameters, is a reasonable description for the asymptotic behaviour in this limit for low-Z plasmas [i.e., a graph of $Q̃$ against $log (k)$ resembles a straight line]. The form of this fit function was taken from work by Bychenkov et al.⁶⁰ The extra factor $3/2$ compared to the formalism of Hammett and Perkins⁴³ (which inspired previous implementations of the NFLF model^6,7) was found to be necessary due to the isotropic definition of the electron temperature used here.^47,60

Again, the LL basis was used, this time due to the nonconvergence of the HL basis explained above; however, at least 40, 40 moments were needed to achieve convergence above $kλei(B)≈1$⁠. As found by the original developers of the model,³ the value of $χ1=1.2/π$ agreed with the value predicted by Chang and Callen;⁶¹ this is exactly 40% less than the value predicted by Hammett and Perkins⁴³ (⁠$χ1=2/π$⁠) because of the quasi-stationary assumption. We have also calculated that an alternative value of $χ1=1.225$ can be obtained by numerically solving for zeroes of the response function.

Calculated values of ϵ and $c∞$⁠, as well as the c₁ referred to below, are summarised in Table III for $Z=1,2,4,6,8$⁠. Simulations with EIC at higher Z require a prohibitive number of moments for convergence at high $kλei(B)$⁠. Both the index ϵ and the coefficient $c∞$ vary weakly with Z and have similar orders of magnitude to those predicted by Bychenkov et al.⁶⁰ The values obtained here should be slightly closer to reality as Bychenkov et al. assume complete stationarity (all time derivatives neglected) in their calculations, but there are large uncertainties in our numerical fit (approximately 10% for the EIC data). The limited numerical results available from the assumingly exact SPRING code^47,60 infer a value for ϵ at Z = 1 within 0.5% of the EIC prediction, but the value for $c∞$ (=1.36) is larger by a factor of 2.2.

Due to the combination of stationarity and diffusion approximations, the SNB model without the phenomenological mfp limitation to include electric fields predicts the collisionless heat flow to decrease as $∼1/k$ to zero as the wavelength decreases⁸ (the thermal conductivity correspondingly decreases as $1/k2$⁠). Incorporating the mfp limitation resolves the issue of insufficiently damping temperature perturbations of finite amplitude (such that $kλeiδT≳1$⁠). This improves numerical stability, but introduces an amplitude-dependence of χ₁ that is not observed in VFP simulations.

While the NFLF will also always predict a $∼1/k$ decay of the heat flow for high enough $kλei(B)$⁠, increasing the number of Lorentzians used to improve the fit can progressively extend the validity into lower collisionality regimes. The fitting function we used interpolates behaviour in both the hydrodynamic and collisionless limits with a similar but slightly more robust method than used by Bychenkov et al.⁶⁰ 

where c₁ differs from $c∞$ by optimising the fit for $kλei(B)≤1$⁠. Using the parameters as defined in Table III for Z = 1, Eq. (25) fits the KIPP and SPRING results to within 6% and 10%, respectively, for $kλei≤1$ and up to $26/20%$ above this; altering the value of c₁ to 1.5 reduces the maximum discrepancy with SPRING results to 11%.

This new fit is depicted in Fig. 1 with the simpler fit $1/(1+akλei(B))$ obtained by Hammett and Perkins⁴³ previously used in the NFLF model,⁶ (a can be related to χ₁ by $a=2κ(B)/32χ1$⁠), which overestimates the thermal conductivity at moderate collisionalities around $kλei(B)≈0.5$ by over 25%. Note that we have observed a recent closure in configuration space (thus convenient for convolution models) suggested by Ji and Held⁶³ to closer fit the EIC results with one more fitting parameter (if the α used by Ji and Held is not considered a free parameter)—tuning of these parameters could probably also achieve an improved fit to the VFP results. We would also like to highlight a recent paper by Joseph and Dimits who have performed a detailed analysis of closures for the response function that connects the collisionless and collisional regime.⁶⁴ 

Three Lorentzians [i.e., N = 3 in Eq. (16)] can approximate our new fit to within 2.5% up to $kλei(B)≈1.6$⁠; using six Lorentzians allows this to be extended up to $kλei(B)≈30$⁠. The coefficients used are given in Table IV and were obtained using the variable projection method,⁶⁵ constrained such that Eq. (23) is obeyed to second-order in $kλei(B)$⁠.

We now investigate the accuracy of the EIC, NFLF, and SNB models in calculating the heat flow in the case where we have a large relative temperature variation. We consider the case of an initial temperature profile consisting of a ramp connecting two large hot and cold regions (“baths”). This has the advantages of allowing simple reflective boundary conditions and not requiring any external heating/cooling mechanisms that would also need to be carefully calibrated between codes. Results and initial conditions are here presented in terms of reference quantities encouraging the translation of the problem to both ICF and MCF relevant situations.

The hot and cold baths had temperatures of T₀ and $0.15T0$⁠; these were connected by a cubic ramp given by

where $nc′∈[−154,100]$ is the cell number counting from the centre of the temperature ramp. Cell size in mfp’s was varied between simulations to scan a range of collisionalities. The initial temperature profile is illustrated in Fig. 5 for the smallest cell-size used. For these simulations the electron density, Coulomb logarithm and ionisation (Z = 1) were all kept constant and uniform.

KIPP simulations showed an evolution of the heat flow from the local (due to initialising as a Maxwellian) to the nonlocal, with a reduced peak, over an initial transient phase (over which the temperature ramp flattened somewhat). The transient phase was considered over when the ratio of the KIPP heat flow to the expected local heat flow stopped decreasing. After the transient phase, this ratio begins to slowly increase as the thermal conduction flattens the temperature ramp and the ratio of the scalelength to mfp increases (i.e., the thermal transport slowly becomes more local). We then took the temperature profile from KIPP and used our implementation of the various nonlocal models to calculate the heat flow they would predict given this profile.

Figure 5 shows that the EIC and NFLF models agree well with each other (to within 10% almost everywhere for the snapshot shown). However, agreement with KIPP is not nearly as good; the models overestimate the peak heat flux by 30%–35% and do not predict the observed preheat into the cold region. The SNB model is shown to perform much better here, predicting the peak heat flux to within 6% as well as the existence of preheat (although this is overestimated).

The wide range of heat flow profiles predicted with different flux-limiters between 0.05 and 0.7 are also shown in Fig. 5. These were obtained using the formula $1/Qfl=1/Q(B)+1/fLQfs$⁠. We find that a flux-limiter of $∼0.25$ best matches the peak kinetic heat flow, but in this case the peak is shifted towards the hot rather than the cold bath (the latter is observed in the KIPP simulation). Similar results are observed at all temperature ramp scalelengths investigated as illustrated in Fig. 6, which depicts the reduction in the peak heat flow compared to the local Braginskii prediction.

While comparisons between the SNB model and VFP codes have previously been performed,^8,45 none have included spatially inhomogeneous ionisation. As inertial fusion experiments involve steep ionisation and density gradients (for example, at the interface between the helium gas-fill and the ablated gold plasma), it is critical that the SNB model be tested in such an environment. Variations in ionisation may also be important in the “detached” divertor scenario where a moderate-Z gas is injected in front of the divertor to radiate excess heat; an investigation of this scenario is left as further work. For evaluating this, the IMPACT⁴⁶ VFP code was used due to its ability to simulate inhomogeneous ionisation profiles.

We performed a HYDRA simulation in 1D spherical geometry of a laser-heated gadolinium hohlraum containing a typical helium gas-fill. A leak source was implemented with an area equal to the laser entrance hole to reproduce the energy balance. Electron temperature (T_e), density (n_e), and ionisation (Z) profiles (shown in Fig. 7) at 20 ns were extracted and used as the initial conditions (along with the assumption that the electron distribution function is initially Maxwellian everywhere) for the IMPACT simulation (which was performed instead in planar geometry). At this point, very steep gradients in all three variables were set up with a change from $Te=2.5 keV,ne=5×1020 cm−3,Z=2$ to $Te=0.3 keV,ne=5×1021 cm−3,Z=39$ across approximately 100 μm at the helium-gadolinium interface.

Spline interpolation was used to increase the spatial resolution near the steep interface for the IMPACT simulations, helping both numerical stability and runtime due to needing a reduced number of nonlinear iterations. For simplicity, the Coulomb logarithm was treated as a constant $log Λei=log Λee=2.1484$⁠. Note that in reality the plasma is only strongly coupled in the colder region of the gadolinium bubble beyond $∼$ 1.7 mm and $log Λei≈8$ up to $∼$ 1.6 mm in the hotter corona. Reflective boundary conditions were used here as in Sec. V A and IMPACT used a timestep of 1.334 fs. The n_e and Z profiles did not evolve in the IMPACT simulation due to the treatment of the electric field discussed in Sec. II that ensures quasineutrality and the neglection of ion motion and ionisation models.

As with the KIPP simulations in Sec. V A, there is an initial transient phase where the IMPACT heat flux gradually reduces from the Braginskii prediction as the distribution function rapidly moves away from Maxwellian. Once again, this transient phase is considered to be over when the ratio of the peak heat flow to the Braginskii prediction stops reducing. This ratio is not observed to change by more than 5% after the first 5 ps of our 15.7 ps simulation. Therefore, we conclude that it is safe to assume the transient phase is over after 5 ps, at which point the temperature front has advanced by approximately 8 $μ$m leading to a maximum temperature change of 41% as shown in Fig. 7.

In comparing the IMPACT and SNB heat flow profiles, we encountered an important subtlety concerning the implementation of the model. While more recent publications concerning the SNB model^9,10 use a formulation similar to that used here in Sec. III C with separate electron-ion and electron-electron mfp’s or collision frequencies, the original paper⁸ used a geometrically averaged mfp $λe=Zλei*λei$⁠. However, this averaging process is only valid for the case of homogeneous ionisation, and Fig. 8 shows the large effect this has on the heat flow when the ionisation varies. While using separated mfp’s provides a very good prediction of the preheat into the hohlraum, the peak heat flow to within 16% and an improved estimate of the thermal conduction in the gas-fill region, the latter is still too large by a factor of $∼2$⁠. This discrepancy could potentially lead to an overestimate of hohlraum temperatures and thus causing the issues similar to those arising with using an overly restrictive flux limiter.¹ 

The capability of the NFLF to closely match the results of EIC for the case of homogeneous density and ionisation is fairly impressive, considering that only 6 Lorentzians were needed for convergence compared to EIC’s 64 moments (16, 4 in the HL basis, chosen instead of LL as convergence is faster for this problem). This implies that the NFLF is about 5 times faster (assuming the NFLF’s second-order ODE’s take approximately twice the time to solve as EIC’s first-order). However, this result should not be too surprising as both models are based on some kind of linearisation procedure, causing them to fail in almost exactly the same way for a nonlinear problem. For example, the lack of preheat or spatial shift in peak location predicted by the models are both features observed in the linear problem studied in Sec. IV. The SNB model requires 25 groups for convergence resulting in an only slightly faster computation time than the EIC model.

Improving performance of the models for large temperature variations would require approaches that did not affect the desirable agreement in the linearised limit. For EIC, a simple method is nonlinear iteration; i.e., updating the right-hand side of Eq. (15) by adding on nonlinear terms such as $eE∥me∂δf∂v∥−∑nAn∂ψn∂s∥$ from the initial calculation and repeating until convergence. However, the computational time to apply the differential operators and separate into eigenvector components would probably increase the computational time by an undesirably large factor on the order of the number of moments used.

Conversely, a correct approach for improving the NFLF model is not immediately apparent and probably requires deeper analysis of the link between the model and the VFP equation. However, it is conceivable that this could be done without additional computational expense; for example, replacing the $a2λei2∇2$ term in Eq. (17) with $a2(λei∇)·(λei∇)$ would affect only nonlinear behaviour.

It is important to investigate the sensitivity of divertor temperature to the errors in these models to confirm whether an accurate treatment of nonlocal transport can reconcile simulation and experiment. Furthermore, the discrepancies observed with the SNB model when ionisation gradients are steep could potentially have critical knock-on effects for integrated ICF modeling; it needs to be determined whether further improvements to the SNB model are necessary to avoid this.

One key neglection in this work is the effect of electron-electron collisions on the anisotropic part of the distribution function $f1$ for the case of spatially varying ionisation. It was shown in Sec. IV A that inclusion of this in KIPP and EIC predicts a noticeably different nonlocal deviation (consider, for example, the value b) than would be predicted by using the phenomenological collision fix ξ [which incorrectly predicts $b(Z)/b(∞)=κ(B)(Z)/κ(B)(∞)=ξ$ as depicted in Fig. 4]. But this did not seem to be the case for the more physically realistic large temperature variation studied in Sec. V A, as using the value r = 2 in the SNB model, derived in the linearised and Lorentz limits, seemed to be preferable to r = 3. Nevertheless, the use of ξ in IMPACT as an ad-hoc substitution for a more complete approach to anisotropic electron-electron collisions could still potentially lead to inaccuracies in the heat flow predictions depicted in Sec. V B, and this should be further investigated.

Less critical to our findings are the inaccuracies experienced by VFP codes in strongly coupled plasmas. While this could play a role in the cooler part of the hohlraum wall studied in Sec. V B where the Coulomb logarithm drops to $∼2$ (theoretically rendering the effect of collisions in this region only accurate to $∼50%$⁠), it does not affect the conclusion that the separated SNB model predicts the same heat flow into the wall as IMPACT while overpredicting that in the corona as both use the same treatment of $log Λ$⁠. We have simply shown quantitatively that reduced models can be an effective stepping stone between hydrodynamic and VFP approaches. However, this does act as a reminder that even a highly sophisticated VFP code could be faced with challenging inaccuracies in certain regions of the plasma (though it would surely still be an improvement to a purely hydrodynamic approach which would experience the same difficulties with strongly coupled plasmas); a potential method in overcoming this and incorporating large-angle collisions in a continuum code could be a Monte Carlo based approach.⁶⁶ Similar points can be made for other deficiencies, such as collisions with neutrals and Fermi degeneracy, although these are probably slightly easier to address and incorporate into models.^32,37

Following on from these basic test problems and sensitivity tests, there are still important questions on predictive modeling of fusion plasma heat flows that could be answered using VFP codes. First, the distribution function predicted by the SNB model should be compared to that of a fully kinetic code to assess the former’s viability in predicting other transport coefficients or parameteric instabilities.⁶⁷ Further modifications of the distribution function to a Dum-Langdon-Matte type super-Gaussian^68–70 due to inverse bremsstrahlung by laser heating in inertial fusion could also significantly alter the transport processes.⁷¹ Furthermore, kinetic effects can still affect perpendicular transport (both heat flow and magnetic field advection rates) for moderate magnetisations;^72,73 this could be relevant to the recent interest in magnetised hohlraums^74,75 or magnetic islands in tokamaks;⁷⁶ and while a few reduced models have been suggested to capture some of these aspects^9,11 they still need to be properly validated with kinetic codes.

In conclusion, we have compared three nonlocal models from ICF and MCF. We have demonstrated their optimal implementations, revealing potential subtleties in the description of the models. We have demonstrated that the SNB model—using the original BGK operator, but scaled according to an analysis of small-amplitude temperature sinusoids (r = 2), along with the modified source term $∇·g1(mb)$ appearing on the right-hand side of Eq. (20)—performs better than NFLF and EIC for the problems investigated with large temperature variations. Ensuring that the electron-electron and electron-ion collisionalities appear separately in this equation further improves agreement with VFP for a problem with spatially varying ionisation. However, the problems studied with large temperature variation only reach a nonlocality parameter of $∼15%$⁠, suggesting that SNB is most likely suitable for modeling hohlraum energetics problems (with the current exception of gas-fill heat flow, which is overestimated by a factor of $∼2$⁠) and mean SOL profiles but could break down at the even shorter scalelengths relevant to transient events.

The NFLF and EIC models have been found to perform favourably against KIPP when predicting the rate of decay of a small-amplitude temperature perturbation over a wider range of collisionalities than the SNB. However, these models overestimate the peak heat flux by up to 35% in the case of a large temperature variation, as well as failing to predict preheat. Additionally, a new analytic fit to kinetic results for temperature sinusoids has been presented in Eq. (25) that could be useful in traditional Landau-fluid implementations.

The authors would like to thank J. Y. Ji for sharing the numerical results of his closure,⁶³ and A. V. Brantov for assistance in understanding his work with Bychenkov et al.^42,77 We would also like to express our gratitude towards A. M. Dimits, I. Joseph, W. Rozmus, and L. LoDestro for interesting discussions and pointing us towards a number of relevant references. We additionally appreciate valuable feedback on our manuscript from O. Izacard and our Physics of Plasmas reviewer. Finally, thanks to M. Zhao for sharing a script to analyse the KIPP results.

All data used to produce the figures in this work, along with other selected supporting data, can be found at http://dx.doi.org/10.15124/3de4e00c-9087-4d95-89aa-9acd2071d3fb.

This work is funded by EPSRC Grant Nos. EP/K504178/1 and EP/M011372/1. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under Grant Agreement No. 633053 (Project Reference CfP-AWP17-IFE-CCFE-01). The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. This work was supported by Royal Society award IE140365.

For long wavelength perturbations, the diffusion term in Eq. (20) can be ignored, and thus, the distribution function and nonlocal heat flow easily computed in this limit. An outline of the derivation is given here; note that using the BGK collision operator and $g1(mb)$ in the source term gives the same $δf0$ as when using the AWBS operator with $f1(mb)$ in the source term if r = 2. The different consequences of choosing each source term $g1(mb), f1(mb)$ are distinguished by the terms, on the left and right respectively, inside the curly brackets ${·,·}$⁠. Note that integration by parts is employed for the AWBS calculation and a change of variables to $u=v/2vT$ is used. Additionally, we define and use a further two dimensionless variables $X=ξZk2λei(v)2$ which is velocity-dependent and $X(B)=ξZk2λei(B)2$ which is independent of velocity. Numerical results of these calculations are summarised in Table II.

A similar analysis can be performed with slightly greater difficulty at arbitrary collisionality. Integration by parts must be used again for the AWBS derivation, along with some mathematical identities. Recall that the electric field correction made by the SNB model is a nonlinear correction and does not come into play if the amplitude of the perturbation is infinitesimal

where γ is the incomplete gamma function. Computing the definite integral numerically with Mathematica shows that the AWBS heat flow can become negative for $X(B)>0.0154$⁠, which corresponds to $kλei(B)>0.124/ξZ$⁠.