Laser-driven ion beams are appealing for many applications in a broad range of research fields ranging from nuclear physics1–3 to cancer therapy.4 Targets can be strongly heated on picosecond time scales by using laser-driven ion beams, enabling fundamental studies of high-energy-density matter,5 ion fast ignition (IFI),6,7 neutron sources,8,9 and other topics.10

Most of the ion beam–plasma interaction models used so far have relied on strong assumptions, such as perfectly focused beams interacting with plasma only by Coulomb collisions.11–13 However, this optimistic view vanishes when new physical effects beyond the standard assumptions are accounted for such as the effects on ion beam energy deposition of the stopping power modeling, the generation of self-generated resistive fields, and the initial beam divergence. The differences between the stopping power models found in the literature,14 especially near the Bragg peak, have been assessed recently in a dedicated proton stopping experiment. The stopping powers measured are in good agreement with the nonstandard stopping theory based on the T-matrix15 and the Brown–Preston–Singleton (BPS) model,16 while there are substantial differences with the standard models widely used so far.17 As the ranges predicted by the T-matrix theory and the BPS model are higher than those predicted by the standard models,18 the optimal ion kinetic energies used in many applications may change. Therefore, experimentally validated stopping power models should be used in ion beam–plasma interaction studies to get realistic predictions.

In some applications such as IFI, the ion beam intensities are high enough to produce strong self-generated resistive fields, which may give rise to “anomalous” energy deposition and beam compression. With the development of hybrid simulation methods,19–22 such “collective effects” have recently been reported in the context of creation of high-energy-density matter. It is found that the large magnetic field growth in collective interactions can result in self-focusing of the beam and much stronger localized heating of the target.23–26 Here, we go a step further and perform a systematic study of the influence of self-generated fields up to the beam current densities typical of IFI by means of two-dimensional (2D) and three-dimensional (3D) hybrid simulations of ion beam propagation in solid-density resistive media. Ion divergence is an important additional effect on beam propagation. In a typical target normal sheath acceleration (TNSA) scheme, ion beams usually emerge with a relatively high divergence angle,27,28 which will have a critical impact on the ignition energies of IFI.29

In this paper, the aforementioned three effects, namely, stopping power model, resistive fields, and ion divergence, will be presented and discussed, highlighting their importance for the generation of high-energy-density matter. The paper deals exclusively with intense proton beams, and for simplicity does not address the complexity associated with an evolving ion charge of the beam. The remainder of the paper is organized as follows. In Sec. II, the model used in our simulations, including the stopping power, is briefly described. In Sec. III, the three effects are studied by means of hybrid calculations. Finally, in Sec. IV, conclusions are presented and the outlook for future work is summarized.

Simulations of ion energy deposition have been conducted by means of a 2D/3D hybrid code similar to that used for electron-driven fast ignition.19,20 The hybrid code solves the Fokker–Planck equation for collisional ion energy deposition29 and includes particle-in-cell (PIC) routines for ion motion in self-generated electromagnetic fields, which are obtained from a simplified version of the Maxwell equations. Cylindrical r–z coordinates are used for 2D simulations and Cartesian coordinates for 3D simulations. In this section, the stopping power and the plasma resistivity used in our hybrid simulations are described.

Ion stopping power calculation in warm and hot dense matter can be decomposed into electron and ion contributions. In the warm dense matter (WDM) regime, both bound and free electrons contribute to the total ion stopping. In fully ionized high-temperature plasmas, the ion contribution due to Coulomb collisions between fast and plasma ions can also be important. Thus, in general, the total stopping power can be written asdEdxtotal=dEdxbound+dEdxfree+dEdxion.(1)

Cayzac et al.14 have pointed out discrepancies between different models for stopping power calculations, specially near the Bragg peak. An experiment to discriminate between those models was carried out, concluding that only nonstandard theories such as that based on the T-matrix15 and the BPS model16 agreed with the experiments.17 Thus, we have used in our calculations the experimentally validated BPS stopping model for free electron and bare ion contributions. The BPS theory relies on the dimensional continuation method to compute the energy loss for a nonrelativistic particle moving through a fully ionized plasma and includes an exact treatment of two-body quantum scattering. To show the differences between the BPS model and classical stopping power models, such as those described in Ref. 29, we compare the stopping power of 5 MeV monoenergetic protons in fully ionized aluminum at solid density and 1 keV in Fig. 1. Note the different penetration depths and Bragg peaks for each model. The stopping power obtained by the classical theory is about 10% higher than that from the BPS model, and therefore the corresponding range is about 10% shorter for the classical case. From now on, we will refer to the model described in Ref. 29 as the “classical” stopping power model.

