A reduced kinetic method for investigating non-local ion heat transport in ideal multi-species plasmas

A reduced kinetic method (RKM) with a first-principles collision operator is introduced in a 1D2V planar geometry and implemented in a computationally inexpensive code to investigate non-local ion heat transport in multi-species plasmas. The RKM successfully reproduces local results for multi-species ion systems and the important features expected to arise due to non-local effects on the heat flux are captured. In addition to this, novel features associated with multi-species, as opposed to single species, cases are found. Effects of non-locality on the heat flux are investigated in mass and charge symmetric and asymmetric ion mixtures with temperature, pressure, and concentration gradients. In particular, the enthalpy flux associated with diffusion is found to be insensitive to sharp pressure and concentration gradients, increasing its significance in comparison to the conductive heat flux driven by temperature gradients in non-local scenarios. The RKM code can be used for investigating other kinetic and non-local effects in a broader plasma physics context. Due to its relatively low computational cost it can also serve as a practical non-local ion heat flux closure in hydrodynamic simulations or as a training tool for machine learning surrogates.


Introduction
The collisional transport of energy is a significant issue for the design and optimization of a wide range of schemes for realizing nuclear fusion in both magnetized and inertial systems [1].In the most commonly encountered context of ICF, where a fuelcontaining pellet is irradiated directly by lasers [2] or indirectly by a laser-produced thermal X-ray bath [3], the influence of heat flow is prominent in the transport of energy absorbed from the driver to the ablation front [4] and the hot spot.Here, the dominant carriers of the absorbed energy are delocalized (or 'free') electrons, which transfer their energy through collisions with ions and other free electrons in the irradiated material.If the temperature gradient is sufficiently shallow, this energy is transported deeper into the unheated ablator through conduction, i.e., the collisional 'local' diffusion of heat; the heat flux then obeys Fourier's law [5].However, with highintensity irradiation the mean free path (MFP) of the high-energy particles that contribute most to the heat flow routinely exceed the scale length of the gradient, resulting in transport with quasi-ballistic features [6].While 'non-local' behaviour may be accurately captured using full kinetic simulations [7], these methods are impractical for integrated modeling of the implosion process.Alternatively, and more pragmatically, numerous simplified models have been constructed which may readily be implemented as sub-cycled operators in radiation-hydrodynamics codes [8][9][10][11][12].At later times, when the fusion fuel is in a state of pre-ignition, heat flow also acts as a major factor in the power balance that determines the pathway and probability of success of fuel burn and the onset of ignition.
The foregoing arguments are usually framed from the perspective of the electrons, as it is commonly assumed that the ionic heat flux can be neglected in comparison to the electronic contribution.The assumption relies on estimating the ratio of the ion to electron heat conductivities by the square root of the electron-ion mass ratio.However, as recently discussed by Chapman et al. [13], for fully ionized fusion fuel one finds parity of the conductivities under modest temperature separations (around T i /T e ∼ 4 for deuterium), which can be vastly exceeded in burning fusion plasmas.Even more significant separations may be realized during the ignition and burn phase as the hallmark of this process is ion temperature runaway until hydrodynamic decompression occurs and the fuel explodes and rapidly cools.Non-local effects are also expected for ion heat flux since αparticle deposition will distort the hot tail of the background distribution function (DF), resulting in an over-abundance of energy-transporting ions relative to a Maxwell-Boltzmann distribution [14].
Hydrodynamic quantities such as density and temperature are populated predominantly by thermal particles in the bulk of the velocity distribution, whereas the heat flux is sourced primarily by suprathermal particles in the tail of the distribution.In a plasma, faster particles have longer MFPs, therefore non-local effects are most apparent in the tail of the DF where particle MFPs may be comparable to length scales on which hydrodynamic quantities vary.This leads to a regime where the bulk is dominated by collisions and near-Maxwellian, leaving hydrodynamic quantities undisturbed, but where the tail of the DF may experience non-local kinetic effects, perturbing the heat flux.A set of fluid equations describing evolution of hydrodynamic quantities in such a regime may then remain a valid description of the plasma, but with a non-local prescription for the heat flux.
In the local limit, ion transport coefficients for a simple plasma of electrons and a single ion species are provided by Braginskii [15].Comprehensive formalisms and practical prescriptions for calculations of local ionic transport processes in an unmagnetized plasma with multiple ion species are given by Chapman and Cowling [16], Ferziger and Kaper [17], Zhdanov [18], with compact formulae provided by Kagan and Baalrud [19], and Kagan and Tang [20,21].
The local result for the heat flux is generally a linear combination of gradients of hydrodynamic quantities such as temperature, but extending to other quantities such as concentration and partial pressure in multi-species cases.For a given density and temperature, there is a maximum possible heat flux corresponding to all particles streaming freely at the thermal speed.For a sufficiently sharp gradient, the local result will exceed the free-streaming limit, thereby become unphysical.Non-local effects then have a tendency to inhibit heat flux to remain physical, which is supported by simulations and Fokker-Planck (FP) codes [22].This observation motivates a flux-limiter approach [23,24], where local transport coefficients are limited in non-local cases by a fitted function compared to full FP methods.Although computationally cheap, flux limiters are physically inaccurate and fail to capture other features of non-local transport such as preheating in cold regions.More advanced models, such as the Schurtz-Nicolaï-Busquet (SNB) approach [25], provide a more physically accurate description of non-local transport for electrons but not ions.For ion transport, this leaves flux limiters as the usual practical option implemented in hydrodynamics codes.Chapman et al. [13] indicates that flux limiters may be insufficient for modelling nonlocal ion heat transport accurately and therefore may be a major contributor to discrepancies in modelling ICF systems of interest, in particular reproducing neutron yields from experiments on uniaxially driven targets [26,27].
In this work, we present a reduced kinetic method (RKM) based on approaches by McDevitt et al. [28][29][30] to describe non-local behaviour of the ion DFs, and therefore ionic heat flux, in a multi-species plasma.The RKM solves only for the tail of the DF where non-local effects are expected to be most pronounced.
Further noting that fewer particles reside on the tail compared to the bulk, a linearized Coulomb collision operator is valid even for orderunity deviations from a Maxwellian distribution in the tail.
Transport quantities such as the heat flux can then be computed from the DFs provided by the RKM solver.
This leaves a novel firstprinciples model for non-local ion transport, which as of now has not been widely studied even for a simple, single ion species plasma.Importantly, this study is conducted in the multi-species context where ionic heat flux includes additional features not present in the single ion species or electron cases.In particular, an additional contribution to the heat flux exists due to relative diffusion of different ion species which is generally of comparable magnitude to the standard conductive component.Multi-species scenarios also allow gradients of partial pressures and concentrations act as drives for transport fluxes along with the standard drive due to the temperature gradient.The impact of non-locality on both of these multi-species effects is for the first time investigated.We consider two components of the heat flux, the reduced heat flux and the enthalpy flux, and find that non-local deviations of the latter can be insensitive to sharp pressure and concentration gradients, therefore increasingly important in comparison to the reduced heat flux in non-local scenarios.Many hydrodynamics codes treat ion mixtures as a single effective species, thereby ignore the enthalpy component of the heat flux which may be particularly significant with sharp concentration gradients and therefore influence fluid simulations.
The remainder of this paper is organized as follows; in 2, the mathematical model as derived from the VFP equation and numerical implementation is presented.In 3, the additional complexity of multi-species heat flux in comparison to the singlespecies case is described.In 4, results from the RKM are presented for a single species and weakly and strongly asymmetric ion mixtures with gradients in temperature, pressure, and concentration, as well as consistency checks, convergence studies, and comparisons to SNB and flux-limited models.

