Heat flux calculation in the semi-collisionless regime for substantial temperature variations including magnetic field

In the first part of the paper, we extend the hydrodynamic approach to improve the accuracy of heat transport calculations without a magnetic field. A new parameter and the corresponding evolution equation are introduced into the regular hydrodynamic model to describe the evolution of the electron distribution tail. The self-consistent E-field is found from the zero current condition. The resulting model is relatively easy to implement and shows a good agreement with kinetic simulations for up to kλei < 0.3. In the second part of the paper, we exploit the fact that the distribution function is close to the Maxwellian and solve the kinetic equation this time including the B-field. The resulting approximate solution allows us to find the expression for the heat flux across a B-field in the integral form. To simplify the flux calculation, we drop the ambipolar field term but include two normalization constants (whose values are defined against Fokker–Planck simulations) that account for the E-field effects and the B-field impact on the nonlocal E-field. This semi-analytical approach closely recovers the results of heat transport kinetic simulations in a wide range of collisionality and magnetization providing the way to account for the nonlocal transport across a magnetic field in existing hydrocodes.


Introduction
In the recent years the advances in the laser technology have made the concept of inertial confinement fusion (ICF) [1]- [3] more feasible. The most promising approach to ICF is fast ignition, in which the set of lasers is used to compress the ∼1 mm deuterium-tritium (DT) pellet into a ∼0.1 mm dot, and then an ultra-high intensity single laser is used to ignite the compressed fuel. The large amount of data from numerous experiments and different theories can be directly compared using numerical simulations. In spite of all the different new approaches to plasma simulation (e.g. based on quasi-particles like in PIC codes [4]- [7], or based on the solution of the Boltzman-Fokker-Planck equation [8]) the old hydrodynamic approach is still widely used. The reason for such resilient popularity is the speed and the simplicity of the code. In this paper, we extend the applicability region of the hydrodynamic approach by introducing new expressions for the heat flux which is more accurate in the semi-collisionless fusion plasma.
In the hydrodynamic approach instead of specifying the position and velocity of each particle with mass m, we use the set of 'velocity moments' (density n, average velocity, temperature T ) of the electron and ion distributions to describe the plasma. These quantities are defined as integrals of the distribution function (normalized to the particle density) multiplied by the appropriate power of the velocity, i.e. (1, v/n, mv 2 /2n) correspondingly, over the velocity space. Usually the description of the plasma assumes that the velocity distribution function is Maxwellian and, therefore, this description of plasma breaks down when the electron distribution function (EDF) substantially deviates from the Maxwellian distribution (from this point on we will be talking only about electrons since their contribution to the transport is (m i /m) 1/2 times larger than that of ions, assuming that ion and electron temperatures are equal; m i and m are the ion and electron masses). The Maxwellian distribution corresponds to the state of a local equilibrium which implies the smallness of the electron free path between collisions compared to the space scale L of any variation. Most of the electrons have the collisional free path close to the mean-free path λ ei = τ ei √ 2T/m, where τ ei = 3 √ mT 3/2 /(4 √ 2πZ e 4 n) is the electron mean collision time; is the Coulomb logarithm; e and T are the electron charge and temperature; n is the electron density; Z is the ion charge number. Therefore if kλ ei< 10 −1 , where k = 2π/L, most of the electrons experience on average more than ten collisions before they reach a region where the value of the distribution function parameters are substantially different. These collisions bring the incoming particles into local equilibrium with others. If however we consider the process of heat transport we find that the contribution to the heat flux comes mainly from 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT electrons with energy ε ∼ 5T (in the absence of magnetic field) whose free path is much larger. That can be understood as follows. The electron diffusion coefficient (for electrons with energy ε = mv 2 /2) is approximately D e (ε) ≈ ετ(ε)/m, with τ ∝ ε 3/2 , and the energy diffusion flux is approximately εD e ∂f/∂x, where f is the EDF whose integral over the whole velocity space gives us the electron density. We integrate the energy flux over all velocities to obtain the integral heat flux q ∝ exp(−ε/T )ε 5 dε , where we assumed that f is Maxwellian. The expression under the integral peaks at ε = 5T . Therefore, we can use the hydro model to describe the heat transport (classical Spitzer-Harm expression for the heat flux [9,10]) only if electrons with ε ∼ 5T can reach equilibrium on the space scale smaller than the inhomogeneity space scale L. The time of energy equilibration is about the time between electron-electron collisions but the distance travelled by an electron within that time is limited by its diffusion through ions and we find that energy equilibration length is λ T = λ ei √ Z. Therefore, if we want to use the hydrodynamic approach to describe the heat transport the energy equilibration length λ T of electrons with ε ∼ 5T has to be much smaller than the temperature variation length L, which amounts to kλ T< 10 −2 . For a very short wavelength inhomogeneity of electron temperature, kλ T > 10 −2 , the 'tail' electrons with ε ∼ T propagate far (λ T > L) without equilibrating and their distribution function may strongly depart from the Maxwellian.
Thus, for the heat flux calculation we only need to know the EDF around ε ∼ 5T (for the case of no magnetic field), where the main portion of the heat flux is carried [11], and the knowledge of the full distribution is rather excessive. Therefore, one can improve accuracy of the hydrodynamic description of the heat flux for kλ T> 10 −2 by using one additional parameter W, along with the electron density and temperature, for a more accurate description of the tail distribution around ε ∼ 5T . Some approaches [12]- [14] to the calculation of the nonlocal heat flux are based on finding the full energy distribution which significantly complicates the solution. On the other hand in the linear Grad's method [15]- [17] only two scalar moments (corresponding to the convolutions of EDF with {1, v 2 }) are used for the description of the isotropic part of the EDF, which, thus, is always Maxwellian. In section 2, we develop the approach to improve the accuracy of a hydrodynamic calculation of the heat flux for 0.01< kλ T< 1. For kλ T> 1 and large temperature variations the EDF is strongly non-Maxwellian and, generally, no conveniently small set of macro parameters can describe the particle transport. However, there are works [18,19] which describe the heat flow in plasma using a moment method for special cases, e.g. two-wall configuration, for arbitrary kλ ei . It is worth noting that the approach of an 'additional parameter' to the calculation of the nonlocal heat flux is not limited to small temperature perturbations T T which have been extensively studied in a number of papers (e.g. see [20,21]). The alternative approach to the heat flux which accounts for the kinetic effects is the nonlocal approach [22,23], where the limit on the temperature profile steepness improves to kλ ei < 1 due to the assumption that the isotropic part f 0 of the EDF is Maxwellian only for ε/T<1, while for the super-thermal electrons with ε/T 1 f 0 is found from the solution of the Boltzmann kinetic equation. This approach has been implemented in a number of works [12,13], [21]- [23] for the case of no magnetic field. However, in many cases including the inertial and magnetic fusion and astrophysics the magnetic field effects play an important role in the electron heat conduction (see for example [20,24,25]). Magnetic field effects on the nonlocal heat flux were considered in [20], but only for small temperature perturbations T T . Therefore, in section 3 we extend the results obtained in [12,13] to incorporate the impact of the magnetic field on the nonlocal electron heat transport.

