Intense laser-generated ion beams in plasmas: the rapid heating regime

Previous work (Robinson 2022 Plasma Phys. Control. Fusion 64 105014) has shown that there is a regime of ion beam transport in plasmas where the electron–ion drag becomes so intense that the electron temperature rapidly rises to temperatures at which the electron–ion drag is strongly reduced. This regime is mainly thought to pertain to short-pulse laser-generated ion beams, because of the ion energy fluences required to enter this regime. The regime for was thought to be determined by a dimensionless number, Ξh . Here it is shown that there is another dimensionless number that needs to be considered. Numerical simulations show that this dimensionless number, χ h , does not affect the half-energy range of the ion beam, but does determine how strong stopping is at the ion beam front, and thus how rapidly the beam energy falls past the half-energy range.

Recently this author pointed out that interpreting these experiments has been fraught with difficulty [2,4,21]. There are many factors that complicate the interpretation of these experiments, and, in a previous paper of Robinson [21] it was shown that, once the intensity of a laser-generated proton beam Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. becomes sufficiently large, the drag heating induced by proton beam can influence the proton range via the sharp rise in electron temperature.
In this paper we seek to expand on the physics of this strong heating regime. Previously [21] we characterised the strong heating regime by a dimensionless number, Ξ h (see equation (3) in section 4). Here we show that this is not the only dimensionless number that needs to be considered, and that one may also need to consider χ h = I i /m e n i v 3 i (see equation (10) in section 4). Through a series of numerical simulations we show that Ξ h > 1 is the primary condition for entering the strong heating regime, however χ h is still important because once χ h >1, the ions are not fully stopped at the beam front, and thus the beam can still propagate at reduced energy beyond the half-energy range. This could have a number of implications for ion beam transport, including both target heating and nuclear reaction yield. Thus we conclude that assessing the importance of strong heating requires considering both the Ξ h number put forward in [21], and the χ h number put forward in this paper.

Notation
The main problem considered in this paper is that of a high energy ion or plasma beam that is incident on a static target plasma. This means that direct reference to the target ions or 'background' ions is seldom required. We have therefore adopted the following convention: when we refer to electron quantities, e.g. n e , this should be taken as a reference to the target or background (i.e. n e is the electron density of the target plasma). Likewise, when we refer to ion quantities, this should be taken as a reference to the incident ion or plasma beam. So terms such as v i ,Z i , and n i refer to the ion velocity, ion charge state, and ion density of the incident ion beam and not that of the target/background plasma. We use this convention throughout this paper, unless stated otherwise. Note that we also work in units of eV for temperature.

Background
The drag force acting on a suprathermal ion in a classical plasma is often given in a form similar that below. Note that Z i and v i correspond the ion charge and ion velocity respectively, and n e is the electron density. The term, L, is a dimensionless 'stopping number', where G is the function that corrects the stopping power due to the finite temperature of the electrons. Physically this effect emerges because the electron temperature affects both binary collisions and plasma wave excitation. In particular, once the electron thermal velocity approaches the ion beam velocity, one expects the stopping power to be significantly curtailed. The correction that is often applied, the Chandrasekhar correction, has the form, where y = v i /v th,e . In the cold limit G → 1, however when y → 0, G → 0 which means that the drag will be greatly reduced in sufficiently hot plasmas. Even if a plasma is initially cool, the electron drag itself must heat up the plasma. If the ion beam energy flux density becomes very large then this raises the question of whether a 'self-heating' regime exists where the drag will be automatically curtailed by naturally raising the electron temperature. The quantitive determination of this is the subject of this paper, and a preceding paper [21].

Theory
In a previous paper [21], the authors characterised the onset of the regime where the electron drag on the incident ion beam results in such rapid heating that the electron-ion drag becomes suppressed for much of the ion beam. This can increase the range of the ion beam in a cold target by around a factor of 10. The onset of the regime was characterised via the dimensionless number, with the strong-heating regime occurring when Ξ h > 1. Let us now both re-derive this and obtain another dimensionless measure, by considering three important length-scales. The length-scales are: (1) the length associated with the ion beam, v i τ i , (2) the length behind the front of the beam at which the critical temperature is reached, λ h , and (3) the stopping distance of the ions in the cold plasma, λ s . The critical temperature is the temperature at which the electron thermal velocity equals the incident ion velocity, Accounting for electron-ion drag alone, the heating equation becomes, and from these two equations, we can immediately obtain an estimate for λ h , In general the electron-ion drag on each ion will be, and if we take the cold limit of this (G = 1) and we substitute into equation (6) then we obtain, One condition for entering the strong-heating regime is that λ h < v i τ i , i.e. the distance required to reach the critical temperature is shorter than the length of the ion pulse. This condition is essentially the same as that shown previously, as it is quickly found that, However there is another stipulation that we might put on ratios of these length-scales. We may be interested in the case where λ h < λ s . In this case the critical temperature is reached over a distance shorter than stopping range of the ions in the cold plasma. This implies that ions at the leading front of the ion beam will not be strongly stopped as the cold plasma is heated up, and it implies that there may be a regime where even greater enhancements in beam-range are possible. This is there regime where χ h > 1, where, which follows from both equation (8), and from taking the cold stopping range to be, which in turn follows from equation (7). Note that we can immediately write the ratio between the two dimensionless numbers as, so unless the ion pulse duration becomes quite short, we will expect Ξ h > χ h over a wide range of parameter space. In order to enter the strong heating regime, where the electron-ion drag is suppressed, we will still require Ξ h > 1. The outstanding question is how important the regime where Ξ h > 1 and χ h > 1? To answer this question we have carried out a numerical simulation study which is reported in the following sections.