Reduced kinetic method
Ideal, weakly coupled plasmas may be described by the VFP equation, which is the Boltzmann equation for a plasma with a Fokker-Planck collision operator as a consequence of long-range Coulomb collisions with multiple species.The VFP equation for a species α in laboratory-frame phase-space coordinates (t, x, v) for particle position x and particle velocity v is where f α is the DF of species α and F α is the force per unit mass acting on species α.The FP collision operator C αβ {f α , f β } represents collisions of species β acting on species α, with the sum over all species β including α. Ion-electron collisions are almost always negligible compared to ion-ion collisions and are therefore neglected, although they may be readily included.
It is convenient to decompose each species' DF into a Maxwellian component and its deviation from Maxwellian as , where n α is the species' number density, T α is the species' temperature in units of energy, and v th,α ≡ 2T α /m α is the species' thermal velocity with particle mass m α .The centre-ofmomentum velocity, or hydrodynamic velocity, is defined as u = α c α u α , where u α = 1 nα d 3 v vf α is the the net flow velocity of a species α.The difference between these different velocities and the hydrodynamic velocity is typically much smaller than their thermal speeds |u α −u| ≪ v th,α for an ion mixture near local thermal equilibrium, therefore the Maxwellian drift velocity for all species may be chosen to be u instead of u α for convenience.Although different species may have their own temperatures in the mathematical framework, ion species with comparable masses and different temperatures is indicative of a system far from local thermal equilibrium, therefore a common ion temperature T α = T is taken.Using bilinearity of the FP operator, and noting that Maxwellians satisfy C αβ {f M α , f M β } = 0 for a common temperature and net flow velocity, the collision operator can be decomposed as The first term on the RHS of (2) which is numerically tractable while still capturing relevant kinetic physics.
It is often the case that the length scale of a boundary in a fluid system is much shorter than its spatial extent in the other two directions parallel to the boundary, therefore the geometry near the boundary may be considered to be near one-dimensional planar (1D).This requires that any scalar quantities of the problem (such as temperature and number densities) only vary in the direction perpendicular to the boundary.
In a one-dimensional spatial planar geometry, there is an azimuthal symmetry around the single spatial direction, therefore the velocity space becomes two dimensional (2V), as depicted in figure 1.
Transport quantities such as heat flux are defined in the centre-of-momentum frame.Appendix A gives the kinetic equation for the deviation of the DF from Maxwellian δf α = f α −f M α in the centre-of-momentum frame in a 1D2V planar geometry as Visualization of 1D2V geometry at a boundary.There is an azimuthal symmetry in the velocity coordinate ϕ, and translational symmetry parallel to the boundary.
where z is the single planar spatial coordinate, w = v − u is the velocity in the centre-of-momentum frame (peculiar velocity), w ≡ |w| is particle speed in the centre-of-momentum frame, and µ ≡ w • ẑ/w = cos ϑ is the cosine of the velocity pitch angle ϑ in the centreof-momentum frame.The diffusive driving force d α is defined in Ferziger and Kaper [17] as is the ratio of the species' mass density ρ α = n α m α to the total mass density of the mixture ρ = β ρ β .The partial pressure of species α is given as p α = n α T , and p = α p α is the total pressure of the ion mixture.The hydrodynamic velocity has been neglected by assuming it is much smaller than all thermal speeds |u| ≪ v th,α .(4) will then be relaxed to ∂ t δf α = 0 since we solve for transport quantities being functions only of the instantaneous hydrodynamic profiles.