4
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT In magnetized plasma, λ ei /ρ Be 1, the electron cyclotron radius ρ Be = (2T/m) 1/2 /ω B replaces the mean-free path as a characteristic transport scale (here ω B = |e|B/(mc) is the electron cyclotron frequency; B is the magnetic field strength; and c is the light speed), and the main contribution to the heat flux comes from electrons with energy around 2T for λ ei /ρ Be 1, because the electron diffusion coefficient gets much smaller for higher ε, i.e. D e (λ ei /(ω B τ ei ) 2 ) ≈ vλ ei /(ω B τ ei ) 2 ∝ ε −1/2 (instead of ε 5/2 ) and the integral heat flux q ∝ exp(−ε/T )ε 2 dε . (Note: if one were to account for the return current the maximum contribution would shift approximately to ε ∼ 6.5 and 2.5T for the cases λ ei /ρ Be 1 and λ ei /ρ Be 1 correspondingly.) The classical approach for the heat flux across homogeneous B-field remains valid as long as kρ Be 1 [24], even for large kλ ei . In section 3, we find an expression for the heat flux across a magnetic field with the limit of applicability extended to min(λ ei , ρ Be ) < 1/k (1.1) without limiting the scope of the problem to small perturbations of the plasma temperature. We summarize the results from the two different approaches and draw the conclusions in section 4.

