The fastVFP code for solution of the Vlasov–Fokker–Planck equation

We describe the fastVFP code for solution of the Vlasov–Fokker–Planck equation for non-local electron transport and the generation of magnetic field, especially for application to laser-produced plasmas. We describe the essential features of the code that make it fast and robust and suitable for inclusion as a transport package in a fluid simulation. We present a few sample results that demonstrate the abilities of the code.


Introduction
Electron energy transport plays a crucial role in laser-plasma experiments, most notably in spherical implosion experiments aimed at the realisation of inertial fusion energy (IFE) as a commercial source of electricity.In the direct drive approach to IFE (Craxton et al 2015) the laser propagates through the low density plasma surrounding the capsule and deposits its energy at densities less than the critical density.Laser energy is absorbed by the plasma electrons which diffusively conduct energy to higher density to generate the high pressure that drives the implosion.
Diffusive heat flow is also important as a process that smooths non-uniformities in laser energy deposition and increases the uniformity of pressure.Thermal smoothing is particularly important for IFE where the maintenance of spherical symmetry is essential in the face of the Rayleigh-Taylor instability that exponentially increases the amplitude of nonuniformities.
For gentle temperature gradients in an unmagnetised plasma the electron heat flow is described by the Spitzer theory (Spitzer and Härm 1953).If the plasma supports a magnetic field, electron transport is described by the more complicated (Braginskii 1965) theory that describes the resulting anisotropic linear transport.
The Spitzer and Braginskii models are accurate when the scale length of the temperature profile L ∼ T/|∇T is very much larger than the mean free path λ t of electrons with the thermal velocity v t = √ eT/m e .The linear Spitzer theory breaks down when L < 100λ because the heat flow is carried by electrons on the high velocity tail, v ∼ 3-4v t , of the electron distribution function f.These tail electrons have a mean free path of the the order of 100λ t since the mean free path is proportional to the fourth power of the velocity.
If L < 100λ t the heat carrying electrons are not collisionally confined to high temperature regions of plasma.They are free to escape and become spread relatively uniformly across the whole plasma.The spatial gradients of the heat-carrying electrons are consequently reduced and the heat flow reduced below expectations from the Spitzer theory.This is known as non-local heat flow.Until this was understood (Bell et al 1981), the reduction in heat flow was named 'inhibited transport' and modelled in simulation codes by a 'flux limiter' F that limits the maximum heat flow to FQ free where Q free = n e eT(eT/m e ) 1/2 is the free-streaming heat flow (Malone et al 1975, Gray andKilkenny 1980).
The physics of non-local transport is described by the Vlasov-Fokker-Planck (VFP) equation, but solution of the VFP equation has proved computationally expensive in more than one spatial dimension on distance and time scales relevant to IFE implosions.Consequently, many hydrodynamic simulation codes continue to use a flux-limited Spitzer model for electron transport.
Various approximate transport models have been developed that improve upon flux-limited Spitzer.Some calculate the smoothed heat flow as a spatial integral over the energy input and temperature profile at distances up to 100λ t (Luciani et al 1985).More recently, the SNB (Schurtz, Nicola, Busquet;Schurtz et al 2000) model has successfully implemented a reduced form of the VFP equation.Sherlock et al (2017), Ma et al (2022) compared the SNB method with their K2 code that solves the full VFP equation.
Approximate models function adequately when the temperature profiles are essentially one-dimensional (1D) and the plasma is unmagnetised, but they are less able to capture accurately the heat flow when it is multi-dimensional or when there is a magnetic field or when non-locality produces strong anisotropy in the electron distribution function making the plasma susceptible to the Weibel instability.
There remains considerable scope for improvements in computational efficiency and in the range of physics accessible with VFP codes.VFP codes can take many forms according to the physics that dominates, for example when the Langdon (1980) effect is important in collisional absorption (Shaffer et al 2023, Turnbull et al 2023).Our aim here is to construct a code (fastVFP) that can be used as a reliable, robust, and sufficiently accurate workhorse electron transport package in sophisticated IFE simulations.Carefully chosen approximations are made and strategies adopted to ensure robust fast operation with accurate representation of a wide range of multi-dimensional physics.Because of its increased efficiency and the inclusion of magnetic field and arbitrary anisotropy, fastVFP is able to explore phenomena on hydro time and distance scales that are beyond easy reach with previous VFP codes.
The improvement in speed and robustness in fastVFP mainly arises from our treatment of the zero velocity boundary where we fix the distribution function to its Maxwellian value with suitable adjustments.This provides a stable anchor for more complicated anisotropic non-Maxwellian transport processes at higher electron velocity.
The temperature T is in eV and SI units are used unless otherwise stated.