Test-particle operator
The test-particle collision operator against a Maxwellian background [32] of species β may be written in the form is the pitch-angle scattering operator with azimuthal symmetry, and the three collision frequencies are defined via where ξ α = w/v th,α , erf(x) is the error function [33], where we take ln Λ αβ = ln Λ as the conventional Coulomb logarithm for simplicity.The latter may take many different forms as determined from theoretical arguments [34] or fits to high-fidelity plasma microphysics simulations, however throughout the remainder of this work the Coulomb logarithm will be treated as a tunable parameter of the RKM framework.With the foregoing considerations the test-particle collision operator ( 5) is then linear in species β, such that C T α = β C T αβ can be written in the more compact form where ν α = β ν αβ for each of the three collision frequencies ν α D , να s , and ν α ∥ .This collision operator appears similar to the single species form (i.e., if the collision operator (5) had only a single species present), but where the collision frequencies may have a more complicated v dependence in asymmetric mass and charge cases due to the error and Chandrasekhar functions appearing in (6).

Field-particle operator
Rosenbluth, MacDonald, and Judd [35,36] give the field-particle operator in 1D2V coordinates as where the Rosenbluth potentials H β and G β and their appropriate derivatives are given in terms of δf β as Legendre series over µ as described in 2.
The test-particle operator contains only velocity derivatives on δf α , therefore is local in velocity space and does not couple the deviations of different species.
In contrast, the field-particle operator contains integrals over velocity space and acts on δf β , therefore is non-local in velocity space and couples the deviations of different species.This makes the fieldparticle operator more difficult to include analytically and numerically than the test-particle operator, yet retaining the field-particle operator is required for agreement with analytic results in the local limit such as Braginskii's heat flux results, as well as for other important features such as momentum conservation.
Particles in the tail mainly experience collisions with the bulk of the distribution since it is the most populated region in velocity space.Since the tail is located far away in velocity space from the bulk, tail particles do not resolve the non-Maxwellian part of the bulk well and therefore effectively collide with a Maxwellian background.Such collisions are represented by the test-particle operator, therefore it is expected that the test-particle operator dominates over the field-particle operator in the tail of the distribution, which is the region of phase space which we are primarily interested in for non-local behaviour.This feature allows for a tractable numerical implementation of the field-particle operator as discussed in 2.4.

Expansion in Legendre polynomials
Similarly to previous work [28,37], δf α may be expanded in Legendre polynomials as where the Legendre components which follows from orthogonality of Legendre polynomials.Legendre polynomials are a natural basis to expand in over µ since they are the azimuthally symmetric versions of spherical harmonics which provide a complete orthogonal basis over the unit sphere.Integrating (4) over 2 2m+1 1 −1 dµ P m (µ) and making use of the orthogonality, recursion relations, and eigenvalue property Lµ P m = − 1 2 m(m + 1)P m of Legendre polynomials yields the set of PDEs where β } is given in 2, and δ l,k is the Kronecker delta.
Truncating the expansion to the first N l Legendre modes, the set of PDEs may be solved numerically.Fine detail of δf α in the pitch angle µ is not expected, so the series should converge for a reasonable number of Legendre modes.This is checked with a convergence study over N l , as well as having been demonstrated in previous work by McDevitt et al. [29] for a similar problem by comparing the solutions of the system of PDEs for the Legendre components to the full solution of the 3D equation without expanding in Legendre polynomials over µ.