Simulation method
To address the questions raised in section 4, we have employed the 1D hybrid code apsis. This is similar to the version used in earlier work by Robinson [21], however certain modifications will be highlighted. This treats the high energy ion beam kinetically, the target ions as a fluid, and the electrons are treated as a neutralizing massless fluid. This has certain similarities to a code previously developed by Sherlock [22], certain types of 'hybrid particle-in-cell (PIC)' codes [23,24], and is related to a simulation paradigm commonly used for 'space plasmas' [25]. The kinetic treatment of the ion beam is resolved by a macroparticle numerical model, as is the case in a traditional PIC code. The kinetic ion macroparticles are thus evolved via the equations of motion: and, where the 'e' and 'i' subscripts in the drag term indicate that both drag due to the fluid electrons and the fluid ions is included. In addition the kinetic ions are subject to angular scattering from the fluid ions. The fluid ions evolve via hydrodynamic equations, i.e.
and, ∂v ∂t as well as a fluid ion energy equation. The electrons are treated as a fluid that moves to maintain quasi-neutrality, i.e.
where n fl i and n ki i are the ion fluid density and kinetic ion density respectively, along with an electron energy equation to evolve the electron temperature. The fluid electron and fluid ion energy equations are coupled via electron-ion equilibration terms. The electric field is determined by the Ohm's Law, where F drag is the net drag per unit volume on the ions. Open boundary conditions are used for the kinetic ions, and for the fluid ions we use open boundary conditions at one end (where kinetic ions are injected) and reflective boundary conditions at the other.

Simulation setup
Using apsis, we have carried out a sets of simulations to examine the issues raised in previous sections. Unless stated otherwise, we have used a grid with n x = 30 000 cells, and a cell size of ∆x = 0.1 µm, giving a standard domain size of 3000 µm. The time step was set to ∆t = 2 fs. A total of 1.2 × 10 6 macroparticles were used for the fast ions, with 120 macroparticles being injected per time step. Throughout we have considered a deuteron beam incident on a pure deuterium plasma at a uniform ion density of n i = 4.85 × 10 29 m −3 , with Z ⋆ = 1, and an initial electron and ion temperature of 100 eV. The deuteron beam was monoenergetic with a velocity of v i = 1.55 × 10 7 ms −1 (2.5 MeV)in all cases. Simulations were run up to 180 ps or longer as required. The simulations were divided into three sets (I, II, and III) with simulations in each set denoted by A, B, C, etc. We thus use terms of the form X-Y to label simulations with X being I, II, or III and Y being A, B, C, etc. Details of the simulations are given in tables 1, 2, and 3 for sets I, II, and III respectively.
The purposes of these different sets can be explained as follows: Set I is the case where the ion energy flux is naturally varied with constant ion pulse duration. Since both dimensionless numbers (Ξ h and χ h ) are proportional to the ion energy flux, both increase linearly across this set. In Set II both the ion energy flux and the ion pulse duration are varied in such a way as to keep Ξ h constant while χ h varies linearly. Finally in Set III we keep the ion energy flux constant and vary the ion pulse duration so that χ h is fixed, but Ξ h varies.