The fastVFP code
The VFP equation for a single species of non-relativistic particles with charge q and mass m is ) where f is the distribution function and the term on the right hand side models collisions.The VFP equation is supplemented by the Maxwell equations which calculate the electric and magnetic fields from the current density j of charged particles We follow the KALOS method described by Bell et al (2006) which expands the distribution function as a sum of spherical harmonics: where f −m l = (f m l ) * θ = 0 in the x direction and ϕ = 0 in the y direction.The VFP equation in this form can be applied in three spatial directions with three components of E and B but for simplicity we limit ourselves to two spatial dimension x, y and a single component B z out of the plane of simulation.Correspondingly there are two components E x and E y of the electric field.
The method described can be extended to full dimensionality if required.
In laser-plasmas energy is carried mostly by electrons.The electrons must be treated kinetically by solution of the VFP equation.In most cases, the ions can be treated as a fluid.Here we set out the equations for stationary ions that are uniform in density with a charge Ze.The ions therefore play no part in the calculation except that they elastically scatter electrons.In laser-plasmas ions move characteristically at the sound velocity which is much smaller than the electron thermal velocity.If required, ion motion is best included by defining the electron distribution in the local ion rest frame and adding small terms to the VFP equation.Electron angular scattering by ions is Z times larger than electron-electron collisions but energy exchange between electrons is crucial in determining the electron distribution function.Angular scattering of electrons by electrons can crudely be included in the collision term by increasing the electron-ion collision term by a factor between 1 and 1 + 0.5/Z.This is adequate for many laserplasma experiments since there are similar uncertainties in the Coulomb logarithm and the state of ionisation that determines Z (Le et al 2019).
The full electron-electron collision term is complicated and requires the calculation of the Rosenbluth potentials (MacDonald et al 1957, Lyman 1965) for every spherical harmonic interacting with every other spherical harmonic (Tzoufras et al 2011).The anisotropic terms in the spherical harmonic expansion at high electron velocity are responsible for the transport processes, but their scattering is dominated by collisions with the bulk of the electrons that are at thermal velocities and mainly isotropic.Hence the most important Rosenbluth potential is that formed by integration over f 0 0 .The Rosenbluth potentials formed by integration over all other f m l have relatively little effect.
In the current form of fastVFP we ignore the effect of Rosenbluth potentials formed by integrating over the anisotropic parts of the distribution.That is, in line with many previous VFP codes (e.g.Bell et al 1981, 2006, Matte and Virmont 1982, Kingham and Bell 2004), we ignore Rosenbluth potentials calculated by integrating over any f m l where l > 0. As explained below, we go further by also using an approximate form for the Rosenbluth potential calculated by integrating over f 0 0 Next steps in improved accuracy would be (i) to use an accurate Rosenbluth potential instead of our approximation for the f 0 0 potential, and (ii) to include the effect of Rosenbluth potentials calculated by integrating over f 0 1 and f 1 1 .The most significant improvement in many cases would be to calculate the Rosenbluth potentials calculated by integrating over f m 1 and apply these to f 0 0 in the collision term in the equation for the evolution of f m 1 .This limitation to the f 0 0 Rosenbluth potential has been applied to many VFP codes going back to Bell et al (1981).The omission of the f m 1 Rosenbluth potentials reduces the Spitzer thermal conductivity by around 30% as we discover when reproducing Spitzer with fastVFP in the linear regime with Z = 4 (see figure 3 in section 5).This reduction could be compensated for by artificially decreasing the Coulomb logarithm but we prefer not to do this because in the non-local regime it may introduce other errors that we could not easily quantify.In any case, non-local effects are probably more important in practice than accurately reproducing the Spitzer conductivity.
For the purposes of this paper we ignore the complexity of calculating the full Rosenbluth potentials by simplifying the electron-electron collision term to where T is the local electron temperature.ν 0 = 8.2 × 10 5 n e log Λ which is equivalent to ν 0 = 3.76(eT/m e ) 3/2 /τ NRL where τ NRL is the electron-electron collision time in the NRL plasma formulary (2018).This approximation works well for the electrons on the tail of the distribution that are mobile and carry the energy.At velocities less than the electron thermal velocity, this term overestimates the the effect of collisions, but this matters little since the approximate collision term forces the distribution towards a Maxwellian which would be the case if the full Rosenbluth potentials were implemented.The effect of using this approximation to the Rosenbluth potentials is explored in section 6 below.If required, any of these approximations to the collision term could easily be removed, but at the cost of making the code slower to execute and more complicated to implement.
In a further approximation we omit the displacement current from the Maxwell equations.We have included it as an explicit term in some versions of fastVFP and confirm that it makes very little difference in transport problems.
In the limited dimensionality of two spatial dimensions and retaining only one component of magnetic field, the VFP equations with these approximations can be separated into a series of equations for each m and l using the spherical harmonic expansion for the distribution function (Bell et al 2006): ) where v x = v cos θ and v y = v sin θ cos ϕ.In this dimensionality, all the coefficients f m l are real and evolve according to the equations ν 0 is a constant that takes the same value at every velocity and everywhere in space when the plasma density is uniform.The magnetic field is updated at each timestep from The magnetic field evolves slowly, so the electric field can be updated from the magnetic field at the previous timestep.However, the calculation of the electric field is performed implicitly to accurately produce the required electron current, j = ∇ × B/µ 0 and avoid spurious electron drifts that might accumulate over long times.The electric field is chosen such that The VFP equation for each f m l is updated at each timestep by separating it into a velocity part and a spatial part.The spatial differentials are represented by a spatially centred first order difference across adjoining cells.This works well in the context of the solution of the full VFP equation where the order l + 1 spherical harmonics act to diffuse the order l harmonics in the way that f 0 1 and f 1 1 act to diffuse f 0 0 .The velocity part is solved implicitly because of the short collision times at low electron velocity.Initially we experimented with solving the spatial part implicitly using an iterative alternating direction implicit (ADI) method.ADI worked well for simple problems but the number of iterations became excessive in more complicated problems and we found an explicit method to be more robust and result in faster execution of the code despite the restrictive Courant condition on the timestep imposed by the highest electron velocity.
The number of spherical harmonics can be chosen according to the problem under consideration.In problems for which the tensor term f m 2 plays a significant role, such as for the Weibel instability, we typically extend to l max = 7, with some runs extending to l max = 15 to show that going beyond 7th order makes no difference.We apply a small amount of spatial smoothing to the highest order spherical harmonic to mimic the diffusive effect of including the term at l max + 1 in a nondisruptive manner.Without the additional smoothing there are rare occasions when structures build up at l = l max and then propagate down the orders to cause problems at lower orders in l.
It has been shown that a smaller number of harmonics is sufficient for some problems (Matte and Virmont 1982), and some codes take the expansion to only first or second order (Kingham andBell 2004, Thomas et al 2009).The required computing resource is roughly proportional to the square of the maximum order, which raises the prospect that a systematic examination of the required order may substantially increase the efficiency of the code, with particular benefits if fastVFP is coupled to a fluid code for large scale experimental simulation.
The above approximations are mainly implemented for convenience and they could easily be removed.However, there is one further approximation that is essential to the scheme and makes it more robust and able to execute quickly.This concerns the boundary condition on the distribution function at zero velocity.The correct boundary conditions at zero velocity are that f m l (v = 0) = 0 for all l > 0 and that ∂f 0 0 /∂v = 0.The difficulty with the ∂/∂v boundary condition on f 0 0 is that there is no correction on f 0 0 (v = 0) if it floats away from its correct value.In the end, there is little to stop errors accumulating and f 0 0 (v = 0) becoming unphysical and even going negative, in which case the calculation fails.A possible alternative strategy would be to artificially change the electric field to draw electrons into the hole in the velocity distribution, but this affects the wider distribution function at all velocities and other problems develop.Our answer to the problem is to choose a boundary condition that sets f 0 0 (v = 0) to what it would be if the distribution function were Maxwellian.If necessary we make minor adjustments to the temperature (T ′ instead of T) and sometimes the density (n ′ e instead of n e ), to correct for distortions to the distribution function, for example due to high velocity electrons escaping and depleting the tail of the distribution.So, our boundary condition at zero velocity is In the specimen calculations presented in section 5, T ′ is updated at each time at each point in space from the equation where T f = 3m e 2e ˆ∞ 0 4π f 0 0 v 4 dv and τ ee (T) is the thermal electron collision time.At the end of the reference run (figure 1) the maximum values of T and T f are 2346 eV and 2325 eV respectively indicating that this procedure works well.At present we set n ′ e = n e .At the end of the reference run, ´∞ 0 4π f 0 0 v 2 dv lies within 4% of the specified n e at all points on the grid.The deviation from the correct density is in part due to discretisation errors when integrating over the velocity grid when the electron temperature is low.
A further aspect of our procedure for the boundary condition is that the temperature T is calculated from macroscopic equations instead of by taking moments of the distribution function.At each timestep, fastVFP integrates over f 0 1 and f 1 1 to calculate the heat flow Q.The energy density at each point is then updated and the temperature determined by where σ is the energy input.The macroscopic update of the temperature adds resilience to the code.