Numerical implementation
The set of PDEs ( 12) is solved by discretizing δf (l) α (z, w) for all species together to a 4D array F s,l pq , where z = p∆z, w = q∆w, l is the Legendre index, and s is the species index.Derivatives are estimated via the finite differences where the cell-face values are interpolated via a linear scheme The coefficients δ z p,q and δ w p,q (which are functions of grid position) specify the linear interpolation scheme.δ z p,q = 1/2 is chosen for a spatial central difference and a 'Chang-Cooper'-inspired [38] scheme δ w p,q = 1/u − 1/(e u − 1), where u = 2(∆w) 2 q/v 2 th,p is chosen for velocity derivatives.This choice improves stability and convergence with respect to a simpler finite difference scheme with fixed δ w p,q .Reshaping the mesh F s,l pq into a vector F , the set of PDEs is discretized to ∂ t F = M (F ) + P , where M is the linear operator representing the discretized right side of (12) acting on δf (l) α , and P is the last inhomogeneous term reshaped as a vector.The stationary solution may be found via an explicit or implicit relaxation method.Here, SciPy's backward-differentiation formulas implicit integrator [39,40] is used to relax to the steady-state solution.
As discussed in 2.2, the field-particle term is negligible in the tail compared to the test-particle operator.This is computationally convenient since the Jacobian of the field-particle operator is nonsparse in the velocity index, whereas the Jacobian of the test-particle operator is.For the solver, the RHS of ( 12) without the non-sparse integral terms appearing in the field-particle operator may then be used as an approximate Jacobian, with the full fieldparticle operator entering along with the full timestep function.Although the field-particle operator requires computing Ψ (l) β,(1,2,3,4) (w) integrals as defined in 2 for each l, species β, and spatial position, the integrals are cumulative in velocity, therefore may be computed with reasonable efficiency (see 2).Alternatively, the Jacobian of the field-particle term can be included, and these two approaches have been found to be in good agreement.
The leading-order Chapman-Enskog (CE) solution [16,17] is chosen for the the boundary conditions for the DF at spatial edges (where spatial gradients become shallow, therefore transport should become local) and at the low tail cutoff (w min = ξ min v th,α , with ξ min ∼ 1).The spatial boundary conditions are motivated physically by arguing that over a sufficient length, hot particles responsible for non-local effects will have had many collisions and will not penetrate much further away from the boundary, so the local limit should be restored.The local solution is exponentially small for the choice of functions of the bulk profiles, so choosing f = f M α as a spatial boundary condition is nearly identical.However, the local perturbative solution does provide a useful boundary condition for the tail cutoff, since δf α would otherwise be zero for w ≤ w min , which would remove possibly important parts of the velocity integrand when calculating fluxes.Since the bulk particles have shorter MFPs, we still expect the bulk to match the local solution.At some maximum velocity cutoff above which contributions to the heat flux integrands are expected to be negligible, δf α = 0 is chosen due to the fast exponential decay.

Heat flux in a multi-species plasma
Considering the heat flux in a multi-species plasma has further complications in comparison to the single species case.A single species may only experience selfcollisions whereas multi-species transport requires also considering collisions with other species.Furthermore, different species in mixtures diffuse in the centre-ofmomentum frame which gives rise to an additional source of heat transport of generally comparable magnitude.These additional complications in the multi-species case are also apparent mathematically since there is a set of coupled kinetic equations to solve instead of one in the single species case.Experimental evidence for various types of ICF implosions indicates that multi-species effects can directly affect the yield [41][42][43][44], further motivating the present study.
For a plasma with N ion species, there are N coupled kinetic equations which evolve each species' DF.From this set of kinetic equations, a set of 5N fluid equations may be derived for each species' mass density, flow velocity, and temperature.Since the timescales on which different ion species exchange momentum and energy are much shorter than other (often hydrodynamic) timescales of interest, the ions typically have a single common temperature and similar mean flow.
There may be cases with highly asymmetric mass mixtures where such an approximation is not valid, but these are not considered in the present work.Such extensions could be considered in the future if the need arises.
The separate temperature evolution equations may be combined to form a single evolution equation for the common ion temperature, where the heat flux enters only by the species-summed total heat flux The heat flux for a species α is defined in the centre-of-momentum frame as Q α = d 3 w 1 2 m α w 2 wf α .In a simple plasma featuring electrons and a single ion species, there is no heat flux associated with the flow since the only ion species velocity is identical to the hydrodynamic velocity.On the other hand, in a multispecies plasma, each species α has an individual diffusive flow defined as V α = u α − u relative to the hydrodynamic velocity.This introduces a component to the heat flux in multi-species cases due to the diffusion of different species in the centre-of-momentum frame transporting energy, which does not occur in the single-species case.The heat flux Q may be decomposed into two parts via Q = q + h, the reduced heat flux and enthalpy flux respectively, defined as Although it is not necessary to decompose Q = q + h, it provides more points of comparison for the RKM with the analytic result in the local limit.
The CE expansion [16,17] of multi-species Boltzmann equations yields local results for transport quantities, in particular the heat fluxes in the centreof-momentum frame in which where λ ′ , D α , and D αβ are the partial coefficient of thermal conductivity, multi-species thermo-diffusion coefficients, and multi-species (ordinary) diffusion coefficients respectively, with the sign convention following Ferziger and Kaper [17].These transport coefficients are determined from the collision operator C αβ .For the full Coulomb collision operator with Debye shielding, these may be provided analytically.
Here the weakly coupled limit results of [19] are used, which provide numerical routines for the multi-species coefficients.
In the local single ion species case, the only two gradients available to induce heat flux are the temperature and pressure gradients.By assuming force balance, the single-species ion heat flux may be expressed purely in terms of the temperature gradient.In contrast, the multi-species case involves more gradients on which the heat flux may be dependent, such as concentration gradients.In addition, for single species cases the force F does not appear explicitly in the kinetic equation for δf due to force balance.For multi-species cases, the force F α may vary for different species, such as different specific charges eZ α /m α , giving different electrostatic forces on each species, therefore these forces will appear explicitly in the kinetic equation for δf α and therefore affect the heat flux.
In 1D2V coordinates, the diffusive velocities in the α , and the heat fluxes are and which vanishes for a single species since u α = u.