Governing equations
We start with the Boltzmann kinetic equation for electrons without magnetic field (which is simply the continuity equation in six-dimensional space (3D + 3V) for the EDF, describing the density of electrons with given velocity) where E is the electric field; C ee and C ei are the electron-electron and electron-ion collision terms, correspondingly. We use the approximate form of C ei and C ee from [24] We approximate the EDF with a truncated spherical harmonics expansion (so-called 'diffusive The anisotropic part of the EDF b · v reaches equilibrium through electron-ion collisions, which are Z times more frequent than the electron-electron collisions responsible for the energy distribution f 0 equilibration. Therefore, for slow time variations ∂(b · v)/∂t (b · v)/ν ei we neglect the ∂(b · v)/∂t term in equation (1), but keep the ∂f 0 /∂t term. Taking the zeroth and first came from C ei term, and after the substitution of equation (2.3) into equation (2.4) we can get the equation for f 0 . We have not yet done anything new; basically up to this point we have followed in Braginskii's [24] footsteps. Now we write f 0 as a sum of f M (t, x, ξ) and the correction function is a some function of ξ , the subscript 'S' points at functions and variables whose values depend on the choice of function P S . The new parameter W allows us to correct the 'tail' of EDF without changing the particle density and the average energy. The equation for the Nth energy moment of f 0 S ≡ e −ξ ξ N P S (ξ) dξ; is obtained by integrating equation (2.4) multiplied by ε N over the velocity space where S ]/(τ ei Z) . The factor Z/(Z + 1) in the expression for C A comes from the electron-electron collisions happening on the timescale of τ ei /Z. These collisions 6 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT increase the effective 'collisional' ion charge by one. To calculate the integral in equation (2.9) with N > 1 we have used the high energy ( /T 1) approximation for C ee from [24] introduced earlier, since the main contribution into the value of the integral in this case (N > 1) is coming from high energies. The vector A (N) is the flux of the (N-1)th energy moment, e.g. A (1) is the electron flux, A (2) is the heat flux. Assuming that there is no electron current we get A (1) = 0, and find the expression for the electric field from equation (2.9) using N = 1 (2.10) To make the particle density and the average energy independent of W, we require that the correction function P S (ξ) satisfies 5nT . This way the transport coefficients and fluxes depend on the parameter W and its derivatives but the values of n and T are independent of W.
The function f 0 is determined by three variables n(t, x), T (t, x), W(t, x), and the three corresponding evolution equations (for zeroth, first and Lth moments of f 0 ) in one-dimensional (1D) case (and in the absence of the electron current) are obtained from equation (2.7) ∂A (2) ∂x , where L > 1; ν W ≡ ν L W is the decay rate of W due to collisions. (Note: the presence of an electron current would add a few terms in equation (2.12) which, if needed, could be found from equations (2.7) and (2.9).) The rate of Maxwellization of super thermal electrons with /T ∼ 5 is much slower than for electrons with /T ∼ 1, ν W 1/τ ee . And, though, we limit the rate of change of W, ∂W/∂t W/τ ee , to have f 0 not far from f M , we keep ∂W/∂t term because it is possible that ∂W/∂t ∼ ν W W . If we set ∂W/∂t = 0 we would obtain the quasi-equilibrium model, somewhat similar to the nonlocal models [12,13], describing the heat flux when the non-Maxwellian hot 'tail' electron 'fluid' is already in equilibrium across the temperature gradient while the thermal electrons with ε/T ≈ 1 are not.
The approximate analytical solution of the system in equation (2.12) for small perturbations presented in the following section reveals the main features and limitations of the model, before we proceed to the numerical modelling of equation (2.12) in section 2.3.