Energy input
As described above, energy input is determined through the above macroscopic equation for the evolution of the electron energy density U e .The energy input is controlled by the quantity σ that is imposed externally on the solution of the VFP equation.Energy input in the VFP equation works as follows.
The choice of a locally positive value of σ increases the local energy density and temperature.This increases T from what it would have been without the energy input.The effect of the increased T in the approximate VFP collision term enhances the collisional transfer of energy to higher electron velocities where heat flow carries the energy away to the colder parts of the plasma.In this process, energy input essentially takes place at low electron velocity.This recipe for energy input is a reasonable approximation when energy is absorbed by inverse bremsstrahlung into the more collisional low velocity electrons.It would not be appropriate if energy were absorbed by collisionless processes such as resonance absorption-in this case a tailored absorption term would need to be added to the VFP equation.The current version of fastVFP would need to be improved if it were being used to investigate the details of absorption (see for example, Shaffer et al 2023, Turnbull et al 2023) or the atomic coupling to radiation and ionisation (Le et al 2019).However, if electron transport is the focus of the investigation and inverse bremsstrahlung absorption dominates, the current version of fastVFP should be adequate since it completes the essential processes of absorbing the energy at low velocity and transferring it by electron-electron collisions to high velocity where transport occurs.
Another way of understanding absorption in the current implementation of fastVFP is that the increase in T produced by energy input is the equivalent of adding an extra term to the electron-electron collision term so that it becomes ) .
This models collisional absorption through the addition of an additional velocity diffusion coefficient D vv /v 3 that is proportional to the velocity dependent collision rate (∝ v −3 ).Adding the term in D vv is the equivalent of increasing T to effect the imposed energy input through σ.As may be verified by integration in velocity, D vv is related to σ by which for a Maxwellian distribution becomes