Single species case
Consider a single species with a number density profile n(z) and a temperature step of the form which transitions from T hot to T cold with increasing z over some spatial scale L. n 0 = n total (z = −∞), T 0 = T (z = −∞), m 0 = min α (m α ), and v 0 = 2T 0 /m 0 are introduced as the characteristic number density, temperature, mass, and velocity to which quantities are appropriately non-dimensionalized.In particular, heat fluxes are normalized to a characteristic heat flux Q c = n 0 T 0 v 0 , DFs are normalized to n 0 v −3 0 , and operators appearing in the kinetic equation are normalized to a characteristic collision frequency ν 0 = n 0 e 4 ln Λ/4πϵ 2 0 m 2 0 v 3 0 .The Knudsen number is introduced as N K = λ 0 /L as a measure of non-locality, where λ 0 = v 0 /ν 0 = 16πε 2 0 T 2 0 /n 0 e 4 ln Λ is the MFP of a thermal particle at z = −∞.The RKM solver is run for a temperature step T hot /T cold = 10 and a variety of Knudsen numbers from the local limit to non-local cases.The electric field has been set to zero to focus on the non-local effects due to the temperature profile, although it can be readily included.Figure 3 shows the reduced heat flux for an isobaric (p = nT = constant) and isochoric (n = constant) temperature step.In the local limit the two cases produce the same result and agree well with the local Braginskii result.The latter is also the same in the isochoric and isobaric case.In non-local cases, the isochoric heat flux forms a pedestal in the cold end, whereas the isobaric heat flux does not.This is expected since the higher density in the cold end of the isobaric case gives a higher collisionality, therefore fast particles streaming from the hot end do not penetrate as far.Figure 4 shows the first two Legendre components f (0) = f M + δf (0) and f (1) = δf (1) of the DF at the hot end and cold end for the single-species isochoric case.f (0) is the distribution function averaged over the velocity pitch angle θ.For small Knudsen number N K = 1/1000, the RKM result agrees with the local CE solution.For increasing Knudsen number, the tail of the DF in the hot end is further depleted and the tail of the DF in the cold end is further enhanced.This is expected physically at higher Knusden number as fast tail particles in the hot end may stream over the temperature gradient into the cold end, thereby reducing the tail population of the DF at the hot end and enhancing it in the cold end.
δf (z, w, ϑ) may be reconstructed via the Legendre series to visualize anisotropy in the DF as shown in figure 5.In the hot end, the tail of the DF is suppressed in the ẑ • w > 0 portion (noting ẑ ∥ −∇T ) as is expected from faster tail particles with longer MFPs streaming out of the hot end, thereby depleting its tail population.In the cold end, the tail of the DF is enhanced in the ẑ • w > 0 portion due to suprathermal particles streaming into the cold end from the hot end.It may be noted that this tail enhancement appears particularly anisotropic with the portion of the tail with ẑ • w < 0 appearing mostly unchanged in the hot end, since the suprathermal particles from the hot end must have been travelling in the ẑ •w > 0 direction originally to enter the cold end, and they are subject to little collisional scattering in the cold end since ν D is small in the tail.

Multi-species cases
For multi-species cases, we define the Knudsen number in a similar manner as the single-species case, but instead using a maximum MFP N K = λ max /L defined as λ max = max α,z (λ 0,α (z)), with the MFP of a species α at spatial position z accounting for collisions with all species via λ 0,α = v th,α / β ναβ .A 50:50 DT mixture and a fully ionized CH 2 mixture are chosen as multi-species examples to investigate both weakly and moderately asymmetric mixtures.Both of these types of mixtures are relevant to ICFrelated systems which often use DT fuel and plastic layers.A CH 2 mixture may be considered to be moderately asymmetric in mass, since the masses differ by an order of magnitude.Although the charge ratio is only Z C /Z H = 6, the charges enter the equations through the collision frequencies ν αβ ∝ Z 2 α Z 2 β , therefore the charge asymmetry is effectively on an order of magnitude.Three cases are chosen for each mixture: a temperature gradient at uniform pressure and concentration, a pure pressure gradient at uniform temperature and concentration, and a pure concentration gradient at uniform temperature and pressure.These cases with isolated ∇T , ∇p, and ∇c are selected to investigate the sensitivity of the locality of heat fluxes to different driving gradients in temperature, pressure, and concentration.Similarly to the single-species case, the electric field has been set to zero to focus on the non-local effects from the aforementioned gradients.These six cases are shown in figures 6 and 7 for DT and CH 2 respectively.
In a two-species ion mixture ignoring electric fields and electrons, the diffusive mass flux of the lighter species in the local limit may be written in the Landau form [21,45] where D is the classical diffusion coefficient, and k p and k T are the baro-diffusion and thermo-diffusion coefficients respectively.It has been indicated that for a two-species mixture, the enthalpy flux h is directly proportional to the diffusive flux V l , therefore plots which show the enthalpy flux are directly related to the diffusive flux.
All of the multi-species cases show the same basic qualitiative behaviour as the single-species cases; increasing Knudsen number inhibits the peak heat flux and enhances the heat flux in some regions away from the peak.However, the multi-species cases do demonstrate additional behaviour of interest.Firstly, non-local effects on the heat flux are demonstrated in the absence of a temperature gradient, induced instead due to sharp gradients in pressure or concentration.In many of the cases shown, non-local effects give rise to large heating pedestals in regions |z/L| ≳ 2 away from the boundary due to sharp gradients in T , p, or c.This is in dramatic contrast to the local heat flux in these regions, which is exponentially small for the chosen ∼ tanh(z/L) profiles.In most cases this enhancement only appears in one direction to the boundary, for example the ∇p cases 6b 7b where the non-local behaviour is only pronounced for z/L > 0. In contrast, the DT ∇c 6c case has significant non-local behaviour in both directions due to the near-symmetry, which is not present in any of the other cases.
In the temperature gradient cases 6a 7a, both components of the heat flux are similarly sensitive to non-local behaviour, and appear similar to the singlespecies case.Since the enthalpy flux is proportional to the diffusive flux, non-local behaviour in the former imply that diffusion can be sensitive to kinetic effects at similar Knudsen numbers at which heat transport is sensitive.In contrast, the DT ∇p case 6b and both ∇c cases 6c 7c have the reduced heat flux q being significantly affected by non-local behaviour, whereas the enthalpy flux h is quite insensitive to non-local behaviour.This leads to cases around the peak flux where q is significantly suppressed and h is not, so the latter may be increasingly important and dominant at the peak in non-local cases.