For bound electron stopping, the Barriga-Carrasco (BC) method30 has been used, owing to its simplicity and accuracy. In this method, the stopping power S is based on the Bethe formula:31S=Zeffeωpvp2Lb(vp),(2)where Z_eff is the effective charge state with ionization, ω_p is the plasma frequency, v_p is the velocity of the projectile, and L_b is the corresponding classical Coulomb logarithm. L_b can be approximated by an interpolation between the high- and low-velocity limits, avoiding negative values at low velocities:30,32Lb(vp)=LH(vp)=ln2vp2I−2Kvp2,vp>vint,LB(vp)=αvp31+Gvp2,vp≤vint,(3)where vint=3K+1.5I, K and I are the kinetic electron energy and the mean ionization potential, respectively, given by Barriga-Carrasco and Casas30 and Casas et al.,33α=1.067K/I2 is the friction coefficient for low velocities,34 and G is determined from the continuity condition L_H(v_int) = L_B(v_int).

The penetration depth of monoenergetic protons as a function of plasma temperature is shown in Fig. 2 for different stopping power models. Here, we compare the penetration ranges obtained by (i) the theoretical model SCAALP based on the local density approximation as described in Refs. 35 and 36, (ii) the results reported by Kim et al.,23 and (iii) the BPS and classical stopping models for free electrons and the BC model for bound electrons (henceforth referred to as the BPS and BC model). Note that there are no relevant differences in proton penetration at low temperatures for all proton energies. This is because the bound electron stopping is dominant at low temperatures, where ionization is not yet important. As expected, substantial differences appear at higher temperatures (>10 eV), especially for low-energy protons (1 MeV). The results of Kim et al.23 overestimate the penetration depth by more than 50%, and the classical stopping power model underestimates it by about 15% at temperatures higher than 100 eV. On the contrary, the results from the BPS and BC model fit the theoretical results from SCAALP quite well for the energies and temperatures shown in Fig. 2. Therefore, we will adopt the BPS and BC stopping power model throughout this work because of its accuracy.

The ranges of 5 MeV protons in several materials at solid density obtained by the BPS and BC model are compared in Fig. 3 for a wide range of temperatures. The proton range varies according to the material density. For heavier material, the proton feels a greater stopping power and the range is shorter. Moreover, range lengthening starts to be relevant for temperatures higher than 0.3 keV, and the range is substantially higher for the BPS model than for the classical model with the aluminum target temperature around 1 keV owing to the significant difference in electron stopping power between the two models. However, with the target temperature continuing to rise, the ion contribution becomes dominant. Especially for the aluminum target, the difference in ion contribution of the total stopping power between the two models is not obvious, and so the range difference decreases with increasing temperature.

As pointed out by Kim et al.,23,24 self-generated resistive fields can play a role in intense ion beam propagation in resistive plasmas, especially at very high current densities. In this case, ions generate an azimuthal B-field at the beam edge, which may compress it and even change its energy deposition pattern.

As we intend to study the ion beam energy deposition for a wide range of beam parameters, from generation of WDM to fast ignition of inertial fusion targets, we have performed our simulations by means of the hybrid code used for heating of solid matter by laser-generated electrons37 and electron-driven fast ignition.19 In this model, resistive fields are calculated from Ohm’s law combined with the simplified Maxwell’s equations without displacement current. Specifically, the equations solved by the hybrid code for field generation can be written as follows:38E=ηjr,∇×B=μ0jb+jr,∇×E=−∂B∂t,(4)where η is the plasma resistivity, and j_b and j_r are the ion beam and plasma return current densities, respectively. Note that field generation depends only on the plasma resistivity, which is computed by means of the Eidmann–Chimier model39,40 including phonon excitation for conducting materials. This model is based on an interpolation between the Spitzer regime for hot plasmas and the electron–phonon interaction in the limit of temperatures below the Fermi temperature. Resistivities for solid aluminum, copper, and gold obtained by this model are shown in Fig. 4 as functions of plasma temperature. The resistivities increase with temperature T up to a few tens of eV and then enter the Spitzer regime, where plasma resistivity scales as T^−3/2. In this regime, heavier materials show higher resistivities owing to their higher free electron densities.