The electric field
All quantities (f m l , B z , E x , E y ) are calculated at cell centres.The magnetic field is updated at each timestep from the curl of the electric field.This poses a computational problem since it requires a first difference across three cells if ∇ × E is calculated from E at the cell centres, possibly leading to chequerboard errors and instabilities.To circumvent this difficulty we also calculate the electric field at the cell vertices from the distribution function averaged from the cell centres.
In each case the electric field is calculated by advancing the velocity part of the VFP equation for f m 1 three times, first for E x = 0 and E y = 0, second for E x = E c and E y = 0, and third for E x = 0 and E y = E c , where E c is some chosen characteristic value of the electric field.In each of the three cases, the currents j x and j y are calculated and the results of the three calculations combined linearly with a weighting chosen to give j = ∇ × B/µ 0 .The weighted combination of the three solutions determines the electric field.This procedure guarantees ∇.j = 0 exactly.The procedure is the equivalent of the method previously used by Tzoufras et al (2013) as defined in the application of their equation ( 19).The computational overhead of repeating the calculation three times for different electric fields is small since only f 0 1 and f 1 1 need be calculated.It is important to set j = ∇ × B/µ 0 , and not j = 0, since the non-zero j is needed for the diffusion and resistive decay of magnetic field.

Some specimen results
In this section we demonstrate the operation of the code.Output from a reference run of fastVFP is shown in figure 1.The spatial grid is 100 × 100 µm with a grid spacing of 0.5 µm.The timestep is 1.33 fs.There are 100 uniformly spaced velocity grid-points with the maximum grid velocity set to ten times the thermal velocity of an electron at temperature of 2 keV.The expansion in spherical harmonics is taken to l = 7.The boundary conditions are reflective which means that equivalent mirror images of the energy input are present on the far side of each boundary.
The plasma is initially uniform with an electron density of n e = 10 22 cm −3 and a temperature of 200 eV.Z = 4.The Coulomb logarithm is fixed at logΛ = 2. Constant heating is applied for 40 ps in a cluster of three spots with the spatial profile given in figure 1.This configuration may be relevant to experiments such as those designed to probe magnetic reconnection (for example Willingale et al 2010) where the plasma is deliberately irradiated with a small number of laser spots.
Transport in such a configuration may also be relevant to smoothing the effects of laser speckle in IFE experiments.Laser speckle occurs mostly on the smaller scale (typically a few µm) of the laser wavelength mutiplied by the f -number of the focussing lens, but cross-beam energy transfer between incident and reflected waves can result in structure on larger scales (Follett et al 2017) which is less easily smoothed.Further results are shown below that are more realistic for small scale speckle.
In our two-dimensional geometry it is difficult to define the energy input in terms of an equivalent laser intensity.A meaningful measure of energy input and an equivalent laser intensity is given by the heat flow supported by the plasma as shown in figure 1.The maximum heat flow in figure 1 is 7 × 10 14 Wcm −2 which can be taken as a guide to an equivalent laser intensity.
After 40 ps, figure 1 shows the growth of a substantial magnetic field exceeding 100 T (1 MG).The heat flow is filamented, probably by the collisional Weibel instability (Weibel 1959, Epperlein and Bell 1985, 1987, Thomas et al 2009).In this variant of the Weibel instability, the anisotropy is driven in the tail of the distribution function by the heat-flow.For example, if a temperature gradient is directed along x−, then heat-carrying electrons that move down the gradient enhance the pressure anisotropy P xx , and this effect is captured by the inclusion of f m 2 .Although other mechanisms can generate magnetic field in the non-local regime (Kingham and Bell 2006), we have verified that the magnetic field in our simulations is dominated by pressure anisotropy.
The bottom-right plot in figure 1 relates the local heat flow Q to the local free-streaming limit Q f = n e eT(eT/m e ) 1/2 for each cell on the spatial grid.Each point on the plot is colour coded according to the local temperature as plotted in figure 1. L/λ on the horizontal axis is the ratio of the local temperature gradient T/|∇T| to the local thermal electron scattering mean free path λ.If the Spitzer theory applied, all points would lie on the line labelled 'Spitzer'.The plot shows that the heat flow is everywhere smaller than the Spitzer heat flow with an equivalent 'flux limit' of 0.03-0.1Qf .The spatial distribution of Q/Q f is also plotted separately in figure 1.
The calculation shown in figure 1 took 6 h on 400 processors of the SCARF cluster at the STFC Rutherford Appleton Laboratory.No systematic attempt was made to optimise the code.The value of ∆t = 1.33 fs for the timestep is determined by the Courant condition for electrons at maximum velocity on the velocity grid.A less restrictive timestep is probably acceptable for the evolution in velocity space since it is treated implicitly.Sub-cycling the spatial calculation with a larger timestep in the velocity calculation may give a gain in computational efficiency.Also, ∆t might reasonably be varied according to plasma conditions.Similarly, the maximum order of the spherical harmonic could be varied according to plasma conditions.In similar calculations, we found very little difference when extending the calculation to order 15 (taking approximately four times longer to run).Similar considerations apply to the spacing and maximum value of the velocity grid.The code is parallelised by domain decomposition into blocks of cells in x, y.The execution time decreases inversely with the number of processors down to 4 × 4 blocks of cells on each core.
For comparison figure 2 shows results from a calculation with exactly the same parameters as for figure 1 except that input energy is reduced by 1/3.The heat flow and the temperature is correspondingly reduced.The lower heat flow is not sufficient to excite the Weibel instability.There is no filamentation and the magnetic field is 30 times lower than in figure 1.
Figure 3 shows results after 10 ps when the spatial box is enlarged to 1 cm × 1 cm so that the temperature gradients are small and non-local effects are weak.The grid spacing is 50 µm and the timestep is 6.67 fs.The heat input profile still consists of three spots, but now they overlap and their radius is increased to 3 mm.The maximum heat flow is 5 × 10 12 Wcm −2 .In the colder parts of the plasma (the blue dots in the bottom right plot in figure 3) the heat flow is 2/3 times the Spitzer flow.The departure from Spitzer is caused by the approximation to the electron-collision term, especially the omission of the Rosenbluth potentials due to integration over f 0 1 and f 1 1 .As noted above, this could be corrected at the cost  of additional computation, but the benefit would be minimal since there are similar or greater uncertainties in the Coulomb logarithm and in the effective Z due to partial ionisation.In the high temperature part of the plasma (the red dots) the larger departure from Spitzer is due to non-local effects which are present even at such small temperature gradients and small values of Q/Q f .The magnetic field plotted in figure 3, reaching 10 −5 T is a numerical artefact.It occurs at the edge of the heating spots, and it would be zero due to symmetry except that the square spatial grid introduces a small discretisation error.
Figure 4 presents results from a calculation more directly relevant to laser speckle.Energy is input in a collection of spots, each with a half-height diameter of 8 µm.The spot centres are randomly spread over a circle of diameter 60 µm.The grid spacing is 0.5 µm and the maximum order of the spherical harmonic expansion is 7.The maximum heat flow is 8 × 10 14 Wcm −2 and the temperature reaches 2500 eV at the end of the calculation lasting 50 ps.The magnetic field reaches 0.8 MG.The heat flow filaments as it flows out of larger circle, but is relatively unaffected between the spots.
The calculation producing the results shown in figure 5 is the same as for those shown in figure 4 except that the magnetic field is turned off.The difference can mainly be seen in the lack of filamentation.The maximum heat flow is similar in both cases, but in figure 4 it is focussed into filaments.The maximum temperature is also very similar indicating that the filamentation has little overall effect on the escape of energy from the hot plasma.
The differences between figures 4 and 5, with and without magnetic field, are clearer in figure 6 which plots the temperature in each case.