Electron transport
While the main focus of this paper is the ion heat transport, the newly developed approach can be applied to the electron case as well by treating electrons as an ion species with Z e = −1 and m e ≪ m i .Since all other non-local heat flux models available to date are for the electron heat flux only, we have run this calculation to be able to benchmark our method against an earlier, electron heat flux approach.

Local comparison to Braginskii
In the local limit, Braginskii's results for the electron heat flux may be reproduced for a quasi-neutral isobaric simple   where R e is the electron thermal force, with the local result including the electron thermal force coefficient B e appearing in (2.9) of [15].The terms including flows, friction and viscosity, have been neglected.In the local limit, the RKM reproduces Braginskii's results successfully as shown in figure 8.

Non-local comparison to SNB
The RKM heat flux may be compared to the heat flux as predicted by other non-local transport models, namely the SNB approach and flux limiters.Figure 9 shows a comparison of the non-local electron heat flux from the RKM solver and SNB approach for an isochoric simple hydrogen plasma with a Gaussian temperature profile.
For the chosen values n e = 10 26 m −3 , T max = 200 eV, T min = 100 eV, and ln Λ = 10, the characteristic electron mean free path at the hottest temperature is λ 0 = 6.16 µm.The SNB implementation used here follows the 'iSNB' algorithm from Cao et al. [46].The RKM and SNB hyperparameters are suitably chosen for convergence.Braginskii's result for the electron conductive heat flux q T e in the local limit chooses an electric field such that electron momentum conservation is satisfied, RKM results for the electron heat flux for various ionization numbers (Z i = 1, 2, 3, 4, 100) matched with Braginskii's results.Each dashed black line corresponds to the Z of the nearest colored RKM line.The RKM solver is run with two species, electron and ion, with m i /me = 1000 and Ze = −1.The temperature profile with Knudsen number N K = 1/1000 has a temperature drop by a factor of 10, and the number densities are chosen such that the mixture is quasi-neutral and isobaric.Similarly to Braginskii, an electric field ensuring electron force balance using the local result for the thermal force is included.
which includes the electron friction and thermal forces determined self-consistently from δf e .In non-local cases, the local expressions for electron friction and thermal forces are no longer valid, therefore the correct electric field yielding zero current is not known before running the solver.Root-finding methods are applied  to find the electric field in non-local cases such that the electron current, and therefore enthalpy flux, is zero everywhere.
In addition to SNB, a flux-limited heat flux is included for comparison.The flux limiter takes the simple harmonic-averaged form where q f.l. is the flux-limited heat flux, q Brag. is the local Braginskii heat flux result, q f.s.= n e T min v th,min is a characteristic free-streaming heat flux for the system, and α is the flux-limiter parameter.Here, a range of values 0.05 ≤ α ≤ 0.15 are chosen.The flux-limited heat flux fails to capture enhanced heating in the colder region in comparison to the SNB and RKM approaches.The first case uses the full linearized Coulomb collison operator with test-particle and field-particle operator.The second case uses a modified Coulomb collision operator where the field-particle operator is This modification is proposed for improving computational performance while retaining some physical accuracy as the expensive integrals in the field-particle operator are not computed each timestep, yet agreement in the local limit is retained.As demonstrated, this approximation yields excellent agreement with SNB.The SNB approach is based on a simplified collision operator, with this particular implementation based on a Bhatnagar-Gross-Krook (BGK) type operator C{δf (1) } = −ν ei (w)δf (1) with ν ei (w) ∝ w −3 .The improved agreement between the SNB and RKM approaches when the latter employs the fixed, rather than first-principle, field-particle operator suggests that the field-particle operator, which is nonlocal in velocity space, may play an important role in non-local behavior that is not captured in simpler collision models such as BGK that are local in velocity space.