We have considered a model problem similar to that described by Kim et al.23,24 for studying the proton beam energy deposition in solid aluminum targets at different current densities. Here, in addition to the aluminum targets, we have considered solid copper and gold as target materials in order to explore the importance of resistive fields for higher atomic numbers. In addition, we have extended the current density range by spanning in our simulations from 10⁹ to 10¹³ A/cm², two orders of magnitude higher than the maximum current density 10¹¹ A/cm² reported by Kim et al. This allows us to extend our study from the generation of WDM (10⁹–10¹⁰ A/cm²) up to fast ignition of inertial fusion targets (10¹³ A/cm²).41 The interaction of such high-intensity beams with solid, high-resistivity targets is directly applicable to fast ignition, where intense ion beams should pass through the high-Z material of the cone tip.

We assume in this paper that a monoenergetic ion beam impinges on a solid target. We have considered protons, not heavier ions, in order to avoid the effective charge uncertainties arising for such ions.42–46 The dynamics of proton energy deposition in solid aluminum targets has been explained in detail by Kim et al.23,24 We assume that protons are injected at the left surface of the simulation box. The comoving electrons are not considered owing to their very slight effect on proton transport along the target.23 The initial distribution of protons is Gaussian in radius, n(r)=n0exp[−(log2)(r/Δr1/2)2], with half-width at half-maximum (HWHM) Δr_1/2 = 14 μm and Gaussian in angle with HWHM Δθ_1/2 as a parameter varying from 0° to 20°, where θ is the polar angle. The ion pulse has a flattop shape with a duration of 3 ps and constant current densities from 10⁹ to 10¹³ A/cm², which correspond to beam energies from 0.13 J to 1.3 kJ, respectively. Electron and ion temperatures and plasma ionization are computed each time step from the SESAME tables.47 We assume that the target is initially at room temperature. The mesh width used in the simulations is 0.5 μm for each coordinate, and the time step is kept constant and set to 1.5 fs.

The hybrid simulations presented in this section are similar to those performed by Kim et al.,23,24 who pointed out the importance of collective effects on ion beam transport and energy deposition, including beam compression and filamentation structures. Here, we extend these results to different target materials and a wider range of beam current densities, and we include initial beam divergence.

The energy deposition of a perfectly focused 5 MeV proton beam in solid aluminum with current density as a parameter is presented here. Resistive field generation has been artificially suppressed. We aim at analyzing the effects of the stopping power modeling on the proton beam energy deposition. The beam parameters are described above.

The penetration depth as a function of current density is shown in Fig. 5(a) for both the BPS and classical stopping models. The differences are small at low and high current densities. At low current densities, the stopping by bound electrons is the most important contribution, because plasma temperatures and ionization are low. At high current densities, the plasma temperatures are high enough (>7 keV at 10¹³ A/cm²) for full ionization of the aluminum, and the stopping by ions dominates. For a solid aluminum target, there is not much difference in ion stopping at high temperature. For intermediate current densities, the beam penetration from the BPS model is about 10% higher than that from the classical model. As shown in Figs. 5(b) and 5(c), these range differences indeed affect the energy deposition in the target.

Resistive fields may play a role in the transport of intense ion beams in a resistive medium such as aluminum. For instance, a beam with a current density of 10¹¹ A/cm² in a target with a resistivity of 10⁻⁶ Ω m generates an electric field of 10⁹ V/m. If the beam radius is 10 μm, an azimuthal B-field will grow up to values of roughly 300 T in 3 ps. Such a B-field is high enough to bend the ion trajectories and change the energy deposition pattern. Therefore, we have conducted hybrid simulations of proton transport and energy deposition in solid targets to study this effect and assess the importance of resistive fields for ion propagation.

Figure 6 shows the azimuthal B-field distribution generated by proton beams propagating in a solid aluminum target for different current densities. The B-field grows at the beam edge owing to the radial variations in current density and plasma resistivity.38 As expected, the azimuthal B-field depends strongly on the current density, reaching peak values from about 300 T for 10¹¹ A/cm² to 1 kT for 10¹³ A/cm². As a consequence of these high B-fields, the protons are subjected to a radial Lorentz force that compresses the beam, as evidenced in Fig. 6(a). Note that for higher current densities, the shape of the B-field changes dramatically in the second half of the range. This is due to the scattering of the beam ions by the fully ionized plasma ions. When the current density is relatively low, 10¹⁰–10¹¹ A/cm², the target is only partially ionized and ion–ion scattering is almost negligible. For higher current densities, 10¹² A/cm², the plasma temperature rises to 3 keV, the aluminum is almost fully ionized, ion–ion scattering starts to be important, and the B-field generation is reduced, as shown in Fig. 6(b). When the current density is even higher, 10¹³ A/cm², the temperature rises to more than 7 keV and ion–ion scattering is significant, producing the beam expansion shown in Fig. 6(c).