An assessment of the approximate collision term
In this section we investigate the effect of adopting the approximate collision term as described in section 2: ) .
We compare this approximate collision term with the more complete collision term for an isotropic electron distribution set out in equations ( 2)-(4) by Bell (1985) which can be written in the form ( ∂f ∂t At high velocities where transport occurs δ T and δ n are both small in keeping with our approximate collision term.At low velocities both δ T and δ n are close to one and it can easily be shown that for a Maxwellian distribution both δ T and δ n converge to the same limit. 1 + δ T and 1 + δ n → 1 n e ˆv 0 4π f (u) u 2 du which retains the correct ratio (1 + δ T )/(1 + δ n ) = 1 as required for a Maxwellian to be the equilibrium distribution.The neglect of δ T and δ n in fastVFP over-rates the strength of both the diffusive and advective parts of the collision term, but crucially does not much change the ratio (1 + δ T )/(1 + δ n )  which ensures that the low velocity part of the distribution is Maxwellian.
We assess the effect of neglecting δ T and δ n by calculating what δ T and δ n would be for a simple fastVFP simulation of energy input in a single spot.The output from fastVFP is displayed in figure 7 after 10 ps.The electron distribution produced by fastVFP is output at the positions of the three white crosses at 5, 15 and 30 µm from the centre of the heated region.The distributions are used to calculate δ T and δ n as shown in figure 8.The magnetic field is switched off for this calculation.Figure 8 shows the depletion of the tail of the distribution function in the hot plasma and its corresponding enhancement in the cold plasma.
Figure 8 confirms that δ T and δ n are small where the heat flow is strong, so their neglect is generally acceptable when calculating the heat flow.(1 + δ T )/(1 + δ n ) stays close to one and thereby performs the important task of keeping the distribution close to Maxwellian at low velocity where transport is negligible.The deviation of (1 + δ T )/(1 + δ n ) from 1 at low velocity within the heated region at 5 µm is related to the excess velocity diffusion needed to transfer energy to high velocities from the low velocity electrons.Instead of  augmenting the diffusion coefficient by including δ T in the calculation, fastVFP instead artificially increases the temperature from T to T ′ as described in section 2, which has a similar effect.
In the cold plasma, from 15 to 30 µm, the energetic electrons escaping the heated region are scattered predominantly by electron-ion collisions and slowed down to join the cold plasma by the energy loss part of the electronelectron collision term.The energy loss rate follows a v −3 dependence (|δ n | ≪ 1), as shown in figure 8, except at low velocities where the balance against diffusion is more important through the maintenance of the correct ratio (1 + δ T )/(1 + δ n ).The dip in the distribution function at v/v ref = 1.3 and 30 µm is related to ohmic heating of the return current.
The acceptability of the δ T = δ n = 0 approximation depends on the questions that the VFP solver is being required to address and its function as a physics package in a larger simulation.If the issue is that of calculating heat flow in cases similar to the examples chosen above or providing a basis for a magnetic field calculation, the approximation is probably adequate and has the benefit of being computationally fast and robust.If the details of energy absorption, such as the impact of the Langdon (1980) effect or energy input directly into hot electrons, are important, then a more complicated electron-electron collision term is needed.The numerical formulation adopted here could easily be adapted to include non-zero δ T and δ n , or a specific heating term, if required.

Conclusions
A few straightforward approximations, especially in the boundary condition at zero velocity, are sufficient to facilitate robust and flexible solution of the VFP equation in two dimensions with arbitrary anisotropy and self-consistent magnetic field.Some of the approximations in fastVFP are made for computational convenience and could be removed if needed.The present code should be fast enough to be made part of a hydrodynamic simulation of IFE experiments and there are ways in which it could be optimised for increased efficiency.The extension to arbitrary anisotropy enables modelling of physical phenomena such as the Weibel instability that are absent in the diffusive ('f 0 + f 1 ') approximation.The inclusion of the tensor term f m 2 is important in two dimensions where it is crucial in the lateral spreading or focussing of fluxes carried by f m 1 .Our aim in this paper is to describe the fastVFP code rather than systematically investigate the physics of electron transport in laser-produced plasmas.Nevertheless, our specimen calculations point to matters for further investigation, notably the role of magnetic field and the Weibel instability in experiments.Even without magnetic field, non-local transport may have an important effect on the response of a plasma to laser speckle (Thomas et al 2009, Sherlock andMichel 2022).It appears that inertial fusion experiments sit in a threshold regime in which the Weibel instability may or may not be important.We caution against drawing conclusions from the few results we present.The full problem is three-dimensional.Fluid flow and the Nernst effect may easily advect magnetic field and temperature non-uniformities out of the region in which they are generated.We hope that fastVFP may offer a basis for more complete modelling, either in three dimensions, or in carefully designed two-dimensional simulations, coupled to full fluid simulation and backed up by comparison with VFP codes that make fewer or different approximations.The specimen results presented here appear to indicate that magnetic field may have little effect on the overall heat flow, thereby having little effect on the overall hydrodynamics of inertial fusion capsules.The main effect of magnetic field may be to seed small scale nonuniformities through heat flow filamentation.Further investigation may be warranted, especially of the timescale over which speckles must fluctuate in order to avoid harmful imprint.

Figure 2 .
Figure 2. Same as figure 1 except the input energy is reduced by 1/3.

Figure 3 .
Figure 3.Results for small temperature gradients after 10 ps.The maximum heat flow is 5 × 10 12 Wcm −2 .The computational grid is 1 cm × 1 cm with a cell size of 50 µm.

Figure 4 .
Figure 4.The effect after 50 ps of energy input in small spots, each with a diameter of 8 µm, representing laser speckle.The spots are arranged in a circle of diameter 60 µm.The full grid size is 250 µm × 250 µm.The maximum heat flow is 8 × 10 14 Wcm −2 .

Figure 5 .
Figure 5.The same as figure 4 but with the magnetic field switched off.

Figure 6 .
Figure 6.Comparison in 3D format of the temperature profiles in figures 4 and 5, with and without magnetic field.

Figure 7 .
Figure 7. 30 µm × 30 µm grid.Data shown in the next figure for the three crosses at distances 5 µm, 15 µm and 30 µm from the centre of the energy input spot.Data after 10 ps.

Figure 8 .
Figure 8. Data shown for the three crosses in figure 7 at distances 5 µm, 15 µm and 30 µm from the centre of the energy input spot.The velocity axis is the same in each plot.v ref is the thermal velocity at a temperature of 2 keV.The collision frequency for energy loss is plotted in the bottom-right graph with v −3 plotted (dashed line) for comparison.