Convergence studies
There are a few hyperparameters of the model which should be chosen for fast computational performance Convergence study for the number of Legendre modes included in the expansion δf ≈ N l −1 l=0 P l (µ)δf (l) .The solver is run identically to the isochoric single species example 3 with ξ min = 1, with the local limit results (N K = 1/1000) shown in dashed lines and two significantly non-local cases (N K = 1/20 and N K = 1/10) shown in dotted-dashed lines and dotted lines respectively.In the local cases, the RKM solutions are difficult to see due to good agreement with the CE solution.The RMS deviation is also shown.Here, q ref. uses N l = 6 results since more Legendre components included in the expansion should improve the accuracy of the solver.
while converging to the physical solution.The two hyperparameters of most interest are the number of terms kept in the Legendre expansion N l and the minimum velocity cutoff ξ min = w min /v th .
The RKM solver is run for an isochoric singlespecies setup similar to the isochoric case in 4.1 for a variety of values of N l with N K = 1/1000, 1/20, 1/10 to investigate convergence in different cases of nonlocality.The root-mean-square (RMS) deviation ) 2 is introduced to measure convergence of the heat flux, where D is the spatial domain and |D| is the length of the domain.The reference heat flux q ref. is chosen from a solution of the RKM solver with identical parameters as q apart from the hyperparameter of focus.
Figure 10 shows the convergence study for N l .Good agreement with the local analytic solution is practically identical in all cases as expected, since the local solution for δf (l) is only nonzero for l = 1.In non-local cases, good convergence is demonstrated even for N l = 4. Larger N l increases convergence time, in particular scaling with O(N l ).The default value N l = 6 is chosen.
For the low velocity cutoff, a value around ξ min ≈ 1 should be taken since a smaller value includes the bulk in the RKM solver domain, which gives rise to large computational cost and numerical instabilities, and a larger value does not account for parts of the tail that have significant contribution to the heat flux.
In the region 0.9 ≤ ξ min ≤ 1.1, there is a variation of σ q < 0.1% for N K ≤ 1/10, therefore the heat flux is mostly insensitive to this hyperparameter for a choice which appropriately defines the tail.This further validates the approach of solving only for the tail of the distribution beyond some cutoff since the RKM and local Chapman-Enskog solutions match smoothly at the low velocity boundary.The default value ξ min = 1 is chosen.

Temperature self-consistency
The RKM solver takes input number density and temperature profiles and solves for the tail of the DF of each species.In extreme non-local cases, a sufficiently non-Maxwellian tail may result in the density and temperature moments of the resulting DF to significantly differ from the input profiles, leaving an inconsistent solution.Since the temperature is a higher moment than the density, the former is more sensitive to deviation due to non-local tail behaviour and more likely to be inconsistent.The RKM temperature is defined as ) Figure 11 shows the RKM temperature for the previously discussed single-species isochoric case, DT ∇T case, and CH 2 ∇T case.In all three cases, the local cases are in good agreement with the local result, with inconsistency developing with increasing Knudsen number.Even for the extreme cases of N K = 0.1, the total temperature is consistent to ≲ 5% deviation.In the multi-species cases, it is the temperature of the lighter species that is more inconsistent with the input profile.

Smallness of field-particle operator
Figure 12 shows a comparison of the l = 1 components of the test-particle operator and field-particle operator at the hot end (z/L = −1.5)and cold end (z/L = +1.5)from the single-species isochoric example in 4.1 for the non-local N K = 0.1 case against local analytic results.The vertical orange dashed lines indicates the local low velocity cutoff, chosen to be the local thermal velocity w min = v th .The solution for the tail above the low velocity cutoff w ≥ w min is determined by the RKM solver.The solution below the low velocity cutoff in the bulk is fixed as the local Chapman-Enskog solution δf (w < w min ) = δf local (w), therefore the test-particle operator is also identical for the RKM and local results away from the cutoff.In contrast, the RKM and analytic solutions of the fieldparticle operator may differ in this region due to the    velocity integrals that sample the solution above the low velocity cutoff.This deviation is small due to the small population of the tail compared to the bulk, therefore the local analytic result should still hold as a good approximation of the bulk.
At the velocity grid point just below the low velocity cutoff there is a spike in the testparticle operator compared to the local solution.This originates from the w −2 ∂ w (w 4 ν ∥ ∂ w δf ) term, suggesting that the numerical scheme as implemented does not constrain continuity of the second velocity derivative of the DF continuous at the cutoff.This numerical artifact lies outside the domain of the RKM solver, and continuity of the second or higher derivatives of the DF are not required, therefore is not a point of concern.
The red highlight indicates the regions around w/v th,local = 5/2 ≈ 1.6 where the local analytic result δf (1) local ∝ (w 2 /v 2 th − 5/2) changes sign.The sudden spike of the ratio in this region is therefore not of concern since both the operators are passing through zero.Apart from these red regions, the field-particle operator remains small in comparison to the test-particle operator.The Jacobian of the fieldparticle operator is then a small difference between the approximate and true Jacobian, which explains why the approach of approximating the Jacobian of the collision operator in the tail without the non-sparse elements of the field-particle operator described in 2.4 works numerically.

