Contemporary particle-in-cell approach to laser-plasma modelling

Particle-in-cell (PIC) methods have a long history in the study of laser-plasma interactions. Early electromagnetic codes used the Yee staggered grid for field variables combined with a leapfrog EM-field update and the Boris algorithm for particle pushing. The general properties of such schemes are well documented. Modern PIC codes tend to add to these high-order shape functions for particles, Poisson preserving field updates, collisions, ionisation, a hybrid scheme for solid density and high-field QED effects. In addition to these physics packages, the increase in computing power now allows simulations with real mass ratios, full 3D dynamics and multi-speckle interaction. This paper presents a review of the core algorithms used in current laser-plasma specific PIC codes. Also reported are estimates of self-heating rates, convergence of collisional routines and test of ionisation models which are not readily available elsewhere. Having reviewed the status of PIC algorithms we present a summary of recent applications of such codes in laser-plasma physics, concentrating on SRS, short-pulse laser-solid interactions, fast-electron transport, and QED effects.


Introduction
Particle-in-cell (PIC) codes for simulating collisionless plasma kinetics are amongst the most widely used computational tools in plasma physics. With algorithmic development since the 1970s they are also one of the most studied plasma simulation techniques. The core of PIC algorithms have been extensively summarised in, for example, [1] and [2], and earlier editions of these texts. These early texts and papers covered the stability, self-heating and accuracy of PIC schemes in great detail. However, since these landmark publications modern PIC codes have moved on to using higher order shape functions routinely, coupling Maxwell's equations to particle dynamics through the flux of charge, rather than directly calculating moments of distribution functions and using filters which preserve the solution of Poisson's equation. In addition Contemporary particle-in-cell approach to laser-plasma modelling PIC codes are now routinely extended to include collisions, ionisation and QED effects. It seems timely therefore to review what is now considered modern PIC and include updates to self-heating rates, convergence tests and test problems which if they are published elsewhere at all are scattered and often overlooked by PIC users. The first part of this paper therefore acts as a review of both standard collisionless PIC algorithms and their more recent physics extensions. The second section of the paper highlights some recent examples of the use of PIC codes in laser-plasma physics. These demonstrate what is now possible with the combination of new algorithms and high-performance computing (HPC) and how PIC codes are entering a realm of increased realism with full ion-electron mass ratios and 3D dynamics.
Section 2 deals with a review of standard collisionless PIC algorithms. Section 3 explains how these are extended to include collisions and ionisation and section 4 how QED effects can be included for high-field physics. A full implementation of these algorithms exists in many PIC codes used in laser-plasma physics. This papers uses the implementation of these algorithms in the EPOCH code developed by the authors. EPOCH was tested to estimate self-heating rates, the convergence of the collisional algorithm to Braginskii transport, tests of the ionisation model and convergence to the continuum QED limit. Self-heating rates have been available for a long time but these are rarely reported for PIC codes with high-order splines and Poisson preserving current density accumulation with smoothing. These are presented here in a convenient form appropriate for laser-plasma studies. Similarly, collisional PIC has become relatively commonplace, but systematic convergence tests are mostly lacking. The EPOCH code will give the same results as all codes based on this algorithm set and is chosen here simply as it was developed by the authors and is freely available for all to test via the UK CCP-Plasma website [3].
Key to any credible use of a PIC code is a detailed understanding of the accuracy and convergence of the core scheme and any additional physics modules, e.g. ionisation or QED effects. In section 5 we present such analysis for the core collisionless numerical self-heating, collisional routine convergence to near equilibrium transport coefficients and tests of the ionisation and QED routines. Each of these physics modules is tested in isolation so, for example, the self-heating rate will be different when collisions are on compared to the collisionless case reported here. This is not just due to the physical heating from collisions, but also that collisions change the statistics of numerical heating.
The tests below are intended as a guide only. When using PIC codes beyond the collisionless limit, or with multiple additional physical effects on simultaneously, PIC code users will need to follow similar techniques to test the accuracy their simulations.
Finally, in section 6 we review recent high-profile PIC simulations which demonstrate what is currently possible with the combination of the algorithms presented here and current HPC facilities. PIC codes are now so embedded in plasma physics research that a comprehensive review of all PIC code results relevant to laser-plasma physics is impractical. Instead, we aim to focus on reviewing publications which both elucidate exciting new physics and help to define the limits of what is technically possible with PIC codes at present. This section also discusses the need for hybrid PIC models for laser-solid interactions.

Core collisionless algorithm
The core of a PIC code is two coupled solvers: one that moves charged particles freely in space under the influence of EM fields and calculates the currents due to the particle motions (the particle pusher) and another that solves Maxwell's equations on a fixed spatial grid subject to the currents calculated from the particle motions (the field solver). Between these two solvers the full collisionless behaviour of a kinetic plasma can be simulated. The core algorithm is standard for PIC codes and is briefly summarised here for completeness.

The finite-difference time-domain method
The finite-difference time-domain method (FDTD) is a standard technique for solving Maxwell's equations numerically. The E and B fields are specified on a Yee staggered grid [4] (shown in figure 1). This grid staggering means that centred, second order accurate derivatives are easily implemented. For example, ∂ E x y is defined as ⎛ where Δx is the distance between cells in the x-direction. The resulting derivative is second order accurate at the location of B zi j k , , . Since ∂ E x y is only used to update B z , and all other derivatives are similarly defined at the correct point in space all derivatives are second order accurate. This staggering also means that the scheme conserves ∇ ⋅ B.
The FDTD scheme implemented in most PIC codes, including EPOCH, uses a modified version of the leapfrog scheme in which the field is updated at both the full time-step and the half time-step. The FDTD scheme used in EPOCH takes the following form: first the fields are advanced one half time-step from n to n + 1/2, using currents calculated at n: (3) At this stage, current is updated to + J n 1 by the particle pusher (as described in section 2.2), and the fields are updated from n + 1/2 to n + 1 to complete the step: By re-writing equations (2) and (5) as an update from n − 1/2 to n + 1/2, i.e. shift n by −1 in equation (5), it is clear that the FDTD update is exactly the same as the leapfrog scheme: The splitting used in FDTD is simply to allow all fields to be defined at the same time and calculate half-step values needed for the particle push. This also makes including additional physics packages, which often require EM fields at either full or half-step values, easier.

Particle pusher
The particle pusher solves the relativistic equation of motion under the Lorentz force for each particle in the simulation. In order to calculate the particle trajectory to second order accuracy the electric and magnetic fields at the half time-step are used after they are calculated in the first half of the Maxwell solver: where α p is the particle momentum, α q the particle's charge, α x the particle position and α v the particle velocity. The particle velocity can be calculated directly from the particle momentum using γ = . Most PIC codes use the Boris rotation algorithm [5] which splits the equation (8) into separate parts responsible for the acceleration of the particle in the E field and the rotation of the particle about the magnetic field.
In the standard leapfrog method, the next step is to update the particle position using the following second-order update: The currents needed for the Maxwell solver can then be calculated from the particles using the mechanism of Esirkepov [6] which is a generalisation of the Villasenor and Buneman [7] current deposition scheme. This scheme has the property that the electric field update resulting from the currents calculated in this way, rather than from moments of the distribution, always satisfies ρ ϵ ∇ ⋅ = E / 0 , where ρ is the charge density. One disadvantage of the leapfrog method when written in this form is that the particle positions and momenta are not both known at the same time. To overcome this issue, we split equation (9) into two parts: We then discard α + x n 3/2 after completing the current update. However, this means that before the update of momentum we only have the particle position defined at α x n so we require an extra update beforehand using the equation: n n n 1 2 (12) The combination of equations (10) and (12) yields a second-order update from α x n to α + x n 1 and the whole process is exactly equivalent to that given using equation (9).

Shape functions
To keep the problem computationally tractable each simulation particle normally represents a large number of real particles, with the number of real particles represented by each simulation particle called the weight. When simulations contain changes in density there is the option to either represent this by changing the weight of each simulation particle to reproduce the density structure or by keeping the weight of the simulation particles constant and changing the number density of simulation particles to match the real density.
To calculate the force on a particle, the E and B fields must be known at the particle location rather than on the fixed spatial grid. Similarly the current has to be deposited at the grid locations to update the E field. Since each simulation particle contains many real particles it is necessary to choose a spatial distribution of particle weighting throughout the volume occupied by a simulation particle. The simplest method is to distribute the particle evenly throughout a volume Δ × Δ × Δ x y z. This is commonly referred to as the 'top hat' shape function.
Since any function which has unit integral and compact basis is suitable for use as a shape function, other higher order shape functions are often used such as triangular shape functions with volume of Δ × Δ × Δ x y z 2 2 2 and splines with a valume of Δ × Δ × Δ x y z 4 4 4 . Formally all of these shape functions are b-splines but we use the terms top-hat, triangle and spline for convenience throughout the paper. Note that the triangle shape function is the convolution of the top-hat function with itself and the spline the convolution of the top-hat with itself four times. A spline over three cells is not used as the one sided nature requires the same computational effort as that over four cells for less accuracy. Note that the shape functions defined above refer to the physical shape of the macro-particles. Interpolation from these shapes to find grid quantities, and the inverse of finding field quantities at particle centres, requires particle weight functions which are the convolution of the shape function with the top-hat function. This nomenclature differs from that used in some textbooks [1].