Simulation results and discussion
In the first instance we have analysed the results of these simulations in terms of R 1/2 , which was previously used in [21]. The metric R 1/2 is defined as the distance reached by the last macroparticle to fall to half its initial energy at the time that it falls to half its initial energy. The results for R 1/2 from Set I are plotted in figure 1.
The single particle stopping length of a 2.5 MeV deuteron in the ambient plasma can be obtained from equation (11) which yields 180 µm. One thus finds from figure 1 that the results from parameters used in this paper in concordance with  those obtained previously [21], that is to say that R 1/2 increases with ion flux energy and quickly becomes greater than the single particle stopping length (reaching values up to around ten times the stopping length). As in [21] this can be attributed to the strong drag heating that raises the electron temperature to the point where the drag is strongly reduced. From this  'intensity scan' it is not clear whether the role of Ξ h or χ h is more important. In Set II, Ξ h has a value that is fixed by varying both the ion pulse duration and the intensity, and χ h is thus varied freely. Plots of R 1/2 for Set II are shown in figure 2.
From figure 2 it can be seen that the variation of R 1/2 with χ h and fixed Ξ h is weak, which indicates that R 1/2 is largely non-dependent on χ h and is primarily determined by Ξ h . The results for R 1/2 from Set III are plotted in figure 3.
From figure 3 it can be seen that the inference drawn from Set II is confirmed: when χ h is fixed, but Ξ h is varied, there is a strong variation in R 1/2 comparable to that seen in Set I. We can therefore conclude that Ξ h primarily determines R 1/2 . Thus the conclusions reached in [21] remain unchanged.
Although χ h does not affect the gross properties of initial ion penetration, it does have consequences for the overall evolution of the ion beam. This can seen by comparing the ion phase space for simulations II-A and II-E. The v x − x phase  space of II-E at 100 ps is shown in figure 4, and that of II-A in figure 5. Figure 4 shows that, in simulation II-E (Ξ h = 4.9, χ h = 1.52) the ion beam is being progressively slowed down and that ions are eventually only retaining less than 5% of their initial energy once they can no longer keep up with the beam front. Note that the initial deuteron velocity in all these simulations is 15.5 × 10 6 ms −1 . It is clear that something radically different is happening in simulation II-A. The ions have not been 'completely stopped' at the beam front, and have thus retained more energy than the ions in II-E which reach the beam front. In figure 6 the ion v x − x phase space in II-A is shown at a series of times. In figure 6 the consequences of partial stopping at the beam front can be seen more clearly. The partial stopping and varying drag across the spatial extent of the beam leads to a much more complex evolution in phase space, with some portions overtaking others leading to what might be described as 'folding' or 'tumbling' of the distribution function.
It is thus reasonable to say that although the values of R 1/2 are similar for calculations with the same Ξ h , the simulations where χ h ≫ 1 also exhibit a slow decline of the beam energy beyond R 1/2 rather than the abrupt stopping observed when χ h is low. This implies that the rapid heating regime (high Ξ h and χ h ) at least has the potential for larger nuclear reaction yields, although this depends on the reaction cross-sections over a wide range of energies.

Discussion of limitations and implications
It is important to discuss the limitations and implications of this study. The first aspect to consider is the ion energy range. If ion drag is the only heating term considered (as it is here), then going to much higher ion energies than those we have considered would reduce the drag heating so much as to put one out of the 'self-heating' regime, and this is reflected in the dimensionless numbers we have derived. On the other hand, the electrons that are co-moving with the ion beam will scatter from the static (target) ions and this should produce an Ohmic heating effect. This is unlikely to be significant in the parameter range we have restricted ourselves to in this study, but could be significant at higher ion energies. It would therefore appear that there are open questions about what precisely happens at higher ion energies than those considered here.
One potential limitation of the analytic theory that should be addressed is that electron-ion equilibration is not included in the derivation of both Ξ h and χ h . The thermal energy transfer between electrons and ions is, however, fully included in our numerical calculations. There are good reasons for neglecting e-i equilibration in the analytic estimates. In the limit where equilibration is arbitrarily fast, this would effectively shift the specific heat capacity (in units of eV for temperature) from 3en e /2 to 3en i (Z + 1)/2, which means that the estimates for Ξ h and χ h can only be shifted by a factor of 2 at most (Z = 1 case), however the shift rapidly becomes progressively smaller once Z > 1. That noted, equilibration takes a finite time which is not very smaller than the ion beam duration in the typical scenarios that we have considered. This makes the error introduced by neglecting e-i equilibration even smaller.
An implication of this study that should be considered is that any focussing or defocussing processes in ion beam transport can have major implications for the overall beam range provided that the beam parameters are close to the strong heating regime. If the beam experiences focussing then by increasing I i this will cause the Ξ h and χ h parameters to increase which will lead to increased beam range. On the other hand, if the beam experiences defocussing then the reduction in I i will cause Ξ h and χ h to fall. One can imagine that this will cause the termination of the strong heating regime in a number of possible scenarios.

Conclusions
In summary, in this paper we have shown that there is another dimensionless number which needs to be considered in the context of intense ion beams propagating in plasmas, namely χ h (equation (10)), in addition to the Ξ h number (equation (3)) introduced in prior work. Here we have used a number of numerical calculations to show that χ h does not seem to have an impact on the gross aspects of the initial ion penetration, i.e. it does not affect the half-energy range (R 1/2 ). It is clear, however, that χ h affects the degree of stopping at the ion beam front. At low χ h stopping at the ion beam front is almost 'total' i.e. ions reaching the beam front effectively lose all their energy. At high χ h , stopping at the beam front is much weaker. Thus the ion beam can propagate well beyond R 1/2 , albeit at progressively reduced energy. This could have implications for calculations of nuclear reaction yields in scenarios involving intense ion beams and would therefore have to be considered.

Data availability statement
The data cannot be made publicly available upon publication because no suitable repository exists for hosting data in this field of study. The data that support the findings of this study are available upon reasonable request from the authors.