We compare the reduced non-local electron transport model developed by Schurtz et al. [Phys. Plasmas 7, 4238 (2000)] to Vlasov-Fokker-Planck simulations. Two new test cases are considered: the propagation of a heat wave through a high density region into a lower density gas, and a one-dimensional hohlraum ablation problem. We find that the reduced model reproduces the peak heat flux well in the ablation region but significantly over-predicts the coronal preheat. The suitability of the reduced model for computing non-local transport effects other than thermal conductivity is considered by comparing the computed distribution function to the Vlasov-Fokker-Planck distribution function. It is shown that even when the reduced model reproduces the correct heat flux, the distribution function is significantly different to the Vlasov-Fokker-Planck prediction. Two simple modifications are considered which improve agreement between models in the coronal region.
In laser-produced plasmas relevant to Inertial Confinement Fusion (ICF), the electron temperature gradients are often so strong that the scale-length of the temperature variation becomes comparable to the mean-free-path of the electrons which transport thermal energy. As a result, the heat flux cannot be determined by the local conditions in the plasma and the finite electron mean-free-path must be taken into account (i.e., the thermal energy transport across the gradient is no longer purely diffusive). In these situations, the heat flux can be calculated with Vlasov-Fokker-Planck (VFP) models,1 which resolve the motion, scattering, and acceleration of electrons, but with much greater computational cost in comparison to local thermal transport models.
A widely used alternative to the full VFP models is the reduced multi-group model proposed by Schurtz, Nicolaï, and Busquet (the “SNB” model).2 Marrocchino et al. have also compared the SNB model to a VFP model which retains higher order terms in the expanded distribution function.3 The purpose of this paper is to compare the SNB model to two further test cases relevant to ICF: (1) a “burn-through” problem in which a heat-front breaks out of a high density region into a lower density gas and (2) a physically realistic simulation of a laser-ablated plasma over hydrodynamic time-scales. These scenarios are typically more challenging for Vlasov-Fokker-Planck simulations codes (than the usual hotspot relaxation problems studied previously) due to the presence of large density gradients.
We are also interested in checking the applicability of SNB-like models to study transport effects other than heat flow, for example, non-local Nernst advection of magnetic field and non-local corrections to electron Landau damping rates (which are important for determining Stimulated Raman Scattering gain spectra). These effects rely on a correct representation of the electron distribution function, so we also present some comparisons of the SNB and VFP distribution functions for the simple case of a linear temperature gradient. We find that even when the SNB model predicts the correct heat flow, it may not predict the correct underlying electron distribution function (somewhat paradoxically).
II. COMPUTATIONAL MODELS
The VFP model used in this paper (“K2”) solves the Vlasov-Fokker-Planck equation for the electron distribution function expanded in spherical harmonics in momentum-space to arbitrary order. The details of the K2 model are described in the Appendix; here, we give a brief overview. The electron VFP equation is coupled to the equations of radiation hydrodynamics which account for ion motion, PdV work, electron-ion equilibration, and radiative cooling. These effects are coupled to the electron VFP equation via heating and cooling operators, and the electron distribution function is advected with the ion background to maintain quasineutrality. The electron VFP part of the model accounts for thermal transport, including the determination of a self-consistent electric field and return current. Laser absorption is modelled by solving a simple ray equation, with the heating of electrons determined by an inverse-bremsstrahlung heating operator (which accounts for the distortion of the electron distribution function due to electron oscillation in the laser field).
We use the iterative/implicit algorithm proposed by Cao et al.11 to solve the SNB multi-group energy transport equations, with only two iteration cycles. We see no significant difference in results when using more iteration cycles in a few test runs, and this is probably due to the fact that the timestep used is significantly shorter than typical hydrodynamic simulation timesteps (as we resolve the electron-ion collision frequency). We have developed two SNB models independently and benchmarked them against each other, for a number of test cases.
Our implementation of the SNB model contains a subtle modification to the mean-free-path which is different to that in the original model. For a particle with energy ϵg (where “g” is the energy group index), we use the geometrically averaged mean-free-path where (in cgs units)
and . Our definition can be compared to Eq. (23) in the original SNB model.2 We find that the use of this mean-free-path improves agreement between the SNB and VFP models in high-Z plasmas. See Ref. 14 for a full discussion of this choice of mean-free-path (in the notation used in Ref. 14, our model is equivalent to the choice of “r = 2”). The reason for this modification stems from the fact that the original SNB model employs a Krook electron-electron collision operator. Since a Krook operator is only an approximation to the Fokker-Planck operator, there is some degree of freedom in choosing the exact form of the collision frequency. The factor , originally introduced in Ref. 15, is a common feature in VFP simulation models, allowing them to reproduce Spitzer's thermal conductivity in the low-Z limit without the need for an anisotropic collision operator. This allows the SNB model to more accurately account for the effect of electron-electron collisions in the non-local contribution to the heat flow.
All of the following test problems use f0 and f1 only (i.e., higher order harmonics are not included), except those in Sec. IV, which also include f2. We verified in each case that the inclusion of higher order terms did not produce any significant differences in the results. Although the code is two-dimensional, the test cases considered in this paper are one-dimensional (and only require a Legendre expansion in momentum-space). Note that the same Coulomb logarithm is used in both the SNB model and the K2 VFP model in each simulation (see below).
III. SIMPLE HEAT-BATH PROBLEM
We are interested in comparing the electron velocity distribution function (EDF) in the SNB model to that given by the K2 VFP model, as this gives us insight into how useful the SNB model could be in determining transport effects other than non-local heat flow in the unmagnetized regime. For example, heat flow in magnetized plasmas involves thermoelectric coefficients which are not the same velocity moments of the EDF as the thermal conductivity coefficients, and the damping of electrostatic waves is determined by the shape of the EDF near the wave phase velocity. These thermoelectric coefficients, along with other moments of the EDF, also contribute to the electric field (which in turn affects the magnetic field) through an effective Ohm's Law.
Our first test case is a simple non-linear heat-bath problem, in which the initial temperature (plotted in Fig. 1) is 1000 eV for , with a linear ramp over down to 100 eV (over 100 computational spatial cells). The typical ratio of the electron mean-free-path to temperature-gradient scale-length is . The total computational box size is . The electron density and charge state are fixed at , and the coulomb logarithm was held fixed throughout, . The relatively high value of was chosen to complement the comparisons in Ref. 3 (which used the values ). We plot the electron temperatures after 80 ps in Fig. 1 and the VFP heat flux is compared to the SNB heat flux in Fig. 2, showing reasonably good agreement for the characteristic reduction of the peak relative to Spitzer. In this test-case, we do not include harmonics above f1 (f2 is set to zero). Simulations including f2 do not show any significant differences in the heat flux. The temperature profiles in Fig. 1 indicate that the SNB model tends to overpredict the preheating. This type of simulation comparison is not new (hotspot relaxation was simulated and compared to SNB in Ref. 3), but we address the question of how well the SNB model predicts the energy distribution of heat-flux-carrying electrons in simple scenarios such as this.
The distribution of heat flux for this test case, , is plotted in Fig. 3 at , where the K2 and SNB heat fluxes are in very good agreement. However, the underlying distributions that give rise to this heat flux are significantly different. This indicates that some modifications to the SNB model may be required in order for it to reproduce other aspects of transport (mentioned above). In particular, we note that the heat flux according to K2 is carried by higher energy electrons on average, and K2 predicts a significant return heat flux (associated with the return current). Although return electrons do not contribute significantly to delocalization, their population is enhanced when the heat flux exceeds the local value due to the enhanced current carried by the energetic forward electrons.
In the original paper,2 no way of recovering f1 is given explicitly. In general, if f1 is known, the heat flux can be calculated with the following relationship
where m is the electron mass and v is the velocity magnitude. This integral is numerically represented as a sum
(where g numbers the velocity group vg and is the velocity bin size of group g). The SNB model expresses the heat flux, , in terms of the Spitzer value and a non-local correction
where Hg is the principle function to be solved for in the SNB model, defined in Eq. (27) of Ref. 2, and is the electron mean-free-path of energy group g, defined in Eq. (23) of Ref. 2 (and adjusted in our formulation, as mentioned in the previous chapter).
The Spitzer heat flux appearing in the above equation can be expressed in terms of the Spitzer expression for f1 (denoted ) as
In the above, Te is the electron temperature, is the Maxwell-Boltzmann distribution, and νei is the electron-ion collision frequency.15
Again, the integral in Eq. (5) is represented numerically as a sum
This expression can be substituted into Eq. (4), allowing us to rewrite it as
and we therefore infer a possible definition for which could be the expression
In other words, we have assumed that the quantity which is integrated to get the heat flux is . This definition for f1 is consistent with the Vlasov-Fokker-Planck derivation of the SNB equations, which give [see Eq. (56) of Ref. 2 and the text preceding it].
IV. BURN-THROUGH PROBLEM
In this test case, we consider a thin “target” with peak density ( for light) which is heated on the RHS by a laser ( wavelength), with a linear 1 ns rise time up to a peak intensity of . On the LHS of the target, the electron density drops to a constant low value , representing a gas fill. The density is held fixed throughout the simulation (hydrodynamic motion is turned off). The initial electron temperature was 50 eV, and the coulomb logarithm and charge state were held fixed throughout the simulation. Harmonics up to and including f2 were used in this test case. The initial density profile is shown in Fig. 4, along with the temperature (according to K2) at 70 ps, showing the hot coronal plasma on the RHS. Again we use 100 computational spatial cells.
The temperatures according to K2 and the SNB model are plotted at 170 ps (after the heat front has penetrated the shell and entered the gas) in Fig. 5.
The SNB model overpredicts the preheat of the shell (similarly to the “linear ramp” test-case), as is evident from the increased temperature at . The initial preheat of the gas, close to the gas-shell interface at , is also overpredicted by the SNB model. The dip in temperature at , apparent in both models, is due to hot electrons traversing this region without depositing much energy. Interestingly, the K2 model predicts enhanced long-range heating of the gas (in the region ), probably indicating a breakdown of diffusive transport. This is consistent with the fact that we find that f2 makes significant contributions to the heat flux in this region, which will allow for a strong departure from equilibrium. As a result of this, the heat flux in this region () is about twice as large as that predicted by the SNB model (at ). Contributions from f3 have been checked-for and found to be negligible.
Although this test case bears many similarities to the ICF scenario of capsule gas preheat, we caution against drawing too many conclusions regarding their similarity because the two scenarios differ in important ways: (a) the ICF shell has a much higher ρR, so the spectrum of hot electrons (generated by LPI in the corona) that reach the gas may be different, and (b) the gas is initially much colder than in our simulation, so a Fokker-Planck/Spitzer transport treatment is not valid (which is the reason we did not choose to tackle the full problem with our models). However, this problem represents an important benchmark for reduced non-local models, and highlights the interesting discrepancy in both short- and long-range heating.
V. HOHLRAUM ABLATION PROBLEM
In this test case, we consider the ablation of a Au-like wall next to a He gas, including hydrodynamics. We use an artificial material to represent the Au wall, which has a reduced solid mass density of (to enable us to run with a shorter timestep) and the ionization state is artificially capped at (which allows us to take into account the fact that the Thomas-Fermi ionization level is too high compared to non-LTE models). The solid Au material initially occupies the region , and He occupies the region . The laser pulse shape is shown in Fig. 6. The hydrodynamic part of the K2 model is run alone (i.e., with VFP and SNB transport turned off and Spitzer transport turned on) for the first 1500 ps to produce realistic density and temperature conditions (shown in Fig. 7), after which we allow the hydrodynamics to evolve with the heat flux calculated from the VFP solver. The SNB heat flux and temperatures are self-consistently evolved separately from the VFP temperature, although both “temperature” models use the same hydrodynamics heating/cooling terms (PdV work, radiation cooling, etc). The previous test cases were somewhat artificial/simplified in their design—here we attempt to give an indication of how the models compare in a more realistic scenario. This scenario is similar to the 2D test case compared against in the original SNB work (taken from Epperlein's work12), except we include the effects of hydrodynamic motion, electron-ion equilibration, ionization and radiation cooling, etc., and simulate the Au-gas interface and longer pulse lengths, which should be more directly relevant to hohlraum transport. Since simulations in this section are designed to be more physically realistic than the previous cases, the Coulomb logarithm is calculated locally at each time-step using the NRL prescription.13
The heat flux profiles are compared at t = 1640 ps (i.e., 140 ps after kinetic effects were turned on) in Fig. 8. The SNB model predicts the peak heat flux very well near the ablation surface (), but overpredicts the heat flux into the corona by a factor of . We find that this type of behavior persists over 100's of picoseconds. The overestimation of heat flux in coronal regions was predicted and commented on by Schurtz et al. in their original paper2 and also observed in Ref. 14.
The coronal (He gas) temperature in this simulation is relatively low (), so we also performed a simulation over a longer timescale and at higher laser intensity ( foot, peaking at ), giving typical peak gas temperatures . The initial conditions for the K2 simulation (i.e., just after the thermal transport is switched from the hydrodynamic model to the Fokker-Planck model) in this case are shown in Fig. 9. In this case, we find reasonable agreement between the K2 and SNB heat fluxes at the ablation surface (the SNB peak ablation heat flux is around 10%–20% higher than K2 throughout the simulation), although not as good as at lower laser intensity (see Fig. 10). Again, the peak heat flux into the corona predicted by K2 is a factor of lower than the SNB prediction.
In all hohlraum-relevant cases, we find that the peak flux is typically of the free-streaming value, as illustrated in Fig. 11 (also at t = 5300 ps). However, we caution against interpreting the value of as corresponding to the value of the flux-limiter that could be used with a flux-limited Spitzer model (necessary to reproduce the actual VFP heat flux). We also mention in passing that we have postprocessed data from fully integrated 2D radiation-hydrodynamics simulations of NIF-scale hohlraums and find there that the peak flux may be a higher fraction of the free-streaming value than the factor 0.035 in these simulations.
Note that in this section we computed the SNB heat flux from the VFP temperature profile, since the aim was to focus on the difference between the predicted heat fluxes for a given temperature and density profile. In Sec. VI, we explore how the heat fluxes actually affect the temperature profiles.
VI. INLINE COMPARISON
We showed in Sec. V that in hot hohlraums, the SNB peak heat flux near the ablation surface is overpredicted by 10%–20%. However, simply computing heat fluxes from given temperature profiles may not be a fair comparison of how well the SNB model performs when run inline (compared to an inline simulation with VFP). The argument can be made that if the SNB model overpredicts the heat flux in a region, the temperature will drop there at a faster rate, which will in turn lead to a reduction in the local temperature gradient and hence a reduction in the heat flux. Therefore, small errors in the model may be “self-correcting.” To explore this possibility, we have made a comparison between fully inline SNB and VFP (by “inline” we mean that the heat flux is used to update the temperature during the simulation). For this case, we freeze the hydrodynamic motion after a significant amount of ablation has occurred, resulting in the density profile shown in Fig. 12. After hydrodynamic motion is turned off, the laser continues to heat the plasma which then evolves to an approximately steady-state via a combination of thermal transport and radiative cooling over . The corona becomes hot enough in this steady-state () to drive significant non-local transport.
The temperature profiles corresponding to this steady-state are shown in Fig. 13 for each of the two models. As can be seen, the temperature profiles are in good agreement in the ablation/absorption region and only differ by in the coronal region. This suggests that typical errors on the order of 10%–20% in the peak ablative heat-flux predicted by SNB can effectively “self-correct.” We expect that the improvements discussed in Sec. VII will further reduce the coronal error.
We have not performed all SNB calculations inline in this paper due to the very large computational demands of resolving the electron collision time in the high density regions of the plasma over the long timescales needed for an inline comparison—the results in this section simply serve to bring attention to the possibility of “self-correction” in a relatively simple scenario. Improvements to the standard computational techniques used in VFP codes will be necessary in order to allow them to overcome the timestep limitation. This highlights the fact that one of the major design advantages of the SNB model is its natural ability to relax to Spitzer transport in regions of high collisionality (which is due to its expansion of f about the Spitzer distribution function). A similar expansion (or equivalent relaxation process) of the VFP equation should be possible. A related problem is the difficulty of maintaining quasineutrality in the high density regions in VFP simulation models without resolving the plasma period: implicit electric field solvers do not guarantee quasineutrality over long (hydrodynamic) timescales, but this problem can be largely overcome by the fully implicit simulation technique employed by the IMPACT model.17
VII. IMPROVEMENTS TO THE SNB MODEL
We have found that it is possible to increase the agreement between the SNB model and the VFP model in the coronal regions by simply modifying the mean-free-path definitions. The SNB model employs a somewhat artificial combination of electron-electron (λee) and electron-ion (λei) mean-free-paths into a single mean-free-path (“λe”). This is only possible when the ionization is assumed to be constant. The use of separated mean-free-paths (considered in detail in Ref. 14 and alluded to in Ref. 16) is a straightforward modification of the original model, and it correctly accounts for ionization gradients. The resulting heat flux is shown in Fig. 14 for the higher intensity case described in Sec. V. The heat flux is significantly improved in the coronal region, with the error in the peak coronal heat flux reducing from about 90% to about 27%.
A very similar improvement in the heat-flux is obtained by simply neglecting the “electric field correction” to the mean-free-path. The original SNB paper2 distinguishes between the standard collisional mean-free-path λg and an electric-field-corrected mean-free-path related via
which serves to further limit the particle stopping distance in regions where the electric field may be the dominant cause of electron stopping (see the discussion in Ref. 2). In the above expression, E is the electric field and ϵg is the energy of the electron (in energy group “g”). The green curve in Fig. 14 shows the heat flux when the electric field is neglected in the definition of the mean-free-path (i.e., by simply setting )—the level of improvement is similar to the use of separated mean-free-paths.
We have compared the SNB non-local transport model to fully kinetic Vlasov-Fokker-Planck simulations in a variety of realistic situations relevant to inertial confinement fusion. In hohlraum-like situations, we find good agreement between the SNB model and VFP simulations at the ablation front (where the heat flux peaks). The heat flow out into the coronal region, however, is typically overpredicted by the SNB model by a factor of . This overprediction can be substantially reduced by ignoring electric field corrections to the mean-free-path or properly separating the electron-electron and electron-ion mean-free-paths. However, the distribution function predicted by the SNB model may be significantly different to the VFP distribution function even when the total heat fluxes are in good agreement, suggesting the improvements need to b-e made in order to make the SNB model more generally useful as a predictive tool for a wider variety of kinetic effects. In the more demanding case of hot electron penetration through a thin shell into a gas, the SNB model significantly overpredicts the short-range heat flux while underestimating the long-range heat flux.
This work was performed under the auspices of the U.S. Department of Energy by LLNL under Contract DE-AC52-07NA27344.
C.P.R. and J.B. contribution to this work was funded by EPSRC Grant No. EP/K504178/1. C.P.R. contribution was also funded by EPSRC Grant No. EP/M011372/1. C.P.R. contribution 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-IFECCFE-01). The views and opinions expressed herein do not necessarily reflect those of the European Commission.
We would like to thank M. Rosen for stimulating discussions regarding flux limiters.
APPENDIX: COMPUTATIONAL MODEL
The VFP model used in this paper (“K2”) solves the Vlasov-Fokker-Planck equation for the electron distribution function expanded in spherical harmonics in momentum-space to arbitrary order. The expansion coefficients evolve as4
where v is the electron velocity magnitude relative to the ion velocity, is the electron-ion scattering frequency,10 are the components of the electric force, are the components of the magnetic field, and the functions and are defined as:
where is the magnitude of electron momentum. The electron-electron collision term (Cee) is composed of isotropic and anisotropic components:
where and , and (using Tzoufras' notation5):
Γee is defined in Ref. 5. The term CIB represents the inverse bremsstrahlung heating term, which is solved in the form given by Weng,6 and can be added into . The laser intensity is modelled in 1D by solving the 1D ray equation:
where is the inverse bremsstrahlung absorption coefficient, which is self-consistently calculated from the energy deposited by the term CIB. Similarly, in 2D the ray equations are solved by actually propagating rays, but without the self-consistent calculation of κ.
The term CMHD represents energy exchange between the MHD model and the VFP model, and includes hydrodynamic PdV work, electron-ion equilibration, ionization energy loss/gain, and radiation emission/absorption. These are computed by the MHD model and assumed to heat the electrons while maintaining a thermal (i.e., Maxwellian) distribution: deviations from this assumption for ionization have been discussed by Robinson,7 and for PdV work by Matte.8 We have found no deviation from a Maxwellian distribution when including the full electron-ion energy exchange terms (except in extreme conditions not relevant to this study). Our implicit assumption is that all hydrodynamic terms operate on time-scales significantly longer than the electron-electron thermalization time. In non-Cartesian geometry (cylindrical or spherical), the equations need to be modified to include the correct geometric factors as in Ref. 9.
We follow the KALOS formulism,4 time-splitting the VFP equation and integrating the advection and acceleration terms to second-order accuracy in time, space, and momentum, including polynomial expansions of the close to p = 0. The electron-ion scattering, isotropic electron-electron collision operator, and magnetic field terms are solved fully implicitly. The anisotropic collision operator is solved semi-implicitly. The irrotational electric field is found from an implicit solution of Gauss' Law and the solenoidal component from an implicit solution of the equations (i.e., a generalized Ohm's Law). We find that the use of Gauss' Law improves quasineutrality in the presence of large density gradients. Faraday's Law is solved for the magnetic field.
In this paper, we do not use the anisotropic electron-electron collision operator, but instead implement the simple electron-ion collision frequency multiplier given by Epperlein,15 since this is used in many VFP models as well as the SNB model. We use 180 energy groups in both the SNB and VFP models, and reflective boundary conditions in all cases.
K2 differs from SPARK10 and IMPACT17 by being fully explicit in space (which allows efficient parallelization), being capable of using an arbitrary number of polynomials in the momentum-space expansion, and by its inclusion of laser propagation and MHD. It also contains a number of other features, such as anisotropic collisions, EM wave propagation, 3D magnetic fields, and coupling to Vlasov ions, but they are not relevant to the work presented here.