Binary collision operator
The PIC method is a robust and reliable approach to the kinetic modelling of plasmas which are dominated by collective effects. However, these codes generally neglect particle interactions over very short (less than grid scale) ranges. At high temperatures (≳1 keV) and relatively low densities (≲10 27 m −3 ) collisional effects in plasmas are generally considered minimal (i.e. the mean time between collisions is comparable to the time scales of interest), and the collisionless approximation used in PIC codes is valid. However, at lower temperatures and/or higher densities the effect of sub-grid scale interactions on the evolution of the system can become non-negligible.
As a means of accounting for discrete particle collisions, many PIC codes make use of some form of collision algorithm [8][9][10][11][12][13][14][15] which stochastically scatters particles in momentum space. The majority of algorithms in use are derived from the binary collision approach in [8] which performs Rutherford scattering on pairs of particles which reside in the same grid cell. Alternative methods using, for example, an additional force term in the Lorentz force calculation (a 'collision field') [9] have also been developed.
A binary collision algorithm, based on the approach of Sentoku and Kemp [14] has been implemented in EPOCH. A summary of this approach is presented below and test cases later in the paper. In order to simplify momentum conservation, collisions are calculated in the centre-of-momentum reference frame of the two particles. Lorentz transformations are required in order to evaluate the particles' momenta in the centre-of-momentum frame, and ensure the collision algorithm is fully relativistic.
It is assumed, without loss of generality, that a particle i of species α is being scattered off of a particle j of species β (with α β = a possibility). The collision frequency is given by: is the reduced mass and v r is the closing velocity of i and j. Scattering angles for particle i in the centre-of-momentum frame are calculated as: where ∈ [− ) Q 1, 1 , ∈ [ ) R 0, 1 and ∈ [ ) S 0, 1 are random numbers. The angle θ π π ∈ [− ) , corresponds to the angle within the scattering plane relative to particle i's initial momentum vector. The other angle, ϕ, is the scattering angle in the plane perpendicular to particle i's initial momentum vector. The form of equation (14) ensures that the distribution of scattering angles within the collision plane obeys θ ν ⟨ ( )⟩ = Δt tan /2 2 [14] in order to favour small-angle scattering. However, the second scattering angle, ϕ π ∈ [ ) 0, , is randomly selected from a uniform distribution. The definitions of the scattering angles are also demonstrated in figure 2.
The post-collision momentum of particle i in the centre-ofmomentum frame can subsequently be calculated: p e e e cos sin cos sin sin , and since momentum must be conserved, correspond to unit vectors defining an orthonormal basis set, with ê 1 parallel to particle i's initial momentum vector.
To simplify particle pairing, the main per-processor particle list is split into an array of secondary lists, one per species per grid cell. For intra-species collisions (when α β = ) collisions are performed by scattering each odd numbered particle in the list off of the next particle in the list. For an odd number of particles, the last three particles in the list collide with each other (two collisions per particle), but using a collision frequency at half of its normally calculated value. Essentially the last three particles each undergo two half-collisions. Thus all particles within the cell have undergone one full collision. When performing inter-species collisions, all particles from the more numerous species undergo one and only one Diagram to illustrate the scattering angles, θ and ϕ, that particle i (red) is scattered through when colliding with particle j (blue) in the centre-of-momentum frame.
collision. While running through the more numerous list, the code loops over the second species list to form the successive collision pairs. The result is that some (or all) of the particles from the second species may undergo multiple collisions. This is corrected for by the application of a common collisional time interval, which also corrects for the effect of non-uniform particle weights. The collision frequency is multiplied through by a correction factor F [11]: where the sum over k is over all collision pairs for species α and β. Thus while equation (16) defines a single frequency for each pair of particles in a binary collision event these are modified, and hence potentially different for the two particles, due to variable weights and different particle numbers per cell. Also the secondary particle list is randomised after all particles in that list have been collided once.

Ionisation
Ionisation has been shown to have some important macroscopic consequences in laser-plasma interactions such as defocussing of the laser [16], injection of electrons in specific regions of the plasma [17] and fast shuttering, the process by which a reflective critical surface is formed in plasma mirrors [18]. EPOCH includes a number of ionisation models to account for the different modes by which electrons ionise in both the field of an intense laser and through collisions. How ionisation is included in PIC codes is often not explicitly explained so here we present the approach used in EPOCH along with test cases later in the paper. In 1965, Keldysh derived formulae describing field ionisation for a hydrogen atom in the low frequency regime where photon energy is beneath the binding energy of the electron. This introduced the Keldysh parameter γ which is used to separate field ionisation into the multi-photon and tunnelling regimes. In Hartree atomic units the Keldysh parameter is given by [19]: where ω is the photon frequency, m e the electron mass, ϵ the ionisation energy for the electron, e the electron charge, and E the magnitude of the electric field at the electron. Multi-photon ionisation ( γ ≫ 1) can occur when an electron absorbs a photon that does not have enough energy to cause ionisation or excitation to a higher energy state. In this case electrons can be excited to virtual energy states and it is possible for electrons to absorb further photons and ionise before the virtual states decay. EPOCH models multi-photon ionisation with a semi-empirical WKB approximation, outlined below, based on that of Ammosov et al [20]. Tunnelling ionisation (γ ≪ 1) considers the deformation of the atomic Coulomb potential by the imposed electric field which can create a finite potential energy barrier through which electrons may tunnel. This is modelled using the ADK ionisation rate equation [21] averaged over all possible magnetic quantum numbers. In a special case of tunnelling ionisation the potential energy barrier can be lowered beneath the electron binding energy such that it may escape classically; this is known as barrier-suppression ionisation (BSI) and is included by a correction to the ADK ionisation rate by Posthumus et al [22] with the threshold electric field used for ADK cut-off being the point at which the self-consistent field strength equals the atomic field strength ϵ = E Z / T 2 [20]. Ionisation events are tested for every particle with a bound electron at every time step; the model used for field ionisation is based on the self-consistent electric field strength calculated at the particle to be ionised. In [23] it is noted that γ < 0.5 is a reasonable transition point from multi-photon to tunnelling models which corresponds to a transitional field strength of To produce a monotonically increasing ionisation rate W, where E B is the turning point in the ADK rate equation such that ( ) , we use the following test to select a model: ionise ; this method may be repeated by adding the sampled ionisation time from U 2 to t ionise until no ionisation occurs or > Δ t t ionise . When no further ionisation occurs the particle is removed from the system and replaced with the ion associated with the multiple ionisations and an electron with the appropriate weight for releasing multiple electrons per ion; these are added to the system at the original particle position and velocity. Energy conservation is accounted for by a current density correction through Poynting's theorem; in tunnelling ionisation and BSI the energy loss from the field is the ionisation energy of the electron, in multi-photon ionisation it is the total energy for the number of photons absorbed. The total energy loss from multiple ionisations ϵ t is summed and a current density correction is weighted back to the grid points: An additional energy correction is also applied to electrons released through multi-photon ionisation. The energy, as measured in the ion rest frame, of such an electron is ω ϵ − K , where K is the numer of photons absorbed. This energy is account for by adding ω ϵ ( − ) m K 2 e to the electron momentum in the direction anti-parallel to the local E-field.
Where the plasma must be considered collisional, EPOCH includes an extension to the collision module for ionisation by an electron impacting upon and exciting a bound electron. Electron impact ionisation cross-sectional data is unavailable for many elements and orbitals and so we make use of the MBELL equation [24] based on analytical fits of extensive empirical data by Bell et al [25]. This is valid for ⩽ q 36 where = − q A N Z n l for A Z , the atomic number, and N nl the total electrons in all suborbitals up to the ionising orbital. For q > 36 we approximate the cross-section using the relativistic modified binary encounter Bethe model (MRBEB) [26]. Ionisation is handled by including a preionisation step before collisions. When species are being collided we test whether one is an ionisable species and the other is an electron.
As before we use inverse transformation to test for ionisation; for each colliding pair we test for where n e and v e are the electron density and velocity respectively and F is the weighting correction from equation (17) for the cross-section σ. It is noted that the model used to produce the ionisation rate or cross-section for field or collisional ionisation can be easily changed within the module. Upon an ionisation event both colliding particles are tagged in the code and after all pairings are tested the particles to be ionised are split into ion and electron particles. The ions are added to the relevant species list whilst the ionising impact electron and the ionised target electron are separated into new lists. The collisions routine is called on these electron lists whilst the ions are left unscattered. After an impact ionisation, but before Coulomb collisional scattering, the ionised target electron is assumed to be at rest with respect to the target ion. Unionised particles are unaffected by the preionise routine and undergo normal scattering.
Energy loss due to ionisation is accounted for via a reduction of the incident electron kinetic energy. This energy reduction is performed in the rest frame of the target ion before being transformed back into the simulation frame to find the final kinetic energy so as to account for relativistic effects. This transformation is simplified by first rotating the simulation frame such that the ion moves parallel to one of the new rotated coordinates. When an ionisation occurs, the momentum reduction is calculated in the ion rest frame by subtracting from the incident electron kinetic energy prior to scattering as:

Partial superparticle ionisation
A potential problem with the Monte-Carlo ionisation scheme described in section 3.2 is that a particle with a high weight, i.e. one which represents more real particles than other macro-particles, can only either be neutral or fully ionised and this transition may be too sudden in some simulations.
An alternative scheme, described here for Hydrogen but easily generalised, is to allow a neutral particle to split into a smaller, i.e. lower weight, neutral plus an ion and electron. The weights are then suitably adjusted so that the total mass and charge are conserved. However, this splitting is not performed for every step as this may lead to a large number of very low weight electrons and ions. Instead a minimum particle weight is M min defined, then on each time-step the effective weight of any macro-particle which corresponds to unionised plasma is calculated. After k time-steps this unionised weight, M k is determined from: In this way the cumulative history of fractional ionisation for any particle is maintained. Once ⩾ M M k min the particle is split into a neutral with weight − M M k 0 , and an ion and electron pair each with weight M k . A comparison of this scheme to the Monte-Carlo scheme in section 3.2, now called the whole ionisation scheme to distinguish from the partial ionisation scheme, is shown in figure 3.
This scheme exhibits the correct ionisation statistics as demonstrated in figure 3 and allows for capturing very small amounts of ionisation that would otherwise require a greater number of superparticles to observe. That is to say that ionisation events when − (− Δ ) < W t N 1 exp 1/ can still be analysed. This comes at the cost of using more memory per superparticle due to requiring as many additional doubleprecision floating point numbers per superparticle as there are ionisation levels for the species. Whilst use of the whole superparticle ionisation scheme shows poor resolution when − (− Δ ) < W t N 1 exp 1/ , if viewed over n time-steps such that − (− Δ ) > Wn t N 1 exp 1/ then W is still correctly recovered. The partial superparticle ionisation scheme is best used for simulations where a low amount of ionisation is expected or if the region where ionisation occurs is small compared to the simulation domain; in this situation the computational cost can be significantly reduced.

