Divergence of fast ions generated by interaction of intense ultra-high contrast laser pulses with thin foils

We propose an analytical model that analyzes the divergence of fast ion beams accelerated at the rear of thin foils irradiated with ultra-short intense laser pulses. We demonstrate the critical role played by the non-stationary character of the side components of the electric field, which is responsible for ion acceleration from the back of the foil. The model predictions are in very good agreement with 2D PIC simulations and with the experiments performed in the ultra-high-contrast regime as well.


Introduction
For more than 10 years, a great deal of scientific effort has been put by the scientific community into the domain of high-energy particle acceleration through intense laser pulses. The mechanism known as target normal sheath acceleration (TNSA) [1] is the main acceleration process for a very large part of currently available laser systems and has been mostly theoretically [2]- [4] as well as experimentally [5] studied. The interest this subject has attracted is due to the particular features (low divergence, duration and emissivity, and high energy and laminarity) these mega-electron volt ion bunches present and the consequent number of applications open to exploit them. Laser-driven ion beams have already been used for producing high-energy density matter [6] and probing dense plasmas [7], whereas proton therapy [8], isotope production for medical purposes [9] and fast ignition [10] are equally pointed out as potential applications, on the condition that some of the beam specifications be further improved. In this connection, the two main trends in current research on laser-driven ion acceleration concern the increase of maximum ion energy and the control of beam divergence. The latter shows practically the same variation, as a function of ion energy, in most of the published works, aside from experimental parameters [11]. A lot of attempts at controlling the spatial characteristics of ion bunches have been made acting on the target geometry [6,12,13] or during the ion propagation [14,15], whereas other works have pointed out the main role of the spatial shape of the electric field accelerating ions on the target's rear side [16,17]. Although ion divergence has been analyzed in theory and simulations [18]- [23] for micrometric thickness foils, so far there has neither experimental evidence, nor a clear analytical model able to explain the ion bunch divergence features using ultra-thin foils as the target. This appears to be indeed of great interest as the recent progress in the laser beam contrast ratio improvement devices [24,25] has made it possible to use such sub-micrometric targets and to get a significant increase in the maximum ion energy attainable.
In this paper, we analyze the divergence of fast ions accelerated from the rear of a thin foil by a short and intense high-contrast laser pulse. We will show how the non-stationarily of the transverse component of the accelerating electric field is responsible for the deviation of the traveling ions from the target normal. Actually, exploiting the two main features of this interaction regime (ultra-high contrast and ultra-thin target), we are able to explain and model this non-stationary character through the lateral expansion of the fast electron cloud during the ion acceleration and determine the time dependence of the electric field components. Our model shows remarkable agreement with experimental data and particle-in-cell (PIC) simulations, allowing us to estimate the divergence of fast ions produced by interaction of a femtosecond laser pulse with ultra-thin foil targets.

Theoretical model
In our model, we consider an ultra-short (<50 fs), p-polarized laser pulse impinging on a submicrometer thin foil at 45 • incidence, in a two-dimensional (2D) Cartesian geometry, and take into account the motion of the hot electron cloud through the target (along the z-axis) as well as its transverse expansion in the same target plane (y-axis).
Due to the extreme shortness of the pulse compared to the typical ionic motion timescales, we will assume that the ionic acceleration starts only at the end of the pulse [26]. Under these conditions, we can reasonably assume that the energy transfer between hot electrons and ions follows the adiabatic law: p e /n γ e = T e0 /n γ −1 e0 where p e and n e are, respectively, the pressure and density of hot electrons, T e0 and n e0 are their initial temperature and density, and the adiabatic parameter γ = (2 + N )/N where N is the number of degrees of freedom. Moreover, as the ion acceleration mechanism is a slow process compared to the hot electron motion in the foil, we can neglect the electron inertia or, in other words, neglect the electron mass in the standard force equation for the density and the mean electron velocity [27]. Considering only the propagation direction, the equation of the force then reduces to ∂ϕ ∂z = 1 en e ∂ p e ∂z where ϕ is the electric potential and e the electron charge. Using the adiabatic equation, we then easily obtain the hot electron density n e as a function of the potential ϕ of the ambipolar field as: n e (ϕ) = n e0 (1 + γ −1 γ eϕ T e0 ) (1/γ −1) . It is worth noting that this expression is correct even for relativistic electron temperature (T e > m e c 2 ). Therefore our model applies even to laser intensities exceeding 10 18 W cm −2 . For a 2D expansion, γ = 2 and the relationship between the hot electron density and the potential becomes linear and can be written as η e = 1 + ψ/2 where ψ = |e|ϕ/T e0 and η e = n e /n e0 . We will see that, through a suitable choice of ρ(y, z), this allows analytically solving the two-dimensional Poisson equation, ϕ = 4π(−ρ(y, z) + en e (ϕ)). As the cloud of hot electrons accelerated by the laser pulse circulates through the ultra-thin target and expands in the transverse direction, a positively charged spot is induced on both sides of the target, keeping most of the hot electrons in close proximity of this latter. The net positive charge density on the target is then given by the difference between the cold ion and the electron densities: ρ(y, z) = Z en i (y, z) − en e cold (y, z). We recall here the main assumptions of the model: the target is supposed to be extremely thin compared with the Debye radius for hot electrons and the laser beam contrast ratio high enough to neglect any amplified spontaneous emission or pre-pulse-driven plasma dynamics before the pulse peak. In that case, it is possible to approximate the net positive charge as ρ(y, z) = en e0 η(ξ )δ(ζ ), where ξ = y/r D and ζ = z/r D are the coordinates normalized to the Debye radius of hot electrons r D = T e0 /4π e 2 n e0 , δ(ζ ) is the Dirac function and η(ξ ) is the dimensionless function describing the transverse profile of the charge density. To mimic the dynamics of this latter, we simply suppose a temporal dependence of the charge density profile η(ξ, t) resulting from the normalized transverse scale of the hot electrons, spot l(t) = L(t)/r D that means η(ξ, t) = η(ξ, l(t)).
The Poisson equation in terms of electron density can now be written as The solution of this equation is known (see for example [29]) and it can be expressed as Here K 0 is the McDonald function (modified Bessel function of the second kind). The advantage of this method is the possibility to consider different transverse distributions of the charge density with various shapes-rectangular, Gaussian or parabolic-all with the same width l(t) (see figure 1). As the Poisson equation does not contain time derivatives, the time variable will appear as one of the parameters in the solution. Any profile distribution in figure 1 showing the same total charge (area) and the characteristic scale l(t) is denoted by the same color.