Discussion
A reduced kinetic method has been presented and applied to investigate non-local multi-species ion heat transport for the first time.The RKM solves for the tail of distribution functions according to the Vlasov-Fokker-Planck equation for a classical, unmagnetized multi-species plasma in a 1D2V geometry.In the single ion species case, the non-local ion heat flux show qualitatively similar trends to its well explored electron counterpart.On the other hand, novel non-local features are recovered for the multi-species ion heat transport due to additional thermodynamic drives such as concentration gradients and diffusion of different species in the centre-of-momentum frame.
Results for DFs and heat fluxes in local cases reproduce the analytic results of the Chapman-Enskog expansion [16,17,19] for multi-species ion heat flux and Braginskii's [15] results for ion and electron heat flux in a simple plasma.Non-local electron heat flux results agree well with an SNB implementation, particularly where the field-particle operator is fixed at the local solution rather than implemented in a first principle fashion.Expected qualitative behaviour of non-local effects are demonstrated, in particular suppression of the peak heat flux and enhancement of heat flux away from regions of steep gradient.In multispecies cases, non-local effects are demonstrated from gradients in temperature, concentration, and pressure, with some cases giving rise to large heat fluxes far away from the boundary where the local result is exponentially small.The multi-species heat flux has been decomposed into two components, the reduced heat flux and enthalpy flux, where the former is the term employed in mainline hydrodynamic codes and the latter is an inherently multi-species feature due to relative dynamics of the species in the center-ofmomentum frame.As demonstrated here for the first time, in cases with pressure or concentration gradients, the enthalpy flux is particularly insensitve to nonlocal effects, thus becoming increasingly important in comparison to the reduced heat flux for high Knudsen numbers.The RKM has been applied here to describe non-local heat flux, but since the approach solves for the DF of each ion species, it may have future utility in investigating other kinetic effects such as reactivity reduction in inertial confinement fusion experiments [24,29,47].Future work relating to the heat flux may investigate additional effects such as magnetic fields and the implementation of a RKM prescription to nonlocal heat flux in a simple 1D hydrodynamics code.
Writing 3 in 1D2V coordinates and neglecting the hydrodynamic velocity assuming u ≪ v th yields the kinetic equation for δf α as

Field-particle operator
The Legendre components of the field-particle operator [36] are

Figure 1 .
Figure 1.Visualization of 1D2V geometry at a boundary.There is an azimuthal symmetry in the velocity coordinate ϕ, and translational symmetry parallel to the boundary.

Figure 2 .
Figure 2. Visualization of solver domain for the set of PDEs for δf (l) α .

Figure 3 .
Figure 3. Reduced heat flux q (rescaled by the Knudsen number, N K ) for a single species isobaric (a) and isochoric (b) temperature step for a variety of Knudsen numbers from the RKM solver.The local analytic result for the reduced collision operator and Braginskii's exact result for the full Coulomb collision operator are included.The enthalpy flux has not been included since it is identically zero for a single species, and the temperature and number density profiles are shown below on the same horizontal scale for reference.
(a) f (0) in hot end (b) f (0) in cold end (c) f(1) in hot end

Figure 4 .
Figure 4. DF comparison at the hot end z = −1.5Land cold end z = +1.5Lfor various Knudsen numbers against local thermal speeds.Note the plots of f (1) have a symmetric log scale to include negative values.The smaller subplots show the ratio of each graph to the local Maxwellian with an appropriate scaling against Knudsen number included.

Figure 6 .
Figure 6.Components of heat flux for a DT mixture subject to an isolated temperature gradient (a), isolated pressure gradient (b), and isolated concentration gradient (c).The shape of the relevant profiles are shown at the bottom of each subfigure, with undisplayed profiles being fixed constant at unity or concentrations at 0.5.Results for N K ≤ 0.01 are not visible due to agreement with the local result.

Figure 7 .Figure 7 .
Figure 7. Components of heat flux for a CH 2 mixture subject to an isolated temperature gradient (a), isolated pressure gradient (b), and isolated concentration gradient (c).The shape of the relevant profiles are shown at the bottom of each subfigure, with undisplayed pressure or temperature profiles being fixed constant at unity or concentrations such that n H /n C = 2.

Figure 9 .
Figure 9.Comparison of electron heat flux from the RKM solver, an SNB approach, and a flux limiter approach for an isochoric simple hydrogen plasma with a Gaussian temperature profile.(a) uses the full linearized Coulomb collision operator, and (b) uses the Coulomb collision operator with the fieldparticle part fixed at the local result from the Chapman-Enskog solution.

Figure 10 .
Figure 10.Convergence study for the number of Legendre modes included in the expansion δf ≈

Figure 11 .
Figure 11.RKM temperature for the single-species isochoric case (a), DT ∇T case (b), and CH 2 ∇T case (c).The upper subplots show the RKM integrated temperature, and the lower subplots shows the ratio of the RKM temperature to the input temperature profile.For the multi-species cases, the lighter species (D or H) temperature is shown in dashed lines, the heavier species (T or C) temperature is shown in dotted lines, and the total ion temperature T RKM total = α nαT RKM α / α nα is shown in solid lines.

Figure 12 .
Figure 12.Comparison of l = 1 components of the test-particle operator and field-particle operator at the hot end (a) and cold end (b) from the single-species isochoric example.The solid lines show results for the N K = 0.1 results, and the dashed lines show the local analytic results.The lower subplots show the ratio of field-particle operator to test-particle operator.

( 5 )F
Substituting into the field-particle operator, we have C

( 6 )
[31]he 'test-particle' operator which describes collisions of species α against a Maxwellian background of species β.The second term C F αβ {δf β } = C αβ {f M α , δf β } is the 'field-particle' operator which describes the effect of species β on the Maxwellian part of species α.It can be argued that the third, nonlinear term C αβ {δf α , δf β } is subleading to the testparticle term |C αβ {δf α , δf β }| ≪ |C T αβ {δf α }| even for order-unity deviations from Maxwellian in the tail of the distribution, since δf β (v ′ ) is smaller than f M β (v ′ ) globally in v ′ space, and the second argument of C αβ is integrated over v ′[31].Hence, the full kinetic equation is reduced to a linear problem with the collision operator