The electron temperature distributions at 20 ps for beam and target conditions similar to those shown in Fig. 6 are plotted in Fig. 7. The right and left panels depict the temperatures with and without resistive fields, respectively. When the current density is relatively small, 10¹⁰ A/cm², resistive fields do not play a significant role: there is only a slight temperature increase near the end of the beam penetration [Figs. 7(a) and 7(b)]. For higher current densities, 10¹¹ A/cm², the beam is compressed by the B-field [Figs. 7(c) and 7(d)], and the target heating is enhanced. Owing to this additional target heating, range lengthening occurs, leading to greater penetration by the protons. When the current density reaches 10¹² A/cm², the beam is compressed, the proton range lengthens, and ion–ion scattering starts to be relevant near the end of the range. For a even higher current density of 10¹³ A/cm², close to that required for fast ignition, the target is fully ionized, ion scattering is relevant over most of the energy deposition zone, and the proton beam is not compressed by the azimuthal B-field. It is worth remarking that in this case, the Ohmic heating by the return current increases the energy deposition along the proton path, and, as a result, the range of protons is slightly lower than the range obtained in the collisional simulations of Fig. 7(g).

It is worth pointing out that the range lengthening due to beam compression shown in Figs. 7(c) and 7(d) is substantially lower than that reported in Ref. 23. In our simulations, 5 MeV protons penetrate 200 μm, whereas in that reference the penetration is about 300 μm. This can be explained partly by the differences in the stopping power model. For instance, the model used in Ref. 23 yields substantially higher ranges than ours at high temperatures (>200 eV) and low proton energies (1 MeV), as shown in Fig. 2.

PIC simulations and experiments have shown that even for the simpler laser-driven acceleration scheme, ions emerge with a non-negligible divergence angle.27 Hence, it is interesting to study the effects of beam divergence on ion transport and energy deposition. It is particularly important to determine the conditions under which beam divergence can balance beam compression due to self-generated azimuthal B-fields.

Figure 8 shows the effect of proton divergence in terms of temperature maps of the aluminum target without (left panels) and with (right panels) resistive fields. For perfectly collimated beams, the beam is compressed by the azimuthal B-field [Fig. 8(b)]. On the contrary, for a significant initial divergence, ion beams start spreading [Figs. 8(c), 8(e), and 8(g)], with their penetration range being reduced owing to the lower target heating. In addition, beam divergence reduces the beam current density within the target, and thus the resistive field generation and the beam compression. Beam compression partially balances the proton divergence for a low initial divergence half-angle, 5° in Fig. 8(d), for which self-generated B-fields are still strong enough to compress the beam. For a large divergence, θ = 15°, the beam compression almost disappears [Fig. 8(h)].

We have checked the effect of the initial beam divergence on the energy balance. Figure 9 shows the fraction of beam energy that is deposited in the target as Ohmic heating for different divergence angles and current densities. Note that Ohmic heating is non-negligible even for large beam divergences. As expected, for a fixed divergence angle, the Ohmic heating increases with increasing current density. For a fixed current density, the Ohmic heating decreases with increasing divergence angles, owing to the beam spread. It is worth noting that the energy losses due to Ohmic heating are relevant for large current densities and high divergence angles and should be taken into account in the energy balance. For instance, the Ohmic heating fraction of a 10¹³ A/cm² beam is higher than 10% for an initial divergence half-angle of 20°.

To understand the effects of resistivity on self-generated fields, we have extended our study to ion transport in heavier materials. These materials have a higher resistivity and higher collisionality, which are manifested in the two opposite trends of a rise and a reduction in field generation, respectively. In this subsection, we analyze which one of these two trends prevails.

The rise in the azimuthal B-field and the reduction in the ion penetration for heavier materials are shown in Fig. 10. The B-field and electron temperature at 20 ps are depicted in the left and right panels, respectively. It is important to note the filamentary structures that appear in the temperature profiles between 50 and 100 μm in the copper and gold targets. These are a clear signature of resistive filamentation of the ion beam, which is more evident in heavier materials. Note that this signature does not appear in Figs. 10(a), 10(c), and 10(e), because the azimuthal B-field plots were obtained 17 ps after the end of the ion pulse (3 ps). At such a long time, the B-fields have been diffused throughout the background plasma sweeping the B-field filamentary structures. It should be pointed out that such structures were also reported in the hybrid simulations described in Ref. 23. As the filamentary structures cannot be properly resolved in 2D settings, a 3D analysis is performed at the end of this subsection.