Dynamics of the charged region of the rear target surface
At t = 0 the solution of the Poisson equation corresponds to a quasi-neutral situation, and thus the electric field is absent at infinity for both ξ and ζ . For t > 0, the hot electron expansion will produce a flow of cold electrons from the target periphery to the charged region, frustrating the quasi-neutral state and leading to variation of the width l(t) (see figure 1). In other words, the tangential component of the electric field E y = β E ξ = β ∂ψ ∂ξ (here β = T eh /r D ) will pull the cold electrons from the spot boundary to its center in accordance with the Ohm law j cold = σ cold E ξ .
The electrostatic energy of the hot electrons depends on the chosen profile of charge density. We use in our calculations the following smooth profile: η(ξ, l(t)) = η e0 /πl(t)(1 + ξ 2 /l(t) 2 ) = η s , whose calculated electrostatic total energy was found to be the lower one for the set of functions we tested ( figure 1(b)). Here η e0 = 4π e 2 N eh T e0 where N eh ≈ Aε L D e T e is the total number of hot electrons per length, D e is the initial electron emission spot size, A is the average absorption coefficient and ε L is the laser radiation energy.
We estimate the conductivity of the solid plasma target as: σ cold ≈ ω pc 4π , neglecting the influence of the magnetic field because in our conditions the cold electrons gyro-frequency is lower than the plasma frequency ω pc = 4π e 2 n ec m e . Transverse dynamics of electron spot size can be obtained from the differential form of the charge conservation law for all particles in the foil at the point ξ = l: Here l f = L f r D is the normalized thickness of the foil, and j c = j cold /en e0 c is the dimensionless density of the cold electron current. It follows that the spot expansion is connected to the target conductivity: first, the spot dimension increase is due to the redistribution of the electron electrostatic energy over a larger area. At the beginning, the solution of equation (3) can then be approximated by a diffusion one. Later, when the expansion time becomes comparable with the effective time of collisions 1/ν ef , the diffusion law can no longer be applied as the expansion exponentially decays. According to that, the expansion law can be appropriately written as Here l 0 = L(0)/r D , κ ≈ a p (σ cold /ω ph )l f , a p ∼ 1 is a constant, which depends on the charge density profile, ν −1 ef is the effective collision time and τ = tω ph . To validate our analytical model, we compare its predictions with the results of our 2D PIC code (modified version of [28]) simulation. In figure 2, we show the calculation of the spot expansion of hot electrons for a thin Al target, as given by the model and the code.
The values of the hot electrons initial energy (T e0 = 1 MeV) and density (n e0 = 2.3 × 10 19 cm −3 ) and the total number of hot electrons inside the target diameter (N eh = 10 8 cm −1 ) that we used for the analytical model come from the simulation results. The values of N eh and T e0 have been used to calculate the constant κ for a charge profile distribution like η s . The agreement is quite satisfactory and this expression and parameters will be used for the calculations in the following.

