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

Three models for nonlocal electron thermal transport are here compared against Vlasov-Fokker-Planck (VFP) codes to assess their accuracy in situations relevant to both inertial fusion hohlraums and tokamak scrape-off layers. The models tested are (i) a moment-based approach using an eigenvector integral closure (EIC) originally developed by Ji, Held and Sovinec; (ii) the non-Fourier Landau-fluid (NFLF) model of Dimits, Joseph and Umansky; and (iii) Schurtz, Nicola\"i and Busquet's multigroup diffusion model (SNB). We find that while the EIC and NFLF models accurately predict the damping rate of a small-amplitude temperature perturbation (within 10% at moderate collisionalities), they overestimate the peak heat flow by as much as 35% and do not predict preheat in the more relevant case where there is a large temperature difference. The SNB model, however, agrees better with VFP results for the latter problem if care is taken with the definition of the mean free path. Additionally, we present for the first time a comparison of the SNB model against a VFP code for a hohlraum-relevant problem with inhomogeneous ionisation and show that the model overestimates the heat flow in the helium gas-fill by a factor of ~2 despite predicting the peak heat flux to within 16%.


I. INTRODUCTION
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 1 and BOUTþþ, 2 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][4][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][9][10][11][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 13 to obey Fourier's law Q ðBÞ ¼Àj ðBÞ n e v T k ðBÞ ei rk B T e ; (1) where j ðBÞ is the dimensionless thermal conductivity, n e the electron density, and v T ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T e =m e p 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 K 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 14 have calculated j ðBÞ to an increased accuracy and Epperlein and Short 15 later suggested that this can be approximated well by j ðBÞ % nðZÞ128=3p,wherenð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 4v T . Due to the v 4 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 15 ) 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. 8 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," 16 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 k ðBÞ ei to the local temperature scalelength L T ¼ 1=r k log T e 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. 17 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. 4 Current state of the art codes, such as SOLPS 18,19 and UEDGE, 20 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 21 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. 22 Another important factor is the effect of the enhanced suprathermal population on Langmuir probe characteristics: [23][24][25][26] Ď uran et al. 27 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. 28 For the case of indirect-drive inertial fusion, steep temperature gradients of approximately 100 lm 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 k ðBÞ ei =L T exceeding 10%-20% in this region have been reported. 8 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. 1 A common approach to incorporate the flux reduction aspect of nonlocal transport is to restrict the local heat flow to some fraction f L of the free-streaming limit Q fs ¼ n e k B T e v T . However, at best the flux-limiter f L 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 andupto3forSOLmodeling 30 -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 K 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 31 and even neutrals 32 (though the latter might require coupling to a neutral Monte Carlo code such as EIRENE [33][34][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; 36 nevertheless, methods to incorporate these have been suggested. 37 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. 38 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][9][10][11][12]15,[39][40][41][42] and to the SOL. [3][4][5][6][7]26,43,44 This paper focuses on three of these models, here referred to as the SNB, [8][9][10][11][12] EIC, [3][4][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., 45 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, 7 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. VA, 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. VB 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.

II. VLASOV-FOKKER-PLANCK MODELING
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 46 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, Ð f v d 3 v ¼ 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 colli- (applying standard Einstein summation over repeated indices). Here, we have defined where Z i ¼ Z is the average ionisation and Z e ¼À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, 51 and all other codes and models assume a cold ion population and neglect terms of order m e =m i , simplifying the electron-ion collision operator to (6) where dðvÞ is the Dirac delta function and d 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 k 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 hÁi represents a gyro-average (an explicit representation of hCi can be found in previous work by Xiong et al. 52 and Chankin et al. 33 ) and k denotes components of vectors parallel to the magnetic field. The KIPP code 33 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 k 2½Àv max ; v max ; v ? 2½0; v max 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 (f e ¼ f 0 ðvÞþv Á f 1 ð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 L T Շ10k ðBÞ ei . 55 The IMPACT code 46 (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 nð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. VA 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 (@f 1 =@t) terms have not been included in the simulations appearing in this paper. The main reason for using IMPACT in Sec. VB is that KIPP has not yet been extended to spatially varying ionisation along s k .
Finally, we also include the results previously obtained with the SPRING 47 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 r 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.

A. Eigenvector integral closure
The first model investigated here was originally proposed by Ji, Held, and Sovinec 3 and is directly derived from simplifications of the VFP equation. Necessarily, the timederivative 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 f e ¼ f ð0Þ þ df , where df is a perturbation to an initial, Maxwellian, guess f ð0Þ . Assuming the perturbation df is small, the nonlinear collision and electric field terms in the gyro-averaged VFP equation are linearised, which yields where C L ðdf Þ¼C ee ðf ð0Þ ; df ÞþC ee ðdf ; f ð0Þ ÞþC ei ðdf ; n i dðvÞÞ; is the linearised collision operator. The next step is to attempt a separation of variables into s k and v=v T , where v is made up of parallel and perpendicular components v k and v ? , by expressing where w n are eigenfunctions of the operator hC L i=v k , which depends only on v=v T , and k n the inverse of their eigenvalue with dimensions of length. Substituting (14) into (12) and assuming that the dependence of w n on space through v T is negligible (only valid when relative temperature perturbations are small globally) yields (15) 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, df can be reconstructed and the nonlocal correction to the heat flux computed through an integral in v k (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 w n that converges rapidly for use in a numerical scheme. In practice, this is done by using an orthonormal polynomial moment-basis to express w n as a vector and the operator C L =v k as a matrix. Ji et al. 3 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. 5 proposed a Hermite-Laguerre (HL) basis, decoupling parallel and perpendicular velocity components, which allows for easier implementation of sheath boundary conditions.

B. Non-Fourier Landau-fluid
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. 5 ) 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 fluxQ in Fourier space. This recovers the correct damping rate of a sinusoidal temperature perturbation in both the collisional and collisionless limits (where the 092309- 4 Brodrick et al. Phys. Plasmas 24, 092309 (2017) 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 6 was to enable direct calculation of the nonlocal parallel heat flux in configuration space by approximating the closure as a sum of Lorentzians whereQ ð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 a j ; b j are fit parameters. Equating the contribution from each Lorentzian to a dummy contribution q j , 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.

C. Multigroup diffusion (SNB)
The final model being investigated is the multigroup diffusion or "SNB" model named after the original authors Schurtz, Nicola€ ı, and Busquet. 8 It is widely used in inertial fusion codes such as Lawrence Livermore National Laboratory's HYDRA, 1 CELIA laboratory's CHIC, 58 and the University of Rochester Laboratory for Laser Energetics' (LLE) DRACO; 12 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 0 is then split into a Maxwell- . The anisotropic part f 1 is similarly split, but the "Maxwellian" part f 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. 8 Here, we have defined k Ã ei ¼ nk ei ¼ nv= ei . Note that the factor of 3 difference from the original paper in f where the local form for 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. 10 have also employed a more realistic operator suggested by Albritton, Williams, Bernstein, and Swartz (AWBS) 59  8 and also Sec. III C of Del Sorbo et al. 10 ) which is equivalent to setting r ¼ 4 (except for the treatment of electric field). However, in a later section of the original paper 8 (III F) as well as Sec. II of a later publication, 9 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 45 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 electronelectron collisions (i.e., n 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 n ¼ð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 092309- 5 Brodrick et al.
Phys. Plasmas 24, 092309 (2017) 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 df 0 ðv ¼1Þ ¼ 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.

IV. DECAY OF A SMALL-AMPLITUDE, SINUSOIDAL TEMPERATURE PERTURBATION
First we consider the damping of a small-amplitude temperature perturbation T e ¼ T 0 þ dT cos ðkxÞ (often referred to as the Epperlein-Short 15 test) with a constant uniform background density and ionisation. Due to nonlocal effects as the wavelength is reduced, the dimensionless thermal conductivity j decreases from that predicted in the local limit, j ðBÞ .I n 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 kk ðBÞ ei 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 dT=T 0 ¼ 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 v max ¼ 7v T (with parameters mmax ¼ 256; EPS ¼ 1:01 as defined by Chankin et al. 33 ). The two methods employed by Marocchino et al. 45 were used to calculate the conductivity reduction j=j ðBÞ : (1) directly from the peak heat flow divided by the predicted local heat flow (j=j ðBÞ ¼Q=Q ðBÞ ) and (2) inferred from the exponential decay rate q of the temperature perturbation (j=j ðBÞ ¼ q=q ðBÞ ,whereq ðBÞ ¼ 2j ðBÞ k 2 v T k ðBÞ ei =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 firstorder 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 j=j ðBÞ as a function of nonlocality parameter kk  60 are also shown in Fig. 2 for comparison against the nonlocal models.
Both the LL 3 and the HL 5 bases for the EIC model were investigated using 40,40 moments to achieve convergence to within 1% for kk ðBÞ ei < 1. Figure 1 shows that both bases agree well with KIPP (within 5% and 10%, respectively) in this limit. For higher kk ðBÞ ei (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 kk ðBÞ ei , 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 parallelperpendicular decoupling, we found that even 160 HL modes were insufficent to achieve convergence to within 10% for kk ðBÞ ei > 2. The LL basis, however, manages to converge to within 1% for kk ðBÞ ei < 50 using 20, 20 modes. The collisionless limit predicted by Chang and Callen 61 is approached as the total number of LL modes is increased (see below and also Figs. 2 and 3 of Ji et al. 3 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 kk ðBÞ ei with the EIC. The diffusion approximation made by IMPACT is shown to break down around kk ðBÞ ei % 0:3. Note that the thermal conductivity reduction at which the SNB begins to deviate from kinetic results is about the same (j=j ðBÞ % 0:2) for both Z ¼ 1 and 8; the lower nonlocality parameter kk ðBÞ ei at which this occurs is due to the reduced importance of electron-electron collisions at higher ionisations. Bychenkov et al. 60 have shown that for long wavelength perturbations (i.e., in the hydrodynamic limit) to second-order in kk ðBÞ ei , where b % 264 in the Lorentz limit (Z ¼1). 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 timedependent 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 Zk 2 k Results using the EIC model are summarised in Table I and Fig. 4, which also includes numerical results 14 and rational approximations 15,62 for the low-Z conductivity correction j ðBÞ ðZÞ=j ðBÞ ð1Þ. We find that the approximation bðZÞ= bð1Þ ¼ Z=ðZ þ 11=2Þ fits numerical results to within 7%, whereas simply using n overestimates b by as much as 43% at Z ¼ 1. However, the implications of this for the validity of using n 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.U s i n gt h e AWBS operator and the kinetic source term rÁf ðmbÞ 1 on the right-hand side of Eq. (16) gives ap r i o r ithe closest value of b ¼ 316:9n (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   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 rÁg  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 whereQ fs is antiparallel to the wave-vector and v 1 , c 1 , and are dimensionless fit parameters, is a reasonable description for the asymptotic behaviour in this limit for low-Z plasmas [i.e., a graph ofQ against log ðkÞ resembles a straight line]. The form of this fit function was taken from work by Bychenkov et al. 60 The extra factor 3=2 compared to the formalism of Hammett and Perkins 43 (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 kk ðBÞ ei % 1. As found by the original developers of the model, 3 the value of v 1 ¼ 1:2= ffiffiffi p p agreed with the value predicted by Chang and Callen; 61 this is exactly 40% less than the value predicted by Hammett and Perkins 43 (v 1 ¼ 2= ffiffiffi p p ) because of the quasi-stationary assumption. We have also calculated that an alternative value of v 1 ¼ 1:225 can be obtained by numerically solving for zeroes of the response function.
Calculated values of and c 1 , as well as the c 1 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 kk ðBÞ ei . Both the index and the coefficient c 1 vary weakly with Z and have similar orders of magnitude to those predicted by Bychenkov et al. 60 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 (¼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 8 (the thermal conductivity correspondingly decreases as 1=k 2 ). Incorporating the mfp limitation resolves the issue of insufficiently damping temperature perturbations of finite amplitude (such that kk ei dT տ 1). This improves numerical stability, but introduces an amplitudedependence of v 1 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 kk ðBÞ ei , 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. 60 where c 1 differs from c 1 by optimising the fit for kk ðBÞ ei 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 kk ei 1 and up to 26=20% above this; altering the value of c 1 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 þ akk ðBÞ ei Þ obtained by Hammett and Perkins 43 previously used in the NFLF model, 6 (a can be related to v 1 by a ¼ 2j ðBÞ =3 ffiffi ffi 2 p v 1 ), which overestimates the thermal conductivity at moderate collisionalities around kk ðBÞ ei % 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 63 to closer fit the EIC results with one more fitting parameter (if the a 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. 64 Three Lorentzians [i.e., N ¼ 3 in Eq. (16)] can approximate our new fit to within 2.5% up to kk ðBÞ ei % 1:6; using six Lorentzians allows this to be extended up to kk ðBÞ ei % 30. The coefficients used are given in Table IV and were obtained using the variable projection method, 65 constrained such that Eq. (23) is obeyed to second-order in kk ðBÞ ei .

A. Homogeneous density and ionisation
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 0 and 0:15T 0 ; these were connected by a cubic ramp given by where n 0 c 2½À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=Q fl ¼ 1=Q ðBÞ þ 1=f L Q fs . 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.

B. Spatially varying density and ionisation
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 gasfill 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 46 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 T e ¼ 2:5 keV; n e ¼ 5 Â 10 20 cm À3 ; Z ¼ 2t o T e ¼ 0:3 keV; n e ¼ 5 Â 10 21 cm À3 ; Z ¼ 39 across approximately 100 lm 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 K ei ¼ log K 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 K ei % 8upto$ 1.6 mm in the hotter corona. Reflective boundary conditions were used here as in Sec. VA 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. VA, 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 lm leading to a maximum temperature change of 41% as shown in Fig. 7.
FIG. 6. Ratio of peak heat flow to that predicted classically for each snapshot against inverse scalelength k ei =L T (calculated at the location of maximum heat flow predicted by each model) for the nonlinear temperature ramp using different initial gradients.
Phys. Plasmas 24, 092309 (2017) 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 8 used a geometrically averaged mfp k e ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi Zk Ã ei k ei p . 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. 1

VI. DISCUSSION AND FURTHER WORK
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 righthand side of Eq. (15) by adding on nonlinear terms such as eE k m e @df @v k À P n A n @w n @s k 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 a 2 k 2 ei r 2 term in Eq. (17) with a 2 ðk ei rÞ Á ðk ei rÞ 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 electronelectron collisions on the anisotropic part of the distribution function f 1 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 n [which incorrectly predicts bðZÞ=bð1Þ ¼ j ðBÞ ðZÞ=j ðBÞ ð1Þ ¼ n 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. VA, 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 n 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. VB, 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. VB 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 K.Wehavesimplyshown 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. 66 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. 67 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. 71 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; 76 and while a few reduced models have been suggested to 092309- 11 Brodrick et al.
Phys. Plasmas 24, 092309 (2017) capture some of these aspects 9,11 they still need to be properly validated with kinetic codes.

VII. CONCLUSIONS
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 rÁg ðmbÞ 1 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.

ACKNOWLEDGMENTS
The authors would like to thank J. Y. Ji for sharing the numerical results of his closure, 63 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