Hybrid Vlasov–Fokker–Planck–Maxwell simulations of fast electron transport and the time dependance of K-shell excitation in a mid-Z metallic target

Using a hybrid Vlasov–Fokker–Planck–Maxwell code coupled to calculations using a modified relativistic binary-encounter Bethe model of ionizing collisions we study the time and space dependence of Kα photon generation by a fast electron beam injected into a solid density copper plasma target. The electron beam is chosen to be representative of that expected to be generated by a picosecond duration, ∼1019 W cm−2 intensity laser. Kα photons are produced as electrons reflux laterally across the target and are slowed and thermalized by collisions with, and Ohmic heating of, the background fluid over a ∼10 ps timescale, which dictates the timescale for Kα emission. The results show reasonable agreement with recent experimental results in terms of both the yield and time dependance. We show how lateral expansion of the electrons can be imaged in the Kα radiation.


Introduction
Fast electron transport has been a vigorously studied area in high intensity laser-solid interactions in recent times. It is motivated by high-energy-density physics applications including, but not limited to: advanced ignition schemes [1] for inertial confinement fusion research; K α x-ray generation for use in probing dense material [2][3][4][5], laser driven shocks [6] and even biological matter [7][8][9]. Understanding how fast electrons propagate in solid targets is of fundamental interest, but is also crucial to optimizing these applications. However, the dynamics of fast electrons in laser-solid interactions is extremely complicated.
For lasers exceeding 10 18 W cm −2 [10], laser energy is absorbed into a non-thermal beam of MeV electrons that enter the solid density target and propagate to the rear, setting up resistive electric fields, a return current and magnetic fields within the bulk of the target [11,12]. On reaching the rear, the lack of a dielectric response results in the formation of extremely strong sheath fields [13] that reflect the majority of electrons back into the target, whereupon they 'reflux' back and forth, expanding laterally. The electrons eventually lose energy by a combination of collisional energy loss-predominantly with background ions-and Ohmic heating of the solid target. Collisions with bound inner shell electrons will result in the emission of spectral lines. As well as being useful for the probing applications mentioned above, the resulting K α radiation produced has been used to diagnose the propagation of electrons within a solid target [5,[14][15][16][17][18][19][20][21][22][23][24]. The resulting diagnostic can provide much needed temporally and spatially resolved data on the transport of fast electrons in solid targets, allowing us to resolve and understand their complex dynamics.
There have been numerous experimental measurements of K α emission in such interactions as a diagnostic of the electron transport for fast ignition applications. These have generated overall yields [5,[14][15][16], spatial distributions [17][18][19][20][21][22], and even temporal information for the K α emission [23,24]. In order to understand this experimental data, numerical simulations are required. There have been a number of attempts to couple kinetic codes with crosssections for K α emission, generally using hybrid particle-in-cell approaches [25][26][27]. Here we take a different approach, using a hybrid Vlasov-Fokker-Planck code [28]. The advantage of using a Vlasov distribution as opposed to discrete particles is that it is straightforward to integrate the distribution function and obtain rate equations, allowing both spatial and temporal K α emission information in addition to including collisions and electromagnetic effects simultaneously. We will then go on to show that such a code can reproduce recent experimental K α measurements [24] and moreover allows us to understand how the radiated spectrum reveals information about the injected distribution of fast-electrons and the relative importance of collisional and Ohmic heating on the transport of the fast-electrons.
The structure of this paper will be as follows. First the electron transport model will be described, as used in [28,29]. Then we will introduce a method for calculating ionizing collisions from the distribution function to generate a rate equation for K α generation. We then describe results for a copper target interacting with a picosecond duration, ∼10 19 W cm −2 intensity laser, which are conditions similar to recent experiments [24]. After concluding, an appendix summarizes fits to data used in the code for elements with atomic numbers spanning 10 Z 80.

Fast electron transport and Fido
To model the propagation of energetic electrons in solid density material we used a hybrid Vlasov-Fokker-Planck code, Fido [28,29]. This code is designed specifically for electron transport problems in solid density plasma targets. The hybrid mode uses a background solid density plasma, modeled as separate electron and ion fluids using classical transport equations, and a fast electron population modeled by solving the Vlasov-Fokker-Planck equation, with appropriate simplifications, in a two-dimensional (2D) Cartesian slab geometry (x, y). Maxwell's equations are solved with a contribution to the electric field from the fluid background through Ohm's law in a standard hybrid method [30]. The system is chosen to have self-consistent electric, E x , E y , and magnetic, B z , field components only. Fido uses a spherical polar coordinate system ( p, θ, φ) for the momentum space for the high energy part of the electron population, allowing accurate, conservative magnetic rotations and electron-ion pitchangle scattering. The resulting relativistic Vlasov-Fokker-Planck equation for fast electron distribution function f eh (x, y, p, θ, φ, t) injected into background electron and ion fluids, with number densities n ec and n i respectively, is given by the following expression in this coordinate system [31]: where I (x, y, p, θ, φ, t) is a function describing the injection of fast electrons into the background, v = p/γ m e , γ = 1 + p 2 /m 2 e c 2 and Y ec, i = 4π(Z ec, i e 2 /4π 0 ) 2 ln ei . This equation describes the evolution of the electron distribution function in the presence of statistically smoothed electromagnetic fields (left hand terms of equation (1)) and momentum-space operators describing the effects of slowing (dynamic friction) and scattering (dynamic diffusion) of particles due to two-body Coulombic interactions (right hand terms of equation (1) respectively) between the fast particles and background (interactions between fast electrons are ignored).
In the left hand side advective terms, fluxes are calculated with piecewise-parabolicinterpolation [32], with the overall algorithm being second-order accurate. The electron slowing down term (first on right hand side of equation (1)) is calculated using a monotonicity, particle number and positivity preserving upwind scheme. The electron scattering terms (second and third terms on right hand side of equation (1)) are calculated using an implicit scheme with second-order center-differencing in momentum space.
The basic computational cycle on a time step t consists of the injection of new electrons through the source term I , followed by calculating moments from the distribution function to obtain macroscopic quantities such as j. Current density j is then used to calculate an electric field after a half time step using a standard hybrid scheme [30], directly from the fast particle current, j f , via E = −ηj f by assuming the fast electron current is exactly balanced by a cold return current arising as a consequence of the resistive electric field in the background fluid. η is the standard scalar resistivity of the fluid background according to Braginskii [33]. The magnetic field is calculated implicitly using a backward differenced version of Faraday's equation and the advanced electric field. In vacuum regions, a standard finite difference time domain algorithm is used to calculate fields using Maxwell's equations and the fast electron current.
The fast electron distribution function is then advected in position space for a half time step, followed by advection in momentum space for a full step along with Ohmic heating of the background fluid. Electron-ion and electron-electron collisions are calculated and the change in energy of the fast electrons is transferred to the background fluid. A second half time step position space advection then occurs, followed by a final field update. Hydrodynamic advection of the background fluid may follow, but here the background density was fixed. The time step size for the each integration step varies, and is calculated to be the minimum of either of the Courant-Friedrichs-Lewy [34] conditions for both position and momentum space advection.

The excitation of K -shell electrons
Excitation of K -shell electrons can occur by a number of different mechanisms, including photoionization and impact ionization by energetic electrons. For mid-Z materials, we expect that collisional excitation by electrons is dominant. The rate at which single hole K -shell electronic excitations occur due to electron collisions for an electron distribution function f e (x, p, θ, φ, t) in material of ion density n i and of atomic number Z is: where R K 1 (x, t) is the rate of single hole K -shell electronic excitations per unit volume, and σ K (Z , p) is the Z and p dependent cross-section for electronic excitation of the first K -shell electron. n i = n i − n K 1 is the ion density minus the number density of atoms with single holes in the K -shell, n K 1 . The factor of 2 is extracted from the excitation cross-section for reasons of convenience, as the cross-section for excitation from a full shell is twice that of a shell with a vacancy. Hence, the rate at which double hole K -shell electronic excitations occur is where n K 1 = n K 1 − n K 2 . The number density of both singly and doubly excited K -shell ions can be calculated from these rate equations, if it is of interest to distinguish the two cases. The generation of new vacancies in the K -shell is balanced by the decay of those excited states by radiative, Auger and Coster-Kronig processes [35,36]. The total rate of change of K -shell vacancies per unit volume n K , due to contributions from both singly and doubly excited ions, is therefore given by where K is the K -shell natural level width, and is related to the lifetime of the excited state τ K by τ K =h/ K [35].

K α production
The rate of K α photon emission density, R K α , is related to the decay rate K of K -shell hole density n K by the relationship where ω K is the K -shell fluorescence yield (the proportion of K -shell decays that are radiative in nature, as opposed to Auger electrons for example) and f K α is the fractional emission rate of all K α lines (the proportion of radiative K -shell transitions that are from the L-shell). f K α is related to the ratio of K α to K β emission intensities by the relationship [37] The cross-section for K -shell ionization that we use here is similar to a theoretical crosssection due to Santos et al [38], which was used in [16]. This is a relativistic impact crosssection based on the Bethe cross-section [39][40][41][42], consistent with the Møller cross-section [39,40,43], and is parameterized in terms of the K -shell binding energy B and the kinetic energy of the incident electron T . This is useful for fast electron transport modeling in particular as it is both relativistic (there is a significant difference between relativistic and non-relativistic cross-sections for energies of interest), and is also a continuous mathematical function. This modified relativistic binary-encounter Bethe (MRBEB) [44] cross-section is expressed as where x = X/m e c 2 and β 2 x = 1 − 1/(x + 1) 2 , with X and x representing either T or B and t or b. m e is the electron mass and c is the speed of light. t b = t/b, a 0 = 5.292 × 10 −10 m is the Bohr radius, h = 6.626 × 10 −34 J s is Planck's constant and α = 1/137.0 is the fine structure constant. This cross-section is used in the place of σ K (Z , p) in the rate equations above, with B = B(Z ) and T = T ( p). C 1s1/2 (Z ) is a scaling factor related to the shielding of the nuclear charge [44].

Numerical calculation of K-alpha cross-sections
The rate equation, (5), requires integration over the momentum space distribution of the electrons. The finite-differenced Vlasov distribution from the Fido code is amenable to this calculation, as it is simply a matter of summing over the array. The finite differenced momentum space form of the rate equation, (5), is where K = 4π and is the isotropic part of the finite differenced electron distribution function { f e } k,l,m , where k, l, m are integers denoting points on the gridded spherical distribution, with l=N θ −1 l=0 sin θ l dθ l = 2 and m=N φ −1 m=0 dφ m = 2π , and the momentum space is defined from 0 to p max . p k , {σ K } k , and v k are the gridded momentum space, σ RBEB , and velocity space quantities respectively. p k is the momentum grid spacing.
In equation (9), n K appears on the right-hand side, and so this equation needs to be time integrated, in conjunction with the time-stepping integration of the fast electron transport simulation with step size t. Because it contains a decay equation it will tend to be unstable in a simple forward or center time differencing scheme as the time-step we wish to take may be large compared with the K -shell lifetime (the stability condition is in general α t < 2, where α is a decay constant). However, it is reasonable to assume that f e does not change considerably over t, as otherwise the electron transport dynamical timescale is not resolved. Hence, over time step t, K is approximately constant, and therefore an exact solution for the decay equation can be found: where γ = K /h + K and n n K = n K (n t) ≡ n K (t). Equation (6) does not have a stability issue and can simply be time integrated using simple Euler integration, but it is trivial to solve the equation analytically. The rate of K α emission, averaged over time step t is

Energy loss by ionizing collisions
The effect of energy lost by ionizing collisions with K -shell electrons on the distribution function can be written as The rate of energy loss due to ionizing collisions is the rate of collisions n i σ K v multiplied by the energy loss per collision, which is the K -shell binding energy B. Hence This can easily be added to the collision terms on the right hand side of equation (1) and has been implemented in the code. Although energy loss through ionizing collisions with L-shell, M-shell, etc electrons should also be taken into account for consistency, the binding energies are significantly smaller than the K -shell, and to include them would greatly increase the complexity of the calculation. Hence, for some degree of self-consistency from the point of view of the collision operator for the Fokker-Planck calculation, the plasma is assumed to be ionized only up to the K -shell, with the K -shell electrons shielding the nucleus so that effectively the plasma has a electron number density of n e = (Z − 2)n i . A characteristic radius of the K -shell electron wavefunction is r K = a 0 /Z . In the Fokker-Planck collision formalism, the cross-section is an integral over limits running from the 90 • scattering impact parameter b π/2 = Z e 2 /12π 0 k B T e [45] (orh/ p, whichever is the larger) to the Debye length, as otherwise the integral diverges. The ratio of r K to b π/2 is or r K /b π/2 k B T e (eV)/9Z 2 . Experimental studies of K α production are generally focused on transition metals, particularly copper. Taking this as a reference element, r K /b π/2 ≈ k B T e (eV)/8 × 10 3 , which means that in situations of interest r K b π/2 and hence it is reasonable to consider collisions with a nucleus shielded by the K -shell electrons. To take into account additional ionized K -shell electrons, the background fluid Z is updated using the calculated change in K -shell vacancies, Z = n K /n i . In principle, this model could be extended to include L-shell, M-shell, etc electrons also, although the collision modeling is then rather complex.

Results
Here our model is applied to a system similar to recent experimental measurements [24] of interactions with 500 2 × 20 µm 3 copper targets. The simulation domain consists of a uniform density copper target with dimensions L y = 20 µm by L x = 500 µm in an (x, y) 2D slab geometry. The initial temperature of the bulk target is chosen to be k B T e = 400 eV, which corresponds to an ionization level of Z free = 18.7 according to a Thomas-Fermi model [46], that is to say that the average atom is Neon like. Hence, the electron-ion Coulomb logarithm would take a value of ln ei = 0.67 [47], which renders the Fokker-Planck operator invalid, since it involves an expansion in powers of 1/ln ei . We therefore introduced a cut-off such that ln ei = 2 if the calculated value would be less than 2. It should be noted, however, that this value for ln ei is only for the background plasma electrons, which are modeled as a fluid. Fast electrons modeled by the Vlasov distribution effectively have ln ei > 10. In addition, the background plasma rapidly increases in temperature through a combination of Ohmic and collisional heating. Incidentally, the maximum temperature achieved in the simulation (2 keV) corresponds to Z free = 27, so the K -shell electrons are always expected to be bound.
The boundary conditions are reflecting at y = 0 and L y , and in the x direction the last few grid points before the boundary have an exponentially increasing step size, so that the boundaries are effectively very far from the interaction. The reflecting boundary conditions on the y boundaries are used instead of calculating the sheath fields, as this drastically slows down the calculation because of the need to obey the Courant-Friedrichs-Lewy [34] condition in momentum space, t < p/E. Particles are specularly reflected back into the domain, as they would be expected to behave in the rapidly formed sheath fields provided the spatio-temporal variation in the fields is small compared with the excursion of the electron [29].
An electron beam is injected into a thin region centered at x = 250 µm and y = 0 with a spatial profile that is uniform in the y direction with a width of I y0 = 0.1 µm and has a Gaussian profile in the x direction with a width of I x0 = 5 µm, i.e.
Electrons are injected with a temporal profile that is a sinusoidal function for a half cycle with a full width at half maximum of 1 ps. The electron beam momentum distribution is defined as in [28]. The beam momentum, p, distribution is a shifted Gaussian with center momentum p 0 = ( 1 + I /1.3 × 10 18 W cm −2 (λ 0 µm −1 ) 2 − 1)2m e c 2 /3 [10]. I is chosen to be representative of a laser with an intensity of 3 × 10 19 W cm −2 , with the injected electron beam number density being calculated using the effective laser energy with an assumed 0.3 absorption fraction. The angular distribution in θ is a Gaussian with a 1/e angle of π/4 rad, although in practice this matters little due to the strong collimating magnetic fields that are formed.   Figure 2 shows integrated time series data for (a) the total system energy, (b) energy in fast electrons, (c) energy lost to collisions and (d) energy lost through Ohmic heating. The total system energy does not include the initial thermal energy of the background fluid, and is equal to the assumed laser energy multiplied by the absorption fraction, i.e. the laser energy is 12 J, so the total energy injected into the system is 3.6 J. Energy is calculated in the 2D system by using the injected beam diameter as the characteristic length for the third dimension. Although there is a large contribution to the energy loss from Ohmic heating of the background fluid, slowing We clearly expect that since the hot electron population is confined to the target, and has an average energy much higher than the inner shell ionization energy for copper, that this timescale also dictates the timescale for K α emission. Figure 3 shows data from a numerical 'streak camera', which shows the temporal evolution of the K α emission in the simulation. Panels (a) and (c) show the K α photon emission rate integrated in the y only, and in (b) and (d) the K α emission has been integrated over the whole domain. In all cases the horizontal axis represents time in the simulation. Panels (a) and (b) show the raw data, with a very asymmetric shape in time that falls back to almost zero by 10 ps from the start of the simulation.
The data in these temporal evolution plots are equivalent to streak camera that could be obtained on a real experiment. Panels (c) and (d) show the same data as (a) and (b), but numerically convolved with a Gaussian instrument function of width 2 ps. This is to approximate the effects of the response of a real streak camera, as in [24]. The full width at half maximum duration of the resulting convolved peak shown in figure 3(d) is 4.1 ps, which is close to the reported width of 5.5 ps in that reference for very similar conditions. Figure 4 shows the total emitted K α energy as a function of time in the simulation. The yield of K α photons was calculated by using the injected beam diameter as the characteristic length for the third dimension. This approaches 0.025 J of total energy, which corresponds to a conversion efficiency from laser energy to total K α photon energy of 2 × 10 −3 , consistent with previous calculations [14]. The efficiency is higher than may be realistic in an experiment due to the high initial temperature of the plasma, which is chosen because of numerical constraints. When the fluid background is initialized with a 2000 eV temperature, there is an order of magnitude lower amount of Ohmic heating compared with the 400 eV case, as would be expected from the difference in resistivity. This also reduces the full width at half maximum time of K α emission to 3.9 ps, and has a different shape, with the emission dropping off more rapidly after the peak.  In (c, d), the data has been convolved with a Gaussian instrument function with a 2 ps full width at half maximum.
The reason for this reduction in the time is likely to be due to the larger fraction of energy loss by collisions as opposed to Ohmic heating. In the lower temperature run, Ohmic heating represents energy loss from the fast electron population, but some energy is also transferred into the strong electric fields, reaching 10 10 V m −1 . In addition to decelerating the fastest electrons, it will accelerate slower electrons and will therefore increase the length of time over which inner shell collisional excitations occur, since these only require ∼10 keV electrons. It should be noted that this is likely not the only reason for the experimentally observed lengthened timescale for emission [24]. As can be seen from figure 1(c), the temperature of the background plasma is very high in the vicinity of the initial injection point of the electron current. This means that the ionization level of the copper in this region is likely to be approaching Z free = 27, which means the radiation emission would be helium-like lines instead of K α . This will result in a reduction of the K α emission in the focal region and therefore an apparent lengthening of the timescale for K α emission will result. However, it is also clear that Ohmic heating will have an even larger effect for a colder target, and therefore cannot be ignored.
From the virtual streak camera (and therefore possibly a real streak camera image), we can also obtain a physical picture of the electrons lateral transport in the target. In figure 5, the same time series data for K α emission as in figure 3 is shown, but overlaid are lines indicating different expansion rates relative to the speed of light c. The fastest electrons expand laterally at 60% of the speed of light, but the bulk of the fast electrons expand more slowly.

Conclusions
Using a hybrid Vlasov-Fokker-Planck-Maxwell code coupled to calculations using a modified relativistic binary-encounter Bethe model for ionizing collisions we studied the time and space dependence of K α photon generation by a fast electron beam injected into a solid density copper plasma target. K α photons are produced as electrons reflux laterally across the target and are thermalized by collisions with the background over a ∼10 ps timescale, which dictates the timescale for K α emission. The results show reasonable agreement with recent experimental results in terms of both the yield and time dependance. In addition, it is clear that Ohmic heating plays a role in the electron transport and subsequently the K α production.
It is interesting to note that agreement with experiment has been achieved despite the lack of any fitting to experimental measurements in these calculations. That is to say; the model is built from (existing) theory of laser energy absorption [28], electron transport in solids [30] and collisional ionization [44], with the only empirically based parameters being the absorption fraction of the laser and the initial cone angle of the injected electrons (which appears to be unimportant). The model adapted here by no means encompasses all the physics going on in the target, but it is clear that it is sufficient to obtain a basic understanding of the mechanisms. This indicates that diagnostics such as lateral expansion of the electrons imaged through the K α radiation can yield a direct insight into the dynamics in the target, provided a sufficiently accurate model is adopted.

Acknowledgments
This work was supported by the Defense Threat Reduction Agency, Basic Research Award # HDTRA1-10-1-0077, to the University of Michigan. AGRT acknowledges Jonathan Davies and Philip Nilson for useful discussions.

Appendix. Useful fits to the tabulated data
To implement these calculations in the Vlasov-Fokker-Planck code, it is useful to create fits to the data in table A.1, for reasons of convenience and also because in the numerical calculation, Z at a point in space may be a volume averaged atomic number and therefore may be fractional. Fits were chosen for a combination of accuracy and also for trends that would not widely diverge outside of the range of the data table, as a pure polynomial fit may do. The fits have been characterized by their relative percentage error f (Z ), which is defined for a parameter X (Z ) with a fit X f (Z ) as f (Z ) = 100|X (Z ) − X f (Z )|/ X f (Z ).  [38,48], ω K is the fluorescence yield, from [36], K is the K -shell natural line width in eV, from [35], K β /K α is the ratio of emission rates of K β to K α lines from [49], f K α is the fractional emission rate of all K α lines. (For Z = 10 to Z = 80.) For ω K , a fit exists due to Wentzel [50] as used in [51], but we find that a better fit to the fluorescence yields tabulated in [36] is ω K = Z 3.8 Z 3.8 + 4.568 × 10 5 . (A.1) The relative percentage error in this fit compared to the tabulated data from [36] is shown in figure A.1. The relative percentage error is small except for small values of Z . This is because the absolute error is actually very small but the value of ω K tends to zero for small Z , so that the relative error increases.
The natural lifetime of a K -shell hole, K , is obtained from tables in [35]. The fit we use to this is an exponential relationship:  The relative percentage error in this fit compared to the tabulated data from [35] is shown in figure A.2. As for the fit to ω K , the relative error for small Z increases dramatically. Again, the absolute error does not increase very much but the size of K itself becomes very small, leading to the large relative error.
The ratio of emission rates of K β to K α lines, K β /K α , allows us to calculate the number of radiative transfers which are K α decays as opposed to K β . The fit to K β /K α is more difficult as there are oscillations in the data. Rather than using a complicated fit, however, a simple logarithmic relation is used, for Z 10, which smoothly continues to higher Z . The resulting  The relative percentage error in this fit compared to the tabulated data from [49] is shown in figure A.3. Using a high order polynomial fit could be more accurate here, but at higher Z a polynomial fit will start to diverge rapidly, and the resulting accuracy with the smooth logarithmic fit is everywhere less than 2.3%. The K -shell binding energy B used in the relativistic cross-section calculation is tabulated in [38,48], from calculations using Dirac-Fock wavefunctions. The data points from Z = 10 to 37 are taken from [38] and data points from Z = 38 to 103 are taken from [48] (only data up to The relative percentage error in this fit compared to the tabulated data from [38] is shown in figure A.4.