Ion beam divergence
It is well known that the limitation of the transverse target size increases the ion energy [13,30].
In our case, the maximum transverse size is obtained from formula (14) as l m = l 2 0 + κ/ν ef . A target transverse size lower than l m allows the electron temperature to be increased, therefore 6 enhancing the ion energy compared to the case of an infinite size foil. On the other hand, decreasing the target size increases the ion angle of departure [13]. Thus, there is an optimal transversal target size maximizing the ion energy flux in a given solid angle centered on the target normal [31].
The number of ions accelerated by the ambipolar field is small compared to the total number of ions even in the case of an ultra-thin foil. So we can treat them as 'test' particles. Let us thus write the equations of ion motion in the potential ϕ(ξ, ς) in order to obtain the ion distribution as a function of the angle of departure: This equation can be solved assuming a zero starting speed and an arbitrary position of ions on the target surface:ζ | τ =0 = 0, ξ | τ =0 = ξ 0 . The solution gives the angle of ion departure from the region covered by the field, as a function of starting ion position for the given profile of charge distribution: One can write out the function of ions distribution across the departure angles as where N i is the number of ions, and n is the number of accelerated ions per unit of the target length. Let us note that formula (4) of the expansion of the spot of hot electrons weakly depends on the specific choice of the initial charge density profile, while the outlook of the distribution function reveals stronger dependence. For example, the use of rectangular-shaped or Gauss-like profiles results in higher slope of angular distribution, exceeding the slope of the numerical curve.
In order to estimate this angle at a large distance from the target, we replace the K 0 function in (2) by its asymptotic value at ς 1 and calculate this integral by the method of steepest descent. We obtain And, finally, from equation (8), we obtain θ θ (°) y 0 (µm) Here, H is a step-like Heaviside function. The angular distribution of ions can be depicted as the dependence of the modulus of the derivative dθ/dξ 0 on departure angle θ (see (7)). The angle of ion departure from the region covered by an electric field as a function of starting position of the ion, θ(ξ 0 ), for the charge profile distribution like η s is shown in figure 3. According to the latter, the approximated expression for the angular distribution we get by inserting η s in formulae (9) and (10) is in quite good agreement with more accurate numerical calculations.
One can see that the smooth distribution provides a more adequate description of the departure angle as a function of the ion position compared to the rectangular one. For both kinds of profile, the longitudinal field variation is only due to the hot electron density decrease, but for the rectangular profile, the density gradient does not vary in time. Hence, the transverse field component fades faster in the case of the smooth profile than for the rectangular one. As a consequence, the ions in the peripheral part of the charged spot will be preferentially accelerated in the longitudinal direction and have a small start angle (see the dark blue curve in figure 4); meanwhile, the ions that initially lie closer to the charged spot edge, where the density gradient is the biggest, will experience a larger deviation angle. The dependence of ion departure angle on ion position is consequently not monotonic and presents a local maximum (see figure 3), suggesting the non-laminar character of the propagation for some ion group.
Let us consider the dependence of ion departure angle on its energy. In our model, both the ion energy and its departure angle are determined by the starting position of the ion ξ 0 (see formulae (8) and (10)). Using these formulae, one can outline the dependence of θ(ε). It is worth noting that, unlike what happens physically, our model does not consider that part of the ionic population which originates from the target skin depth. As these ions belong to the low-energy side of the spectral distribution, we will consider the results of our model, using a smooth profile, as better suited for the high-energy ionic population (close to the cut-off energy).
In figure 4, the black curve represents the results of numerical calculation of θ(ε) for the above outlined parameters at t = 133 fs, whereas the blue line depicts the results of the model for the smooth charge density profile η s . One can see the correspondence for the smooth profile is observed only in the range of high energies, while in the low energy range, the model and the PIC calculation reveal opposite tendencies. In contrast, the rectangular-shaped profile gives larger angles at low energies (see the red line in figure 4), behavior that qualitatively corresponds to the numerical simulation results. Such a special feature of the rectangular density profile results from the large value of the tangential component of the electric field, which is found to be more than one order of magnitude higher than that of the smooth distribution. As the profile stands rectangular during the spot expansion, so does the high value of the tangential electric field component at the spot periphery, where the low-energy part of ions is produced. As a consequence, these ions leave the target with a large departure angle.
In the simulations, we see that the rectangular profile of the charge spot appears only at the initial stage of the interaction even for rectangular laser beam shape and later it transforms into some smoother distribution. Therefore the angular distribution given by the numerical simulations actually represents the combination of an initial rectangular distribution further followed by a smoother distribution. Once more, this is an intrinsic drawback of our model that cannot take into account the temporal evolution of the profile. Rectangular and smooth distributions have to be considered as two limiting cases. 9 contrast beam. The experiment was performed at the Saclay Laser Interaction Center Facility, using the UHI10 laser, which delivers 10 TW ultra-short pulses (65 fs) at 10 Hz repetition rate at a central wavelength of 790 nm. The intrinsic 10 6 contrast of the beam is raised to 10 10 thanks to a 'double plasma mirror' [32] (DPM), which preserves the spatial focal spot qualities but halves the laser energy. In the experiment the laser beam was focused to a spot size of 8 µm (full width at half maximum (FWHM)) using an off-axis f = 300 mm parabola, under an incidence angle of 45 • and p polarization, on thin CH foils with thickness varying between 0.1 and 100 µm. The laser intensity using the DPM was about 3 × 10 18 W cm −2 . The spectra of protons generated from the back side of the targets were recorded using a Thomson parabola spectrometer and a radiochromic film (RCF) based spectrometer. The latter is composed of a set of RCF dosimeter media (Gafchromic HD-810) filtered by aluminum foils of different and properly chosen thicknesses spanning from 0.8 to 40 µm. Each aluminum foil was placed in front of a couple of RCF placed one after the other. This allowed us to record selected proton energies ranging from about 450 keV to 1.9 MeV and from about 3 to 3.8 MeV, with an energy resolution of about 300 keV. Each couple of RCF, placed 19 mm from the target, collected the protons produced by one laser shot, before being automatically displaced and the shot was repeated on the following couple of RCF. This procedure, allowed by the high level of reproducibility of our system (about 5% rms observed variation on shot-to-shot maximum energy and particle flux), was repeated for six different aluminum foil thicknesses. Data collected through the RCF give an absolute estimation of the spectral and spatial distributions of the proton beam. The reliability of the RCF-obtained spectral distributions was verified by comparing them to the spectral distributions collected by the micro-charged plate-based Thomson parabola. In the present case, we focus our attention on the divergence of the proton beam at different energies. We scanned each exposed RCF and we measured the FWHM of the proton spot. Due to the oblique incidence of the laser pulse, the produced charge spots have an elliptical shape extended along the incidence plane. To compare the experimental data with the results of our axial-symmetry 2D model as well as the results of numerical simulations, in the following we rather denote as 'divergence' the average value of the divergences along the incidence plane and the normal direction. In figure 5, we report the experimental values of the divergences obtained using 800 nm thick CH target and the related results of the model and of the 2D PIC code.
There is a marked difference between our data and the almost linear behavior of divergence as a function of energy usually observed in the past (see for example [11]). In spite of that, we found large divergence values for low-energy protons, followed by a kind of plateau for intermediate energy values with an abrupt decrease of divergence close to the cut-off energy. The small value of divergence in the plateau-like part of the curve means that, for related energies, the proton beam is well collimated.
In the cases [11], the sheath field accelerates protons for a long time [33] and proton angular spread gradually increases with a decrease in their energies. In the short, high-contrast laser pulse case, most energetic protons are accelerated by the maximal electric field, which continues longer than the pulse duration. At this time, the accelerating field front shows a curvature, resulting in larger angular spread for higher energetic protons [34].
In the light of what was discussed before, the observed divergence trend can be explained as the combined action of an initial rectangular-like spot profile (responsible for the low-energy part of the spectrum) evolving towards a smooth-like profile. The picture shows quite good agreement with the results of the PIC simulations. We used the numerical results we obtained as input parameters for the analytical model, testing different density profiles. The best agreement θ θ (°) ε (MeV) Figure 5. Experimental dependence of proton beam divergence on proton energy in two perpendicular directions. Open circles are the results of averaging between the two curves, in the incidence plane and in the perpendicular one. The line represents the result of the analytical model and the squares are the results of our simulations.
between the model and the experimental data (line in figure 5) was obtained using the smooth density profile with a charge spot size of 5 µm and a temperature of T e0 = 1.2 MeV, as given by the PIC simulations (the initial density of the target was assumed to be 6 × 10 22 cm −3 and the carbon average charge Z = 4). Once again we observe good agreement between the model and the experiment, all along the proton energy range, showing the correctness of our analytical approach.

Conclusion
In conclusion, we have presented an analytical model and simulations confirmed by experimental results allowing us to find the dependence of the hot electron spot size at the rear of thin (less than the Debye radius of the hot electron) foils with time. We show how the angular distribution of accelerated ions and its dependence on ion energy are functions of the characteristics of the non-stationary transversal electron density profile. In particular, we demonstrate how the presence of a smooth electron density profile can lead to a non-laminar ion expansion. We estimated the maximum target size for efficient ion production on interaction of a short laser pulse with a size-limited target. Based on our simulations and experiments, we developed an analytical model of ion beam divergence for thin foil targets. The results of our model are in good agreement with 2D PIC simulations and with the already published experimental results.