Linear analysis
In this section, we apply our model to the decay of a small temperature perturbation to find the linear decay rate. The (N-1)th moment flux A (N) given in equation (2.9) can be rewritten as (2.14) In the linear approximation the coefficients [C A T N+3/2 K (N) i ] (i = n, T or W) in equation (2.13) are assumed to be constant and that simplifies equation (2.12) to where X i and D i are the appropriate diffusion coefficients derived from equations (2.12) to (2.14) We also assumed that the amplitude of W is small so K (N) n (W ≈ 0) ≈ 0 (no contribution to the flux from the density gradient) and X T (W ≈ 0) becomes the classical heat diffusivity X classic T = 256Tτ ei /(9πm). We look for the solution of equation (2.15) in the form of a decaying periodic perturbation T = T 1 exp(ikx − γt) and W = W 1 exp(ikx − γt). We find that the linear decay rate γ for a temperature perturbation with given k is where In the classical limit kλ 1 collisions neutralize W and equation (2.17) recovers the classical heat diffusion decay rate γ SH = k 2 X T . In the opposite limit kλ 1 we get γ/k 2 X T = const ≡ C γ (see figure 1), 0 < C γ < 1, which corresponds to the maximum heat flux reduction in the absence of collisional neutralization of W . For equation (2.12) to be linearly stable one has to have 0 < C γ and that requires X T D W > X W D T (which amounts to . This expression limits the subset of possible values of L for a given correction function P S (ξ) due to a linear stability. Though we used small temperature perturbation approach to find the linear decay rate of perturbation, the model itself (based on equation (2.12)) is not limited to small perturbations. In the next section, we use a test problem featuring large temperature variations and compare our results to Fokker-Planck (FP) simulations.
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 1. Linear decay rate of a temperature perturbation γ/γ SH (ratio of the effective heat conductivity to the classical Spitzer-Harm [9,10] heat conductivity) as a function of kλ ei , where k is the wave number of the perturbation and λ ei is the thermal electron mean-free path. Curves A and K are the delocalization models from [12] and [13] respectively; filled circles relate to FP modelling from [26]. Curves P and D are the linear decay rates of the 'additional parameter' model with correction functions constructed on a polynomial basis or a δ-function basis, respectively.

Test problem: comparison with FP simulations
To start the simulation we have to choose the correction function P S (ξ) which determines the value of the coefficients J (N) S (N = 0, 1, L) introduced after equation (2.6). The choice of P S (ξ) is quite arbitrary, and we choose it to be the sum of three different functions (which have finite integrals with exp(−ξ)ξ N+1/2 , N > 0) with two free parameters p 1 and p 2 (which we need to make P S (ξ) orthogonal to 1 and ξ with weight function W = exp(−ξ)ξ 1/2 (see equation (2.11)) where for any given set of functions F i (ξ), i = 1, 2, 3.
If we use δ-functions to compose the correction function P S (ξ), F i (ξ) = δ(ξ − S i ), then on figure 1 but rather changes the level where it becomes flat (minimal value of γ/γ SH ). This set of values of S i and L (when P S is composed by -functions) seems to provide the least value of γ/γ SH in the limit kλ ei 0.1 while remain stable in the simulation of the test problem described below. If the correction function is a polynomial F i (ξ) = ξ S i , then k i = (S i + k + 3/2). The linear decay rate γ for this case is also plotted on figure 1 for P S (ξ) = ξ 1.5 + p 2 ξ + p 1 ξ 0.1 for L = 3. In both cases, for kλ ei> 0.5 the model reduces to the Spitzer-Harm model [9] with the flux limiter being equal to about 0.3.
We applied this heat flux calculation approach to the test problem, where we imposed a constant sinusoidal heating profile S(x) = (3/2)S 0 [1 + cos(2πx/L H )] on an initially flat temperature profile T 0 allowing the system to develop its 'natural' temperature profile as time goes on. We followed the time evolution of the amplitude of the relative temperature peak (T max − T min )/ T , where T denotes the average temperature. The following values of other parameters were used (typical for a laser plasma): T 0 = 1 keV, n e = 10 21 cm −3 , Z = 1 (in the high Z limit approximation), L H = 2 mm, S 0 = 4.5 keV ns −1 . The value of kλ ei calculated using the average temperature increases from about 0.1 to 1 following a T 2 law, due to the quadratic temperature dependence of the mean-free path (see figure 2).
For the heating to be Maxwellian (i.e., to obtain a distribution of the added energy according to a Maxwellian function) values of both T and W have to change, since For the test problem simulation we insert these heating terms into our system of equations given in equation (2.12).
Using the heat transport model given by equation (2.12) we calculated the evolution of the relative temperature peak (T max − T min ) T in the described test problem with the spatially modulated heating. The calculation results are plotted on figure 3 for two different forms of correction function (polynomial and δ-function). As one would expect from the analysis of figure 1 the polynomial and δ-functions descriptions are very similar and provide a rather good fit (less than 15% difference) to the FP simulation data (obtained using IMPACT code [8]) but only for kλ ei < 0.3. For kλ ei > 0.3 the model behaves as the classical hydro model with the flux limiter 0.3 (see figures 1 and 3). In the test problem kλ ei reaches 0.3 at t ≈ 0.15 ns (see figure 2) and according to the linear analysis our model starts to depart from the FP data [8].
The space distribution of T and W are plotted on figure 4 for t = 0.1 and 0.2 ns. Initially T is set to 1000 eV and W is set to be zero everywhere. As time goes on W diffuses from high to low temperature region establishing gradient of W in the opposite direction to gradient of T . Amplitude of W initially grows for ∼0.1 ns and then fades away to 0 (as T profile flattens). The gradient of W creates an additional heat flux (see equation (2.13)) which is going against the heat flux created by T gradient. The relative amplitude of the 'counter flux' is limited to the correction function P S (ξ)) and L (determines the equation for W(x)) seem to provide the maximum possible value of the 'counter flux'for given functions P S while keeping the simulation of the test problem stable.

Integral expression for the nonlocal heat flux across magnetic field
In the approach discussed in this section we assume that the super-thermal electrons (which are responsible for the most of the heat flux) reach quasi-equilibrium within the space scale exceeding the mean-free path much faster than thermal electrons do. Therefore if we are working on the timescale exceeding super-thermal electrons equilibration time (but less than the timescale of the thermal electrons equilibration) we can use quasi-steady state kinetic equation for the super-thermal electrons εT 1 where the expression for C ee and C ei are given in section (2.1). Assuming E = −∇ϕ we fuse two first terms on the left-hand side of equation (3.1) into one by introducing a new variable ε T = mv 2 /2 + eϕ (i.e. in the potential field the distribution function is the function of total energy To solve equation (3.2) we use the same approach as before and approximate f with a truncated spherical harmonics expansion f(r, v) = f 0 (r, |v|) + b(r, |v|) · v and take the 0th and 1st angular moments (multiply by 1 and v/|v|, and integrate over the full solid angle) of the equation (3.2) which provides us with two equations (where we neglected low order energy terms): From equation (3.10) one can easily recover equation (1.1) describing the applicability of our nonlocal approach. Finally, using the Green's function of the diffusion equation in the WKB approximation we find from equation (3.10) 3 ) −1 dp, and S 3 = 16(ω Be τ ei ) 2 /9π. The expression for the particle and heat fluxes comes from the integral over the anisotropic part of the EDF, which gives: In equation (3.12), the x-derivative ∂f 0 /∂x is calculated from equation (3.11), after which we can switch the derivative variable from x to −x , and then integrate by parts the product of the Green's function (which depends on both x and x through ξ and ξ ] by ∂f M /∂ε [which depends only on ξ (x ) ). After the substitution of equation (3.11) into equation (3.12) we get where eϕ = eϕ l + eϕ nl ; eϕ l = T(d ln n/dx − (5/2)d ln T/dx) is the electrostatic potential energy in a local equilibrium and eϕ nl is a correction due to nonlocal effects [12]; P i (x, x ) are the appropriate kernels derived from equation (3.12) (3.14) We can proceed further by setting the electron current j x to zero, which gives us the integral equation for the nonlocal electric field, which enters the expression for the nonlocal heat flux. However, the solution is expected to be very cumbersome, and for practical usage we can try to simplify equation (3.13). One of the possible ways is to approximate the heat flux q x (x) only with the first term under integral in equation (3.13) and account for the electric field effects through the normalization against the FP simulation. This assumption is justified if we would be able to closely reproduce FP results in a wide range of kλ ei (0.01 < 1kλ ei < 1) with only two normalization constants. We would need a second constant to account for the B-field impact on the nonlocal E-field. With this simplification we obtain the following expression for the nonlocal heat flux where The kernel P given in equation (3.16) has been simplified by calculating the integral over ε and it is a function of only two parameters P = P(G, S B ). Therefore, P is a constant surface in the space of the parameters ω B τ ei and G which can be interpreted as the effective magnetic field and distance correspondingly. Constants C A < 1 and C B are introduced to account for the neglected ambipolar electric field, where C A accounts for the impact of the nonlocal E-field and C B accounts for the change in the nonlocal E-field due to B-field. In general C A and C B are the functions of kλ ei but we expect these functions to be weak enough to allow the approximation with a constant in the wide range of kλ ei values.
We applied the expression for nonlocal heat flux given in equations (3.15) and (3.16) to the test problem described in section 2.3. In that problem T grows linearly with time (constant heating and no losses) from T (t = 0) = 1 KeV to T (t = 0.5 ns) = 3.2 KeV and, therefore (for a given constant B-field), the average magnetization ω B τ ei ∼ T 3/2 increases six times from its initial value at T (t = 0). We calculated a set of relative temperature peak (T max − T min ) T curves for different values of B-field. The values of the magnetic field were chosen in a way that at temperature T (t = 0) the values of ω B τ ei would be 0, 0.1 and 1. The time evolution of the 'nonlocality' parameter kλ ei /(1 + ω 2 B τ 2 ei ) with and without a magnetic field is shown on the figure 2. We tested our model against the FP simulation data obtained by using the IMPACT code [8] which also uses the Lorentz approximation. The best fit to the FP simulation for ω B τ ei ∈ [0; 6] and kλ ei ∈ [0.06; 0.6] was achieved with C A = 0.11 and C 2 B = 0.34 (see figure 5). These values of the constants are not too far from the values one would get by matching the classical [24] heat flux model in the limit kλ ei 0.01 at ω B τ ei ∼ = 0 and ω B τ ei 1(C * A = 0.167 and (C * B ) 2 = 0.154). The value of C A (the heat flux multiplier) corresponding to larger kλ ei is somewhat smaller (than C * A ) accounting for the presence of the additional nonlocal E-field (dropped out in equation (3.15) which is the simplification of equation (3.13)) inhibiting the heat flux. The value of (C * B ) 2 (the magnetic field multiplier) is less than C 2 B because the contribution from the additional E-field gets smaller as ω B τ ei increases. By slightly adjusting only two constants in our model we were able to closely reproduce three nontrivial curves (see figure  5, for (ω B τ ei ) t=0 = 0, 0.1, 1; and kλ ei ∈ [0.06; 0.6]) from the FP simulation [8]. We recommend using this model with C A = 0.11 and C 2 B = 0.34 in simulations including magnetic field where kλ ei> 0.01. We are not aware of any other works which provide the analytical treatment of nonlocal heat transport across B-field for strong temperature variations.

Conclusion
In this work we considered two different ways to improve the accuracy of the heat flux calculations based on the hydrodynamic model in the limit of weak collisionality. In the first Figure 5. Time evolution of (T max − T min)/ T according to classical [10] (Epperline) theory (curve CL), nonlocal model (curve NL with C A = 0.11 and C 2 B = 0.34) and the results from the IMPACT kinetic code [8] (curve FP) for different initial values of ω B τ ei (ω B τ ei = 0 corresponds to curve '1'; ω B τ ei = 0.1 to curve '2'; ω B τ ei = 1 to curve '3').
approach the addition of a one scalar moment of EDF to the set of classical hydro moments increases the maximum value of kλ ei in the applicability range by 30. The linear decay rate provided by the presented model is in a good agreement with kinetic simulation results up to kλ ei< 0.3 as oppose to kλ ei< 0.01 for the classical model (see figure 1). This approach given by equation (2.12) provides a rather straightforward (compared to delocalized models) and more rigorous (compared to flux limiters) way to account for the effect of the collisionless tail of EDF on the heat transport for kλ ei< 0.3. For kλ ei> 0.3 the model reduces to the regular hydro model with a flux limiter set to 0.3. For the cases when kλ ei can exceed 0.3 (or when the magnetic field is important) we consider another approach to the heat flux.
In the second approach we find the equilibrium distribution of hot electrons in the presence of B-field and obtain the expression for the heat flux containing unknown nonlocal E-field. Assuming that the relative contribution from the nonlocal E-field to the heat flux is approximately constant over a wide range of kλ ei values (but allowing it to vary with B-field imposed) we neglect the ambipolar field term, but include two normalization constants (whose values are defined against FP simulations) that account for the electric field effects and the B-field impact on the nonlocal E-field. The corresponding expression for the heat flux across B-field given in equations (3.15) and (3.16) is applicable as long as min(λ ei , ρ Be ) < 1/k and reproduce the FP simulation results quite well in a wide range of parameters. Overall this approach provides a way to account for the nonlocal heat transport across magnetic field in hydrocodes even when the perturbation approach [20] is not applicable.