The ion–ion scattering mentioned in Sec. III B is only seen for the aluminum target [Fig. 10(b)], which is almost fully ionized for a beam current density of 10¹² A/cm² and an aluminum plasma temperature of 2.8 keV. On the contrary, ion–ion scattering is not evidenced in the copper and gold targets, because they are not fully ionized at the peak temperatures of 2.5 and 2.3 keV, respectively.

We now analyze the energy balance for heavier materials. Ion energy deposition can be split into drag heating due to Coulomb collisions and Ohmic heating due to return currents. The relative importance of Ohmic heating for ion energy deposition as a function of beam current density is shown in Fig. 11 for five different materials of increasing atomic number, from deuterium–tritium (DT) to gold. For small current densities of 10⁹–10¹⁰ A/cm², drag heating dominates and Ohmic heating can be neglected. For higher current densities, the Ohmic heating fraction increases dramatically, except in the case of DT, for which it remains negligible owing to the very low resistivity of this material. For instance, the Ohmic heating fraction for gold rises to about 30% at 10¹³ A/cm². It is interesting to note that for current densities around 10¹⁰–10¹¹ A/cm², although the resistivity of aluminum is lower than that of gold, its Ohmic heating fraction is higher. This can be explained by the higher stopping power of gold, which reduces the generation of resistive fields. In addition, the beam is compressed in the aluminum target by the azimuthal B-field, whereas such an effect is not as evident for copper and gold targets. This is an additional factor that further increases the Ohmic heating in the aluminum target.

To check the robustness of the proton transport in solid targets, 3D simulations of 5 MeV perfectly collimated protons with current density 10¹² A/cm² impinging on an aluminum target are performed. Figures 12(a)–12(d) show cross-sectional distributions of electron temperature at 20 ps. Four different transverse cuts are presented. It can be seen that with the proton transport in the target, the ion beams are gradually compressed by the resistive fields [see Figs. 12(a) and 12(b)] and are strongly compressed at z = 140 μm with an aluminum plasma temperature of 3 keV [see Fig. 12(c)]. Meanwhile, it can be clearly observed that the filamentary structure becomes more pronounced [see Fig. 12(b)], which may be due to filament instability. Dozens of filamentary structures with sizes of 1–3 μm appear in the transverse direction and gradually converge together [see Figs. 12(b) and 12(c)]. As the protons propagate farther, the target is almost fully ionized, and then the ion scattering start to take effect and ion beams spread into space [see Fig. 12(d)]. The electron temperatures on the four different cross-sectional distributions at the central axis in the 3D simulations are shown in Fig. 12(e) and are consistent with the 2D simulation results at the central axis shown in Fig. 10(b). It is important to note that although 2D simulations could provide a general understanding of the physics involved, they still suffer from limitations compared with 3D simulations. 2D simulations are always uniformly distributed in the transverse direction, whereas 3D simulations may have a certain distribution in the transverse direction, which may have some effects. For example, the intense proton beam heats up the entire target transversely in the 2D case, whereas in the 3D case it only heats up the cylindrical region.

The effects of resistive fields and stopping power modeling of the ion beam-plasma interaction have been examined using hybrid simulations. The experimentally validated BPS free stopping power model and the BC bound stopping power model have been applied. It is found that the self-generated resistive field can compress the ion beams, which is quite important for ion beam transport in a solid target. When the target is almost fully ionized, ion scattering effects start to come into effect, competing with compression by the resistive field and finally leading to divergence of the ion beams. Moreover, hybrid simulations show that Ohmic heating is important for higher current densities and materials with high resistivity, for which the energy fraction deposited as Ohmic heating can reach 20%–30%. It is found that the initial ion divergence can substantially modify proton energy deposition and subsequently affect proton transport in solid materials. Regarding the importance of self-generated fields in the proton fast ignition scenario, these can play an important role in proton beam cone tip interaction. The cone tip will be an intermediate- or high-Z material compressed by the shock coming from the target center. Resistive fields may scatter and slow down the proton beam, affecting the energy deposition in the dense DT. This will be studied in detail in our future research.

Acknowledgment. This work was supported by the National Natural Science Foundation of China (Grant Nos. 12005298, 12275356, 11774430, U2241281, and 12175309), Research Grant No. PID2022-137339OB-C22 of the Spanish Ministry of Education and Research, the Natural Science Foundation of Hunan Province (Grant Nos. 2021JJ40661 and 2022JJ30656), and a research project of the NUDT (Contract No. ZK19-25). The authors used the HPC resources and technical assistance from the BCS and CeSViMa centers of the Spanish Supercomputing Network.