QED effects
Next generation 10PW-class lasers, such as several of those comprising the Extreme Light Infrastructure, are predicted to generate electromagnetic fields so strong that nonlinear quantum electrodynamics (QED) emission processes play a crucial role in the plasma dynamics [27][28][29]. The parameter determining the importance of QED effects for an ultra-relativistic electron in a strong electromagnetic field is: Here E RF is the electric field in the rest-frame of the electron and = ℏ= × E m c e / 1.3 10 s e 2 3 1 8 Vm −1 is the critical field for QED-the Schwinger field ( μν F is the electromagnetic field tensor and ν p is the electron's 4-momentum). QED effects become important when η ≳ 0.1 [30]. This can be parametrised in terms of the laser intensity I as η ∼ ( × I 0.1 /5 10 22 Wcm −2 ). QED effects therefore become important at ≳ × I 5 10 22 Wcm −2 , well within the range of 10PW-class lasers. The model used for including QED processes in EPOCH is described in detail in [31] and will be summarised here.
The inclusion of QED emission processes in a PIC code is dramatically simplified by relying on the large separation of scales and splitting the electromagnetic fields into low-frequency and high-frequency components. The low-frequency component is due to the electromagnetic field of the laser (and collective plasma processes). The laser field is a coherent state [32] and therefore the low-frequency fields may be treated classically. The high-frequency component, henceforth referred to as 'gamma-ray photons', is that emitted by the electrons on acceleration by the laser fields. This emission is incoherent so emissions may be treated as independent and the emitted gamma-ray photons as particles following null geodesics. The interaction of electrons and positrons with the electromagnetic fields is treated using the quasi-classical model of Baier and Katkov [33], whereby electrons and positrons follow classical trajectories in the laser and plasma (lowfrequency) fields between discrete emission events. Emission is described in the 'furry' representation [34]: the basis states are 'dressed' by the laser/plasma fields. We then make a perturbation expansion and keep the lowest order Feynman diagrams: single gamma-ray photon emission by an electron or positron; pair production by a gamma-ray photon.

The important effects: photon emission and pair production
The determination of the rate of photon emission and pair production is simplified by making the following two assumptions on the low-frequency laser and plasma fields: (i) the formation length of the photons and pairs is much less than the scale of variation of the field. In this case the low-frequency fields can be treated as constant during the emission process. The ratio of the photon formation length to the laser wavelength is 1/a 0 in the field of a plane electromagnetic wave (where ω = a eE m c / 0 0 e L , E 0 is the peak electric field of the laser and ω L its frequency) [35]. The QED processes are usually only important for ≫ a 1 0 so this assumption is reasonable in most cases. (ii) The lowfrequency plasma and laser fields are much weaker than the Schwinger field. At the peak intensity likely to reached on next generation laser facilities (∼10 24 Wcm −2 ) the electric field is only ∼10 −3 times the Schwinger field. In this case the emission rates are approximately independent of the field invariants and depend only on the invariants η (for photon emission) and (for pair production). Emission in any general low-frequency field configuration is then approximately the same as that in any other field configuration with the same η or χ with small F and G. We choose emission in a constant magnetic field where photon emission is referred to as synchrotron emission and pair production as magnetic pair production [36]. The rate of synchrotron emission by an ultra-relativistic electron is: (25) where the emitting electron has Lorentz factor γ.
The rate of magnetic pair production from a gamma-ray photon of energy ν γ h is: 16 2/ 3 / 1/3 2 (K 1/3 is the Bessel function of the second kind).

Including QED effects using a Monte-Carlo algorithm
An efficient way of capturing the stochasticity of the quantum emission processes is by using a Monte-Carlo algorithm (a) The partial superparticle ionisation scheme correctly reproduces the ionisation rate even when − (− Δ ) < W t N 1 exp 1/ , but (b) it can be shown that under the whole superparticle ionisation scheme the correct ionisation rate can still be recovered over enough time-steps. [37,38]. The cumulative probability that a particle emits after traversing a region of plasma with optical depth τ em is = − τ − P 1 e em . The point at which a given particle emits is decided by assigning it a value for P at random between 0 and 1. The equation for P above is then inverted to yield τ em , the optical depth of plasma it traverses before emission occurs. For each particle the optical depth evolves according to where N t d /d is the appropriate emission rate as given in equation (25) or equation (26) above. When τ τ = em emission occurs. In the case of photon emission the emitting electron or positron recoils, conserving momentum. The change in momentum of the emitting electron or positron is set equal to ν −( )γ h c p / i , p i is its initial momentum and ν γ h the energy of the emitted photon. This recoil gives the quantum equivalent of the radiation reaction force [39]. The probability that the emitted photon, emitted by an electron or positron with energy parametrised by η, has energy parametrised by a given , / and the cumulative probability that it has a given energy is . Each emitted photon is assigned a value of χ P randomly between 0 and 1. The value of χ to which this corresponds is determined by linear interpolation from tabulated values for χ P . For pair production the pair-producing gamma-ray photon is annihilated and the constituents of the pair share its energy. The probability that a fraction f is taken by one particle in the pair , after creation from a photon with energy parametrised by χ, is given in [40]. For each emitted pair f (which we arbitrarily define as the fraction of energy taken by the electron) is determined in the same way as the energy of emitted photons. The cumulative probability is assigned a random value between 0 and 1, which is then used to obtain f . The emitted photons and pairs are added to the simulation as additional macro-particles.

Additional numerical constraints
The inclusion of the above QED processes in a PIC code adds additional numerical constraints. It is important that the timestep be small enough that the emission is well resolved. The constraint is that Δ ≪ Δ t t QED where: is the inverse of the maximum possible rate of photon production. As well as this time-step constraint it is also important that sufficient macro-particles be used to resolve the emission. If the standard deviation in emitted energy due to quantum fluctuations is σ, then the standard deviation in the energy emitted by N macro-particles is σ σ = N / N . We require that this be much less than the total emitted energy if the emission is to be well resolved.
The possibility of pair cascades-where an emitted photon generates a pair, that emits further photons, which go on to generate further pairs-means that the number of macroparticles in the simulation can grow exponentially. This can necessitate the use of particle merging if the number of macroparticles becomes infeasibly large [41].
The method for calculating radiation reaction, i.e. electron and positron recoil, outlined above, does not conserve energy. The relative error in the energy for an electron or positron with initial Lorentz factor γ i and final Lorentz factor Only ultra-relativistic particles will have sufficiently high η (as to have a reasonable probability of emission thus γ i is always much larger than unity and the error in energy conservation is small. The method for sharing the energy of the gammaray photon between the electron and positron in the pair also introduces an error in energy conservation. For a photon of energy ϵ ν = γ γ h m c / e 2 generating an electron and positron with Lorentz factors γ − & γ + respectively, the relative error is: Only photons with energy much larger than m c e 2 will generate pairs so this error is also small.

Stability and self-heating
It is well known that PIC codes are prone to a phenomenon known as self-heating. This is a stochastic heating which, possibly after an initial thermalisation stage, leads to linear heating of the plasma. This has been studied in detail in [42]. However, the analysis presented there concentrates primarily on the case in which particle forces are assigned to nearest-neighbour gridpoints (NGP) and is also only applicable to momentum conserving PIC codes and not charge conserving schemes. A full parameter study of self-heating for NGP and Top-hat shape functions for electrostatic problems showed improvements when moving to larger and smoother shape functions [43,44] although again these tests were not based on Poisson preserving current accumulation [6] or smoothing [45].
Current smoothing can also have a significant impact on self-heating and a series of test cases have been presented in [46]. EPOCH has both high-order shape functions and current smoothing although care must be taken in order to ensure that charge is still conserved. A simple method for achieving this is to weight the current over neighbouring cells such that the total sum over all cells is preserved. The implementation of chargeconserving current smoothing used in EPOCH is based on that presented in [45]. Figure 1 of the paper [46] presents results obtained using the OSIRIS PIC code for a periodic box containing a plasma with initial density of × 1.11 10 23 cm −3 (100 times the critical density of the plasma in a 1 μm laser field) and initial temperature of 1 keV. The number of grid points is kept fixed and the length of the domain varied in order to alter the ratio between Debye length, λ D , and cell width, Δx. The number of particles per cell is also varied. figure 4 illustrates the same set of test cases performed using the EPOCH code. The results obtained are broadly similar to those presented in [46] and are shown here using the same scales for comparison.
From [43,44] the estimated heating rate for a Top-hat shape function, once the heating has become linear in time, can be written as: here T 0 and λ D 0 are the temperature and Debye length defined at a time once the heating is linear in time. We have assumed that the Debye length λ D 0 is not resolved by the grid spacing Δx, as is common in laser-solid simulations. n ppc is the number of computational particles per cell. This can be re-arranged into a form more useful for estimates in laser-solid interactions: where T eV is the temperature in electron volts, t ps is time in pico-seconds, n 23 is the plasma electron number density in units of 10 23 cm −3 and Δx nm is the cell width measured in nanometres. Using this formulation the results in [43] for a Top-hat shape function correspond to α ≃ × 5 10 were not used in these averages as these numbers were an order of magnitude larger than The reason for these anomalous results is due to the rapid selfheating, before the linear heating stage, which is observed only for Top-hat shapes (figure 5). In these cases, the heating is so rapid that the temperature becomes relativistic before linear heating is reached. Clearly Top-hat shape functions should never be used for laser-solid simulations. The averages are therefore only quoted to the first significant figure and are meant purely as a rough guide to self-heating for laser-solid interactions. Nonetheless, using these numbers, along with equation (30), gives an order of magnitude estimate for the self-heating to be expected with EPOCH or similar PIC codes.
Finally, we show a detailed comparison of heating curves for a single case in which each particle shape function is tested along with current smoothing. Here it can be seen that by far    the greatest reduction in heating effects comes via the use of high-order particle shape functions (figure 7).

Thermalisation of non-equilibrium distribution via selfcollision
A basic requirement for any collision operator is that it asymptotically approaches a unique equilibrium distribution in a homogeneous plasma. That is to say that it must satisfy the H-theorem. A simple test of the efficacy of a collision operator is, therefore, its ability to relax from a non-equilibrium distribution to a Maxwellian distribution via self-collisions. We initialise the system with a waterbag type distribution, whereby the the distribution of particles in = ( ) p p p p , ,  [47], where v te is the electron thermal velocity. In these tests we use = × − p 6.04 10 c 24 kg m s −1 (equivalent in energy to a 50 eV thermal distribution); a single species of electrons at a uniform density of = × n 1 10 e 30 m −3 ; 4 cells in x; periodic boundaries; a total system length equal to the Debye length; and a Coulomb log of (Λ) = log 4. The electric field solver was disabled in order to eliminate the possibility of numerical collisions, although we find their impact to be negligible in the collisionless case under equivalent conditions. Figure 8 shows the electron energy distribution at t = 0 and = t t c , along with the equilibrium Maxwellian distribution for runs conducted with 1024 and 8192 particles per cell. In all cases the distribution rapidly converges toward a Maxwellian distribution with some residual particle noise, particularly in the tail.
We measure the convergence of the electron distribution using an L1 norm on the difference between the distribution and the target thermal distribution; defined as: where f 0 represents the equilibrium, or target, thermal distribution and f represents the particle distribution at time t. Both distribution functions are calculated on a discrete energy grid with N 1 eV bins and normalised to the total electron number. Figure 9 shows the convergence in norm of the electron distribution over time for various particle numbers. This convergence is initially dominated by the rapid convergence of the lower energy components, and is consistent with the results of reference [48]. The asymptotic value is determined by the particle noise inherent in any PIC model, as demonstrated by the scaling with particle number shown in figure 9. Further tests (omitted for brevity) were performed for an initially mono-energetic distribution with equivalent energy, in which the particles are initialised in a shell in momentum space at | | = p p / 5 c . These demonstrated similarly rapid convergence to the target thermal distribution.
These results are a demonstration of the ability of the collision operator implemented in most PIC codes to model self-collisions, as well as a confirmation of the relative insignificance of numerical collisions in the case of systems with fully self-consistent fields.

Relaxation of an anisotropic particle distribution
In order to test the angular scattering properties of the collision operator, we perform a test on temperature isotropisation. An initially anisotropic temperature distribution will be isot- Here v te is the thermal velocity calculated using the final equilibrium temperature of T = 40 eV. The lowest energy band shows little angular dependence, while the higher energy band shows increasing anisotropy with the maximum around θ = 0. In the highest energy band, electrons only exist around the central maximum in the v x direction. As time progresses the anisotropy reduces. We expect this reduction to take place on the time scale dictated by the collision frequency of the electron-ion collisions. Because the collision frequency behaves like v −3 with the relative velocity of the colliding particles, we expect the isotropisation of the higher energy bands to proceed more slowly.
The right panel of figure 10 shows the decrease in the difference between the maximum θ After about τ 0.05 c the distribution has isotropised. This is the timescale expected for electrons around v 0.36 te . The relatively large noise after this time is caused by the normalisation and the small initial anisotropy in this energy band. The next energy band, < < v v 0.5 / 1.5 te , reaches the isotropic state after about one collision time, τ c . This is the timescale expected for the energy band around the thermal velocity. The next higher energy band, < < v v 1.5 / 2.5 te , again isotropises slightly faster than expected on a timescale of approximately τ 5 c . The central velocity, = v v 2 te , would suggest a timescale of τ 8 c . This is due to the fact that the distribution function falls off rapidly with v in this energy range. Therefore, the isotropisation is mainly due to the slower electron with higher collision frequency. This trend is also seen in the highest energy band, > v v / 2.5 te , although full isotropy has not been reached by the end of the simulation.
In conclusion, the simulations of the isotropisation of an anisotropic Maxwellian distribution produces the expected results. Isotropisation for the bulk of the distribution takes place on timescales comparable to the collision time and the v 3 scaling has been demonstrated for other energy bands.

Spitzer resistivity
The current drawn by a constant, uniform electric field applied to an infinite plasma provides a means of testing the ability of a collision algorithm to reproduce the Spitzer conductivity, σ τ = q n m / 2 e e . This is itself a prerequisite for modelling the electron transport and Ohmic return current heating which dominate the transport processes in high intensity laser-plasma interaction. The current density for such a system will increase with time, asymptotically approaching the In all cases the distribution rapidly converges to a Maxwellian distribution. Energy is normalised to the average energy (equivalent to a 50 eV Maxwellian distribution), the particle distributions were generated using 1 eV energy bins. steady-state 'Spitzer-limited current' given by Ohm's law: σ = J E. The variation of the current density with time can be obtained by seeking a time-dependent solution to the Drude model for electron transport: For this test, the system consisted of a periodic 1D box containing a uniform hydrogen plasma of density = n 10 e 29 m −3 and temperature T = 300 eV. Fields ranging between × 5 10 8 V m −1 and 10 10 V m −1 were applied across the system. In order to ensure that the field was constant, and that stopping was only due to collisions, the field updates were disabled. A fixed Coulomb log of (Λ) = log 10 was used throughout. The current densities, as functions of time, from simulations performed with × 5 10 4 macro-particles (50% electrons, 50% ions) are shown in the left-hand image of figure 11, and are in reasonable agreement with equation (32) (dashed lines). The current densities, averaged over the last 10 fs, are plotted as a function of applied field in the right-hand image of figure 11.
Performing a linear fit to these results allows for an estimate of the conductivity. This was found to be × Ω − 9.09 10 6 1 m −1 , which is comparable to the analytically expected value of × Ω − 9.98 10 6 1 m −1 [47].

Thermal conduction
The use of periodic boundaries and averaging of the current density over the entire system in the Spitzer conductivity test, tends to increase the particle statistics. In contrast, the thermal conduction test detailed in reference [49] requires accurate modelling of particle collisions on an individual grid cell basis, and as such is a more demanding test of a collision algorithm's accuracy. A uniform density plasma is initialised with a temperature profile which increases from 100 eV to 400 eV over λ 200 D (where λ D is the electron Debye length). A fixed Coulomb log of (Λ) = log 10 was used throughout. The results for a density of   For some values of L, the heat flux is double-valued. This is due to the temperature profile being such that each value of temperature gradient occurs twice. The lower band of heat flux values correspond to the region near the high temperature region, and exhibit flux limiting as observed in reference [49]. The upper band correspond to the cooler region, and are artificially inflated by high energy particles which have moved through from the high temperature region; hence exceeding the Spitzer-Harm heat flux (solid line).
Additional simulations were also performed with lower particle statistics. Figure 13 shows the same data as figure 12, but with 500 particles per cell instead of 50 000. Comparing the two results demonstrates that in order to accurately model electron transport in high density plasmas, very large numbers of particles are required in each grid cell. This is one application where the weaknesses of the PIC model, and the strengths of the continuum Vlasov-Fokker-Planck (VFP) approach [50][51][52], are most apparent. Figure 11. Left: current densities as a function of time resulting from applied electric fields of × 5 10 8 (blue), 10 9 (green), × 2 10 9 (orange), × 5 10 9 (red) and 10 10 V m −1 (black). Right: late time current densities (averaged over the last 10 fs), as a function of applied field (crosses). Squares denote equivalent values for the Drude model (equation (32)). The error bars indicate the standard deviation in the simulation data over the averaging period. The dashed line represents a linear fit to the simulation data, and was used to estimate the plasma conductivity.

Ionisation injection for laser wakefield acceleration
The use of a laser to accelerate electrons in the wakefield of a plasma wave was first suggested by Tajima and Dawson [53]. This a novel alternative to conventional radiofrequency particle accelerators that may serve to reduce the cost and space requirements of laboratory electron accelerators. The mechanism by which wakefield acceleration occurs is now well understood and verified by numerous experimental results [54][55][56][57]; modern investigation focusses upon how to increase the population of electrons being accelerated by this mechanism.
Umstadter et al suggested selecting the gas and laser intensity used for laser wakefield acceleration (LWFA) such that field ionisation of tightly bound electrons would occur in such a way that the ionised electrons would be added exclusively to the accelerated electron bunch [58]. McGuffey et al presented the first experimental results of this ionisation injection as a means for enhancing the high energy electron population for LWFA [17]. They used a neutral helium gas mixed with 1-5% additives of various high-Z gases and laser intensities in the bubble regime [56] for low electron plasma densities <1% n c . McGuffey demonstrates an order of magnitude increase in high energy electron density with 1% nitrogen added to helium for a 0.8 nm, 30 fs laser pulse of focussed intensity × 3 10 19 W cm −2 and spot size 10 μm. We reproduce these results here to demonstrate the field ionisation module.
The wakefield is a fairly stable structure and is readily produced in 2D PIC simulations; initially this was produced with EPOCH in a preformed plasma. The domain was 64 μm wide, with reflecting boundaries top and bottom to simulate a gas capillary; it was 64 μm long with a moving window that would follow the laser pulse for 0.75 mm. The laser had a wavelength of 800 nm focussed to a 10 μm spot size and a Gaussian temporal profile with a 30 fs full width at half maximum. The domain was divided into × 1024 1024 cells; it was found that any less than 8 cells per wavelength caused the wakefield structure to break up after a relatively short propagation distance. The capillary was uniformly filled with a hydrogen plasma using 64 particles per cell and a density approximately 1% that of critical density. The wakefield was seen to form after a relatively short distance and remain stable for the full propagation distance with a population of electrons accelerated into the rear of the bubble. As reported in previous works, the relativistic laser intensity sees the electrons almost entirely vacating the wake of the laser, and the electric fields within the bubble are large and mostly longitudinal.
In McGuffey's results, the 1% nitrogen case showed as much as a three orders of magnitude difference in electron density within the relevant region of the energy distribution compared to the pure helium case; this is illustrated in figure 14. The mixed gas was preionised in all cases up to the inner s-shell of nitrogen; so to match this, the helium was fully preionised from a neutral gas density of 10 19 cm −3 , whilst the nitrogen preionised to the fifth ionisation state (i.e. + N 5 ). The + N 5 was added at 1% the density of the helium and the preionised electron plasma density was therefore increased by 5%. Preionising the gas in this way relaxed the requirements on the number of ion superparticles, as the majority of the electron superparticles were already present at the start of the simulation. 1024 electron particles per cell, and 16 + N 5 and + He 2 superparticles per cell were used for the simulation.
The ionisation module allows the specification of multiple electron species; the electrons from the ionisation of + N 5 and + N 6 were given their own species which allowed easier tracking and separation of the energy distributions for comparison with the preionised electrons. By separating the electrons in this way, figure 14 reveals that these deeply bound electrons are freed near enough to the trapping region that all are trapped and contribute to the accelerated electron population. The energy distribution of the helium and weakly bound nitrogen electrons compared to that of the strongly bound nitrogen electrons are found to be in reasonable quantitative agreement with McGuffey's results. The slightly lower maximum electron energy is attributed to the reduced self-focussing in 2D.

Collisional ionisation of carbon
In [59], 1D simulations from a collisional Fokker-Planck code are presented for a solid carbon target at initial densities between × − × 4 10 3 10 28 29 m −3 , with an incident 10 16 Wcm 2 , 350 fs, λ = 0.25 μm laser pulse by Town et al who make use of BED ionisation cross section modelling whilst neglecting recombination. This is reproduced using the collisional ionisation module presented in this paper. Carbon ions are chosen to be once ionised and their density is chosen to provide a neutralising background to the electrons. We use the initial electron density and temperature distribution provided in [59]; these are included as simple piecewise linear approximations for the electron temperature T e (in eV), the density ρ for both the electrons and + C ions (in − m 3 ) and the position x (in μm): ) an initially very high ionisation rate to + C 4 followed by ionisation to + C 6 over 1 ps which is in good agreement with Town et al; note that the much faster rise to the fourth ionisation state is simply a result of more frequent sampling of the plasma density [59]. The ionisation rate for laser-plasma interactions with carbon will be correctly predicted in cases where the electrons do not gain sufficient energy to cause ionisation before the laser field ionises the carbon up to + C 4 , since solutions for collisional ionisation with and without recombination converge as demonstrated in the results of Town et al [59].

Review of recent laser-plasma PIC simulations
Here we summarise some important recent applications of PIC codes in laser-plasma interactions. Within laser-plasma instability studies the most widely studied is stimulated Raman scattering (SRS). This is of central importance to both NIF-scale plasmas and shock ignition. Also being on the fast electron timescale requires shorter runtimes than stimulated Brillouin scattering (SBS). Nonetheless the simulations discussed below do include ion dynamics, e.g. through the Langmuir decay instability (LDI), which are essential ingredients of the full description of the effect of SRS on reflectivities and hot electron generation. Such simulations in 3D or 2D with hundreds of speckles represent the current state-of-art for LPI simulations. We do not cover recent work on the two plasmon decay instability (TPDI), which is crucial for conventional direct drive and important in shock ignition, as the aim here is to represent to complexity and realism acheivable by the combination of modern optimised PIC codes and HPC. This is well covered by a survey of just SRS where full 3D single speckle and 2D multi-speckle simulations are at the forefront of what is possible with PIC for laser-plasmas where collisionless, relativistic physics is all that is required. As examples of PIC code applications where additional physics is added, we discuss recent advances in laser-solid interaction using hybrid PIC schemes and QED processes for high intensity laser-plasma interaction. Preionised helium at neutral gas density of 10 19 cm −3 with a 1% nitrogen additive preionised into the + N 5 state; plots are shown 250 fs after the laser pulse enters the plasma. The top plot shows the energy distribution for the electrons produced from ionisation of + N 5 and + N 6 (upper curve) and that from preionised electrons (lower curve). The bottom plot shows the density of the + N 5 and + N 6 electrons plotted over the preionised electron density.

SRS
Detailed numerical studies of stimulated Raman scattering (SRS) have been ongoing for decades. Here we concetrate on the advancements in kinetic regime simulations afforded by the increased computing power available over the last decade. Single hot-spot experiments at the TRIDENT laser facility established that the transition from fluid to kinetic regimes was controlled by λ k De [60] where k is the wavenumber of the SRS daughter Langmuir wave, or electron plasma wave (EPW), and λ De is the electron Debye length.  [62] completed the first at scale 3D speckle simulation by assuming a ω 3 NIF laser, i.e. a wavelength of 351 nm, was focused into a diffraction limited beam of width 1.4 μm for an f /4 speckle and the simulation domain was × × 35 6 6 μm 3 . The initial plasma had = T 4 e keV, T i = 2 keV and a density = n n / 0.14 e cr which corresponded to λ = k 0.34 De for the initial EPW coupled to SRS. Simulations used 500 particles per cell on a × × 2048 256 256 grid and were only possible due to the use of the optimised VPIC PIC code [63]. These 3D single-speckle simulations were in qualitative agreement with previous 2D simulations and theory in that they reproduced the sequence of events expected in this regime. First the damping rate of EPWs due to electron trapping lowers the threshold for SRS compared to linear estimates. This electron trapping results in a lowering of the frequency of the EPW causing wave-front bowing of the daughter EPW. This bowing begins the process of SRS saturation by introducing a transverse (in the y-z simulation plane) phase variation in the EPW reducing the SRS drive term and allowing trapped electron transverse losses increasing the EPW damping. Eventually the EPW amplitude exceeds the threshold for the trapped particle modulational instability (TPMI) and it begins to self-focus, increasing further the transverse EPW phase variations and loss of trapped electrons which terminates the SRS pulse. Much of this can be see in figure 16 which is reproduced from figure 2 of Yin et al [62]. Importantly these simulations were able to test numerical convergence and compare to previous 2D work thereby allowing estimates of the accuracy of such 2D work which otherwise exhibited the same qualitative growth and saturation mechanisms. For example the onset threshold for SRS is increased from ≃ × 6.5 10 15 W cm −3 in 2D to ≃ × 1.5 10 16 W cm −3 in 3D, for f/4, and the saturation level dropped from 12% to 5.5%. These SRS results have also been confirmed by Rousseaux et al [64] and Mason-Laborde et al [65] who also considered the role of the Weibel instability on the trapped particle current. Collisions have also been shown in 1D simulations [66] to increase the threshold for SRS by detrapping electrons from weak EPWs and thereby increasing the Landau damping rate.
More recent PIC simulations [67] of SRS relevant for NIF using the OSIRIS PIC code [68] concentrated on intense laser speckles. For intensities of × 3 10 15 W cm −2 , temperatures of ∼3 keV and densities of ∼0.1 critical, a mechanism was highlighted by which electrons could be acclerated up to energies in excess of 100 keV. These energies are reached by the electrons interacting with a succession of EPWs of different phase velocity. This broadband EPW spectrum resulting from SRS rescatter of both forward and backward SRS and the subsequent LDI of the SRS EPWs.
The realistic scale 2D and 3D SRS simulations discussed above concentrated on initially uniform plasmas as appropriate for controlled single hot-spot experiments or perhaps the large-scale plasmas one might find in NIF hohlraums. For shock ignition (SI) both the compression lasers and the ignitor pulse are incident directly on the fuel pellet and necessarily tranverse significant density variations along a single speckle. 1D simulations in such configurations, e.g. [69,70], show that ∼70% of the laser incident energy can be converted to fast electron through SRS and cavitation leading to a hot-electrons with temperature less than 30 keV. This laser absorption takes place in the cavitation initiated by a resonator formed by SRS at the 1/4 and 1/16 critical density surfaces. Such a collisionless laser absorption process in the hot coronal plasma envisaged in SI may be advantageous to driving the final ignitor shock. Extending this work to 2D using the EMI2D PIC code presents an even more complex mixture of LPI [71]. Here SRS and SBS are restricted to the region of plasma below 1/4 critical density. For high intensity speckles the convectively unstable SRS leads to electron trapping, reducing Landau damping and an increased burst of SRS, as seen in [62], on sub-ps timescales. These inflationary SRS bursts are strongly space and time-dependent. LPI at the 1/4 critical surface generates rapidly varying density structures which act as a self-induced dynamic random phase plate inhibiting SBS beyong the 1/4 critial surface. For lower intensity speckles the inflationary Figure 15. 1D simulation of collisional ionisation of + C based on assumed initial electron temperature and density profiles over a μ 1.5 m domain, as given by [59]. The domain is split into 256 grid points with carbon ions at T = 15 eV and initial density matching the electrons. The simulation was carried out with 2000 particles per species per cell for μ < ( ) < x 0.25 m 0.7 over 1 ps. The average ionisation state is taken from 30 cells around the outermost cell.
SRS is absent and SRS is dominated by absolutely unstable growth at the 1/4 critical surface.
Whether for NIF scale hohlraums or SI direct drive all descriptions of SRS for laser systems using RPP smoothing must deal with the effects of multi-speckle interaction. Yin et al [72,73] have investigated the interaction of two speckles in 3D and hundreds of speckles in 2D in the kinetic regime relevant to NIF. Here it was found that the most intense speckles evolve in isolation, i.e. as in [62], but the side and longitudinal losses of trapped electrons reduce the damping rate of EPW in neighbouring speckles thereby triggering SRS in speckles which would otherwise be below the threshold for the onset of SRS. Overall the threshold for SRS is reduced, the saturated amplitudes scale only weakly with laser intensity but the fast electron flux does scale with laser intensity. Across the parameter regime studied the net SRS induced reflectivities were found to scale with λ ( ) − k De 4 suggesting that either slightly higher background electron temperatures or applied magnetic fields, to supress trapped particle speckle migration, may help to mitigate SRS reflectivity in some cases.

Short-pulse, laser-solid interaction
With the advent of chirped pulse amplification (CPA) [74] and the subsequent development and proliferation of high-intensity short-pulse laser systems [75], many of them petawatt class, the interaction of extreme fields with solid matter has become a lively research topic at the forefront of both experimental and computational plasma physics [76].
Short-pulse laser-solid interactions have found applications in: the study of materials in extreme conditions [77][78][79]; relativistic particle accelerators [80][81][82][83] and neutron sources [84]; ultrafast imaging [85], including x-ray sources [86,87]; high harmonic generation [88][89][90][91]; as well as alternative routes to Inertial Confinement Fusion, such as Fast Ignition [92][93][94][95]. In the presence of extremely high intensities, matter becomes rapidly ionised, placing free electrons directly into the laser field. The laser intensity is often expressed as the normalised laser amplitude ω = a eE m c / 0 0 e 0 , at intensities where a 0 > 1 (i.e. beyond = × I 1.37 10 0 18 W cm −2 for a 1 μm wavelength laser) the motion of these electrons in the laser field will be relativistic-quivering with energies in excess of their rest mass. All of the applications for high intensity short-pulse lasers mentioned above rely on the laser coupling effectively to the electrons near the target surface and driving a flux of energetic particles into the target. The ions in the plasma respond very slowly compared to the laser pulse duration. Ion motion is generally in response to large electrostatic fields generated via charge separation either at the front of the target, as electrons travel away from the laser interaction region, or at the rear surface, as they exit the target. Furthermore, many of the key processes are collisionless; the energetic electrons they produce travel a significant distance away from the laser interaction region, and heating occurs mainly via Ohmic heating due to a return current. This is in contrast to long pulse laser interactions where electrons excited in the laser field primarily transfer energy to the plasma in the vicinity of the interaction via collisional processes such as inverse bremsstrahlung, and heating over greater distances occurs predominantly via thermal conduction [47]. Although non-local effects are important in this context, they can be dealt with as a perturbation to the background population.
There are a number mechanisms which may be responsible for the absorption of high-intensity lasers and the generation of fast electrons: • Resonance absorption: A p-polarised electromagnetic wave obliquely incident on a plasma with a density gradient can propagate up to a density of θ = n n cos e c 2 due to refraction, where θ is the angle of incidence of the laser measure relative to the target normal, and a linear density ramp has been assumed. Beyond this density the electric field decays evanescently, but is still non-zero at the critical surface. At this point the plasma frequency, by definition, matches the frequency of the incident laser, and thus a plasma wave is resonantly driven [96,97]. The resonance acts to drive the plasma wave to very large amplitudes until it can no longer be supported and wave breaking occurs, resulting in the ejection of energetic electrons from the resonance region. The process of resonance absorption requires moderate (10 14 -10 17 W cm −2 ) laser intensities in order for the evanescent field to be strong enough to drive the electrons at the critical surface, but avoid deformation of the plasma density profile due to the ponderomotive force (discussed later). This process also requires the plasma conditions to be relatively uniform over the absorption region. • Vacuum heating: If the scale length of the plasma density ramp is very short compared to the laser wavelength, vacuum heating occurs [98]. Electrons in the overdense region (densities such that the electron plasma frequency is greater than the laser frequency) are directly exposed to the laser's electric field. Electrons near the surface of the plasma can be pulled, by the electric field of the laser, out of the plasma, beyond the Debye sheath. As the field reverses polarity the electrons are accelerated back into the plasma. By the time the field reverses again these electrons will have travelled beyond the skin depth ) of the target, and thus are no longer influenced by the laser fields. • Ponderomotive acceleration: The finite size and duration of the laser pulse tends to result in a net force on the plasma, driving charged particles away from regions of high laser intensity. An electron initially at rest in the middle of a Gaussian laser pulse will move under the influence of the laser's electric field. When the field reverses, however, the electron will have moved to a region of lower intensity (and thus field amplitude), and so experiences a weaker restoring force. Over time this leads to a net drift away from the peak electric field amplitude. An expression for the time-averaged ponderomotive force can be derived [99]: where ω is the frequency, and E 0 the electric field amplitude, of the laser. This effect, parallel and perpendicular to the laser pulse, can lead to channelling (in underdense plasma (see, for example figure 19) or hole boring (in overdense plasma) [100]. Note that the ponderomotive force is independent of the sign of the charge, and so will also act to expel ions from regions of high laser intensity. • × J B force: at high laser intensities ( ≳ I 10 18 W cm −2 ) the force due to the magnetic field of the laser acting on electrons oscillating in the electric field is non-negligible compared to the electric field force. For a lone electron this results in a lemniscate orbit rather than simple oscillation. If the laser is incident upon an overdense plasma, however, a similar effect to vacuum heating occurs, except with a longitudinal driving force, due to the magnetic field, at twice the laser frequency [101].
Although recent theoretical work [102] has been successful in bounding laser absorption, and the various processes are individually well understood, the dominant effects, the interplay between them and, most importantly, the characteristics of the fast particle populations which carry the energy away from the absorption region, are not dictated by theory or suitably constrained by experiment, and are sensitive to a broad range of laser and plasma conditions. It is often difficult to diagnose high intensity laser matter interactions, which are characterised by small spatial scales and fleeting timescales, despite the dramatic effect they ultimately drive: shocks; isochoric heating; x-ray emission (continuum and line); fast particles; and scattered electromagnetic radiation (see figure 17).
Fine scale physical processes can have a material impact on the macroscopic system. For example, the strong currents of fast particles leads to the generation of 10 T 4 fieldsstrong enough to oppose laser field and re-direct the fast electron flux. The ω 2 0 bunching of electrons, characteristic of × J B absorption, has been shown to drive Langmuir turbulence [103] which contributes to plasma heating. Progress in understanding short-pulse laser solid interaction necessitates detailed modelling [104], alongside experimental and theoretical advances. However, the level of complexity outlined here presents a significant challenge for computational physics.
Laser-solid interaction typically takes place over a relatively small volume (on the order of a few 100 μm 3 ), but the density contrast in this region gives rise to extremely short scale lengths. Resolving these density features requires a fine spatial mesh, this is often a more severe constraint than the need for accurate dispersion of EM waves, for example. In a PIC code, the presence of such steep profiles makes the problem of numerical self-heating particularly acute. Left unchecked, self-heating can lead to the expansion of the target surface prior to the arrival of the laser; effectively altering the plasma conditions encountered from those requested. Macroparticle weights are also a potential issue for laser-solid interaction. Figure 18 shows the fraction of laser energy absorbed into fast electrons (defined as those with an energy greater than ten times the average electron energy) for a number of equivalent 1, 2 and 3D simulations of the interaction of a 100 J, 0.53 μm pulse with a pre-ionised carbon target at a density of 64 n c preceded by an exponential density ramp with a scale-length of 1 μm. There is a clear correlation between the absorption fraction and the mean fast electron weight, which appears to converge as the particle weight falls (i.e. as the total number of particles in the problem is increased). May et al have recently reported anomalous macro-particle stopping [105] (also discussed in [92]). In the macro-particle description adopted in PIC codes the charge-to-mass ratio, q/m, is the same for all particles of a given species, irrespective of the their weight. This allows the simulation to run with a range of particle weights (often the only practical way to resolve density features) since the particle push itself is dependent only on q/m and not the particle weight. However, some processes are sensitive to the macro-particle weight. For example, slowing by plasmon emission (i.e. the wakefield generation of electron plasma waves) is dependent on q 2 /m. Consequentially, the macro-particles with difference weights will behave differently. This means that the simulation as a whole will display a sensitivity to particle numbers which is not simply due to numerical self-heating.
It is clear then that the choice of mesh and particle parameters, separate from the physical system, can be critical to the efficacy of PIC modelling of problems in short-pulse laser-solid  The observed absorption appears to be dependent on the macro-particle weight and the dimensionality of the system. The absorption converges as the average particle weight falls, suggesting that poor particle counts can contribute to an underestimate of the absorbed energy. Two data points from 'typical' 2D simulations, such as those used in [116] are also given. In this case it was not feasible to run the 3D simulations at an equivalent level of fidelity.
interaction. Fine spatial meshes, high particle counts, and high order particle shapes are a necessity. Modelling using direct Vlasov solvers [106][107][108] avoids some of these problems, the fidelity of the continuum description employed in an Eulerian Vlasov solver is not sensitive to the local particle density. However, the computational expense of such an approach currently rules out 3D modelling (as this would necessitate advancing the solution on a 6D phase space mesh), and 2D modelling is often prohibitively expensive, particularly at the highest intensities due to the need to service a momentum mesh extending up to high βγ while still resolving the cold background Maxwellian. 1D systems are tractable and offer a mechanism to study laser absorption under certain conditions (see [106] and [104]), but this rules out a number of the absorption processes discussed above, and conservation of canonical momentum ensures that the electron divergence, a key factor in the onward transport of fast electrons, cannot meaningfully be derived from such an approach. The noisefree properties of direct Vlasov solvers, and the promise of collision operators [109,110] freed from the restrictions of the macro-particle approach, make them a desirable tool in select cases, but the associated computational cost can make them impractical for routine use in the design and evaluation of short-pulse laser-matter experiments. The robust, flexible and (by comparison) computationally tractable PIC approach is the only practical option for pursuing a modelling capability to support a broad range of short-pulse experiments.

Fast-electron transport
In practice many of the applications of short-pulse laser-solid interaction span a range of spatial and temporal scales. While the key physics of laser absorption and fast particle generation might be localised close to the relativistically corrected critical density surface (i.e. γ = n n e c ), the location, relative to the solid material, and morphology of this surface as well as the local plasma conditions (e.g. ionisation, scale length and temperature), and those along the beam path, are generally determined by effects on longer time-scales and over greater volumes [111,112].
An example of the PIC code EPOCH forming part of an integrated modelling suite to simulate a material properies experiment of the type discussed in [78] is shown in figure 19. An aluminium microdot buried in a diamond foil is compressed with a third harmonic long pulse and heated with 500 J of first harmonic (i.e. λ = 1 0 μm) light. The plot of electron densities shows that the main short pulse has been able to channel through the low density preplasma and deposits it's energy close to the solid target. The hot electron flux generated is used as a source term for a hybrid MC-VFP model of electron transport at solid densities (following the approach in [116]). This model is subcycled within the radiation hydrodynamics model to capture effects such as radiation losses, expansion of the heated sample and, over longer timescales the dissasembly of the target. A density in excess of 10 g cc −1 is achieved but despite delivering ∼150 J in hot electrons into the target, the temperature only reaches ∼250 eV-this is due to the high energies of the fast electrons generated, due in part to the interaction of the laser with the low density pre-plasma before it reaches the solid material. This serves to highlight the complexity of short-pulse laser-solid interaction under realistic experimental conditions.
The 'pedestal' that typically precedes a CPA pulse can itself be of considerable intensity, generating significant preplasma (see figure 19) ahead of the solid material. Progress has been made to reduce the pre-pulse in the current generation of glass lasers, leading to extremely high contrast beams [113,114], and the conversion to higher harmonics has been able to deliver high contrast for some time, albeit at the cost of total energy [115]. However pre-pulse effects remain an issue in most deployed high power short-pulse systems-particularly those at the highest energies. In some situations, such as fast ignition (in the absence of a re-entrant cone) [92,95,112] and material properties experiments [78], the target conditions which the short-pulse laser encounters are be determined by a low-intensity, long-pulse interaction over many nanoseconds.
In these cases, a complete description of the problem naturally requires an extensive range of physics from radiation hydrodynamics effects for the long pulse interaction and target response to short timescale kinetic effects where the laser is absorbed and collisional effects in the transport of the fast particles generated. One approach to this problem has been to integrate a number of codes which capture the detailed physics for each spatial and temporal scale [116,117]. Here the PIC code forms part of a simulation framework which, at present, offers the only tractable solution for the complex multiscale modelling problems posed by experiments on the current generation of high-power lasers, see figure 19.
However, this approach introduces its own potential sources of error. It introduces an artificial 'transport gap' between the laser interaction region and the fast electron probe plane in which collisional effects may begin to have an important effect, and beam-plasma instabilities may be seeded. Since each code is suited to a different range of densities and temperatures the minimal overlap in densities for which collisionless PIC codes and hybrid codes are valid may result in phenomena such as magnetic field formation, beam filamentation and seeding of instabilities, which might occur at the intermediate densities, being neglected. Furthermore, any short-pulse driven hydrodynamic motion (for example ponderomotive hole boring) which may have a significant effect on the region described by the hybrid model at the highest intensities [118] cannot be captured self-consistently without the complexity and computational expense of sub-cycling the PIC model within a host hydro-code together with the hybrid fast electron transport algorithm. The lack of feedback from the transport model to the PIC code prevents transport phenomena affecting the laser absorption via the return current.
An alternative approach is to broaden the capability of the PIC code to cover fast electron transport through dense material. This implies the inclusion of some collisional physics which, as we have seen, can require extremely high particle counts. Progress can be made by altering the algorithm at the highest densities to relax the requirement to resolve high-frequency phenomena such as Langmuir waves, which would be expected to be critically damped, and EM waves, which should not penetrate much beyond the relativistically corrected critical density surface [14,119,46].
The changes made by Cohen et al [119] to the PSC, and also incorporated into OSIRIS [46], allow for a more consistent approach to laser absorption and transport by employing a conventional, fully electromagnetic PIC model in the laser interaction region (i.e. electron densities ≲100 n c ) and a reduced field solver, with many similarities to the hybrid approach used in the transport models discussed above, in the high density region. A equivalent algorithm has been developed for the PIC code EPOCH and is outlined briefly here.
The need to resolve the Debye length and electron plasma frequency can be avoided by updating the electric fields via Ohm's law, instead of Ampère's law: where the b, i and f subscripts denote thermal background electron, ion and energetic electron properties, respectively. This does not permit electromagnetic or electron plasma waves to propagate. This approximation is valid at high densities, since such waves would normally be heavily damped by Debye screening and collisional effects, respectively. To distinguish fast electrons from the background, the electron population is partitioned into two species. Background electrons are promoted to the fast electron population if their velocity satisfies the condition: , where T b is the local background electron temperature and α is a free parameter (typically ∼5-10).
By retaining a kinetic model for all of the particle species throughout the simulation domain, rather than following the standard hybrid approach of a fluid model for the background plasma, PIC codes can easily employ Ohm's law in high density grid cells concurrently with a standard Maxwell field The use of an Ohmic field solver alongside a standard Maxwell solver provides a short pulse laser modelling capability which is less prone to non-physical numerical effects, such as self-heating; is capable of modelling high density plasma kinetics with more relaxed time constraints; can efficiently model collision-dominated regimes using currently available computing resources; provides a means of checking that no significant transport phenomena are neglected when using collisionless PIC simulation data as an energetic electron source in Monte-Carlo transport codes; can model the effect of energetic electron transport phenomena and target heating on the laser interaction physics; and is capable of simulating both the production of energetic electrons via laser interaction and their subsequent transport through the high density plasma within a single simulation code.

QED effects
The development of Monte-Carlo models such as that outlined in section 4 [37,38] and their coupling to PIC codes [27][28][29] occurred rapidly after the initial calculations showing that QED processes will be important in next generation laser plasma interactions [30,120]. The first coupled QED-PIC simulations demonstrated the coupling of QED and plasma processes in these interactions. Nerush et al performed simulations of a pair cascade induced by two counter-propagating laser-pulses and showed that the cascade is limited by plasma effects [27]. Similar work in the astrophysics community showed that the same effect can occur in pulsar magnetospheres [121]. Ridgers et al performed the first simulations of a 10 PW laser-solid interaction including QED effects and showed that: (i) the complex plasma processes occurring in a dense plasma strongly affect the rates of the QED processes [29]; (ii) Conversely, QED processes can strongly alter the energy spectra of the plasma constituents and so be expected to modify the plasma processes [122]. Figures 21 and 22 are reproduced from [122] and show an example of the strong gamma-ray and pair production seen in these simulations and its effect on the electron and ion energy spectra. Recently coupled QED-PIC simulations have been used to determine quantitative scaling laws for the conversion laser energy to gamma-ray photons [123][124][125]. Quantitative predictions of the effect on both radiation-pressure ion acceleration [126] and relativistic transparency [127] have also been made.
Although higher intensities than those available on current PW-class systems are required to enter the regime where QED processes affect the plasma dynamics, it is possible to observe the important QED processes in the relevant very non-linear regime at intensities achievable today. If the electrons are accelerated to much higher energies than they acquire from the laser pulse alone, for example by a conventional particle accelerator or laser-wakefield acceleration, the resulting increase in the electric field in the electron rest frame leads to a dramatic increase in η. Experiments of this type have measured the important processes in the weakly non-linear regime [128], but there has been recent interest in using state-of-the-art high intensity laser systems to explore the very non-linear regime. Coupled QED-PIC codes have been used to simulate this type of experiment [129] which would allow benchmarking of the current QED model used to make predictions about next-generation laser-matter interactions.

Conclusions
Over the decades, PIC codes have become ubiquitous tools in many areas of plasma physics. They offer flexibility and stability, but their greatest strength is their underlying simplicity. While direct Vlasov methods [50][51][52]108] offer a much higher fidelity solution, they are often too costly, in terms of compute, for routine use in many problems. The ratio of the computational cost of direct Vlasov methods relative to PIC for collisionless problems was estimated by Besse et al [130] where N v is the number of points in momentum space in the Vlasov code, d v is the number of degrees of freedom in v and n pic is the number of particles per cell in the PIC description. Noise levels in a PIC code are inversely proportional to n pic . With one degree of freedom, direct Vlasov is clearly the best approach; for two the choice would be problem dependent. With three degrees of freedom in velocity, and two or three spatial dimensions, PIC is often the only practical approach for all but the smallest systems and shortest temporal scales. The results presented above show that this picture changes dramatically once collisions are included, and that for PIC codes to accurately reproduce Braginskii Figure 20. Lineouts of background electron temperature along y = 0, at t = 100 fs, from a simple 2D diffusion test problem. The system consisted of a uniform density background plasma with an initial temperature of 100 eV and an additional population of energetic electrons (at 10% the background density and 100 keV) initially located in the x < 0 region. Simulations using collisional EPOCH with a Maxwellian field solver throughout the domain (labeled EPOCH-C, dot-dashed blue) and with the Ohmic field solver in the x > 0 region (labeled EPOCH-H, solid red) were compared with the results of a simulation using the Monte-Carlo electron transport code THOR [116] (dashed green). The electron source injected in the THOR simulation was sampled from the energetic electrons recorded crossing the x = 0 plane in the EPOCH simulations. 2D simulations of prolific gamma-ray and pair production in next-generation laser plasma interactions, taken from [122]. (a) Gamma-ray photon (2D blue) and positron (red contours) production in the interaction of a 12.5 PW laser pulse with a solid aluminium target (ion density-3D grey). (b) Equivalent plot for a 320 PW laser pulse striking solid aluminium. Reprinted with permission from [122]. Copyright 2012, AIP Publishing LLC. Figure 22. Energy spectra of the plasma components with and without the inclusion of QED effects in both the 12.5 PW and 320 PW lasersolid simulations in figure 21, taken from [122]. (a) and (c) show electron spectra for 12.5 PW and 320 PW respectively, (b) and (d) ion spectra for 12.5 PW and 320 PW. Reprinted with permission from [122]. Copyright 2012, AIP Publishing LLC. and non-local transport, thousands of particles per cell may be needed. In such cases VFP techniques, e.g. [50,51], are likely to be quicker, and noise free, for most laser-plasma problems. Collisional routines in PIC codes will still be useful as a means of including some thermalisation processes and collisional ionisation routines but if detailed modelling is sensitive to recovering transport then currently the best option for PIC is to adopt a hybrid approach, as in [119], or use VFP. The hybrid approach [119] does still require an accurate collision operator for the particles although these collisional processes are no longer required to reproduce transport correctly by themselves.
Recent algorithmic developments and the availability of large scale computing now mean that the domain of applicability for PIC methods in plasma physics is much extended. Realistic sub-grid physics models mean that many properties of a fluid-like background medium can be included in the simulation, for example, detailed elastic (Coulomb) and inelastic (ionising) collisions, without the requirement to resolve the Debye length and while maintaining fully consistent treatment of high and low energy particles. Analytically intractable problems such as the transport properties of non-Maxwellian magnetised plasmas will be modelled correctly by the current PIC methods.
The basic explicit PIC algorithm has been refined with the addition of higher order weights and interpolation schemes, and augmented with additional physics, broadening the range of physical conditions which can be considered. Here we have detailed the approaches taken in the development of the PIC code EPOCH, and used it to demonstrate the robustness of modern PIC codes and the resilience of algorithms based on high order particle weights and charge-preserving particle pushes to numerical self heating. Furthermore, we have demonstrated the efficacy of binary collision and ionisation models, including a novel partial superparticle ionisation scheme. The algorithms adopted here, and the tests used to underwrite them, while broadly applicable, are geared toward problems in laser-plasma interaction where a daunting range of physical conditions can be explored in a single laboratory experiment. The challenge this posses for numerical modelling is driving the continued development of PIC codes, including; novel hybrid algorithms [119], the inclusion of QED effects [131], and the coupling of kinetic and hydrodynamic models [116]. With efforts to exploit next generation computer architectures already underway, e.g. [132], it is clear that PIC codes will remain a vital tool for many years to come; but their suitability, stability, and validity should not be taken for granted and must be routinely tested as the basic PIC code continues to develop at the heart of an evolving, complex, and sophisticated, numerical toolkit.