Ion energetics in electron-rich nanoplasmas

For three examples involving different cluster sizes and laser parameters, figure 1 shows time-resolved key quantities of the ionization process: the average ion charge q_av, the number n_p of nanoplasma electrons per atom inside the cluster, the relative cluster radius R/R₀ of the expanding cluster (taken as the distance of the most distant ion from the cluster center-of-mass, with R₀ being the initial cluster radius) and, for a comparison of time scales, the oscillating laser electric field. The three examples represent three categories of incomplete outer ionization:

• (A)  
  Vertical moderate outer ionization for short pulses, Xe₁₀₆₁, I_M = 10¹⁶ W cm⁻², τ = 10 fs (figure 1(a)). The nanoplasma electron population is in the range 1/2–3/4 of the average ion charge (n_p/q_av = 0.73 in this example), and the cluster does not expand appreciably during the laser pulse. Depending on the cluster size, this category is found for moderate laser intensities 10¹⁵–10¹⁷ W cm⁻².
• (B)  
  Non-vertical low outer ionization for short weak pulses, Xe₁₀₆₁, I_M = 10¹⁴ W cm⁻², τ = 10 fs (figure 1(b)). Outer ionization is very low, only a few per cent of the average ion charge. Despite the negligible cluster expansion during the laser pulse, outer ionization is non-vertical, which will become apparent in the next section of this paper.
• (C)  
  Non-vertical moderate to high outer ionization for long weak pulses, Xe₆₀₉₉, I_M = 5 × 10¹⁴ W cm⁻², τ = 230 fs (figure 1(c)). The cluster expansion during the laser pulse is large; at the end of the laser pulse, t ≈ τ, the cluster radius is increased by a factor of 10. For the long pulse of τ = 230 fs, outer ionization is very efficient via resonance heating [3, 49]. In this example, the cluster size and laser intensity were chosen such that the outer ionization level n_p/q_av = 0.69 at the end of the trajectory is comparable to case (A). (As a general rule, outer ionization increases with increasing laser intensity and with decreasing cluster size.)

Zoom In Zoom Out Reset image size

Before analyzing the ion kinetic energies, it is instructive to consider the total energy of electron-rich nanoplasmas. Figure 2 exhibits for the examples (A) and (B) the time-resolved total energy E_tot of the entire system consisting of n ions and the total number nċn_e of electrons inside and outside the ionic framework (n_e = n_p + n_oi, with n_oi being the number of electrons per atom outside the ion framework). Shown also is the partition of E_tot, E_tot = T_n + T_e + V_tot, with T_n and T_e being the ion and electron kinetic energy, respectively, and V_tot the total potential energy as the sum of ion repulsion, electron repulsion and ion–electron attraction. In both trajectory calculations all electrons were included; the standard procedure of removing distant electrons >10R (see section 2) was not applied. E_tot is measured with respect to a total decomposition of the nanoplasma into single particles. E_tot is not constant after the termination of the laser pulse due to the continued formation of ion charges and nanoplasma electrons by EII, as the expended ionization energy is not counted in this balance of E_tot.

Zoom In Zoom Out Reset image size

In the case (A) (figure 2(a)) E_tot is positive, allowing in principle for a total decomposition of the plasma into single ions and electrons. The largest contribution to E_tot at long times is the ion kinetic energy T_n. Another large part of E_tot is carried away in terms of electron kinetic energy; more than 96% of the long-time T_e value is attributed to electrons that were stripped by outer ionization. The total potential energy V_tot converges to a negative value despite the ongoing cluster expansion, indicating ion–electron trapping. Further analysis showed that 40% of all nanoplasma electrons remain trapped by ions at long times, orbiting the ions at distances r_cut < 2 Å. Considering only the effective ion charges q_i^(eff) , i.e. the bare ion charges q_i reduced by the number of trapped electrons, and the remaining n⋅n_e^(free) free electrons, the resulting potential energy V'_tot (dotted line in figure 2) becomes positive and converges to zero at long times. (Regarding notation, we shall distinguish energies, which were calculated from effective ion charges and free electrons, by a prime.) Of particular interest is the potential energy V'_ions which is available for the conversion into ion kinetic energy T_n. When the spatial charge distribution is spherical and the nanoplasma sufficiently diluted, V'_ions is approximately given by (in atomic units)

$$\begin{equation}V'_{{\rm{ions}}} \approx \mathop \sum \limits_i \mathop \sum \limits_{j(r_j < r_i )} \left( {\frac{{q_i^{({\rm{eff}})} q_j^{({\rm{eff}})} }}{{r_{ij} }}} \right)-\mathop \sum \limits_i \mathop \sum \limits_{k(r_k < r_i )} \left( {\frac{{q_i^{({\rm{eff}})} }}{{r_{ik} }}} \right),\end{equation} \tag{ 5 }$$

where k in the second double sum runs over all n⋅n_p^(free) free nanoplasma electrons. For comparison, V_ions is also included in figure 2, i.e. the corresponding quantity with the full ion charges q_i and nċn_p nanoplasma electrons. While the long-time value of V_ions is negative, V'_ions assumes a small positive value and goes to zero at long times, driving the ion acceleration.

In the case (B) (figure 2(b)) E_tot converges to a negative value at long times; total decomposition of the nanoplasma into single particles is therefore precluded. Some 70% of the nanoplasma electrons were found to be trapped. V'_tot and V'_ions assume small positive values.

Another interesting issue of the energy analysis pertains to the question of the dominant driving force of the cluster expansion, CE or hydrodynamic expansion (HE), being attributed to the conversion of total potential energy and electron kinetic energy, respectively, into ion kinetic energy [36]. Under the condition that the total energy E_tot of the system is constant in the course of the trajectory, the CE part is given by the difference ΔV_tot = V_tot(t_L) − V_tot(t_I) of the total potential energy of the system and the HE part by the difference ΔT_e = T_e(t_L) − T_e(t_I) of the electron kinetic energy, t_I and t_L being instants close to the beginning and at the end of the trajectory. In the ideal case, t_I is chosen at a time when the inner ionization and energy absorption of the electrons has terminated and E_tot has reached its final constant value. On the other hand, the kinetic energy attained by the ions at t < t_I cannot be assigned to CE and HE, so that t_I should be chosen as close as possible to the beginning of the trajectory. In practice, t_I is a trade-off between both requirements. For the same reason, the energy analysis is not meaningful for long pulses (case (C)).

Table 1 contains the results of the analysis for examples (A) and (B). For the example (A) we chose t_I = τ = 10 fs, when the nanoplasma build-up is nearly complete and the variation of the total energy in the course of the trajectory, E_tot(t_L) − E_tot(t_I) = 93 keV, is small compared to the final ion kinetic energy, T_n(t_L) = 2079 keV. Also, at t = 10 fs the non-assignable part of the ion kinetic energy is still relatively small, T_n(t_I) = 155 keV. With this choice of t_I and t_L, ΔV_tot = 1702 keV and ΔT_e = 317 keV, so that the relative contributions of CE and HE to the sum ΔV_tot + ΔT_e + T_n(t_I) are 78 and 15%, respectively. The close agreement of the sum ΔV_tot + ΔT_e + T_n(t_I) = 2174 keV with T_n(t_L) = 2079 keV reflects the good accuracy of the energy partitioning. Clearly, example (A) pertains to the CE regime.

The analysis of example (B) is fraught with much larger inaccuracies. Until t = τ, only 73% of the inner ionizations occur. We therefore chose t_I = 112 fs; when T_e becomes maximal, E_tot(t_L) − E_tot(t_I) drops to 20 keV and T_n(t_I) = 2.6 keV is still relatively small. With this choice of t_I one obtains ΔV_tot = 30 keV and ΔT_e = 36.3 keV, corresponding to contributions of 44% and 53% CE and HE, respectively. Although ΔV_tot + ΔT_e + T_n(t_I) differs from T_n(t_L) by 20 keV and E_tot(t_L) − E_tot(t_I) = 20 keV, indicating a large error of the energy partitioning which is 40% of T_n(t_L), the trend is visible that a large, if not a major part of the cluster expansion, stems from HE.

Proceeding from an analysis of the total energy of the system to the individual ions, their long-time kinetic energies can, in principle, be obtained directly from a trajectory calculation. In practice, because of computational workload, trajectories are mostly not long enough, so that the potential energy residues have to be estimated. There are several methods for an energy correction:

• (1)  
  A pairwise interaction formula based on V'_ions , equation (5). The corrected energy T_i^(corr) of ion i is given by its instantaneous kinetic energy T_i and the sum over all Coulomb interactions with the effective ion charges q_j^(eff) and free nanoplasma electrons k closer to the cluster center than ion i:

  $$\begin{equation}T_{i}^{({\rm{corr}})} = T_{i} +q_i^{({\rm{eff}})} \mathop \sum \limits_{j(r_j < r_i )} \left( {\frac{{q_j^{({\rm{eff}})} }}{{r_{ij} }}} \right)-q_i^{({\rm{eff}})} \mathop \sum \limits_{k(r_k < r_i )} \left( {\frac{1}{{r_{ik} }}} \right).\end{equation} \tag{ 6 }$$
• (2)  
  A central charge formula with effective ion charges. The energy correction for ion i is a single pairwise Coulomb term of the effective charge q_i^(eff) and a charge Q_i in the cluster center, which is the sum of all the ion and electron charges closer to the cluster center than ion i:

  $$\begin{equation}T_i^{({\rm{corr}})} = T_i +\frac{{Q_i q_i^{({\rm{eff}})} }}{{r_i }},\end{equation} \tag{ 7 }$$

  $$\begin{equation}Q_i = \mathop \sum \limits_{j(r_j < r_i )} (q_j )-\mathop \sum \limits_{k(r_k < r_i )} (1).\end{equation} \tag{ 8 }$$
• (3)  
  The pairwise interaction formula, equation (6), with the full ion charges q_i and the full number nċn_p of nanoplasma electrons k:

  $$\begin{equation}T_i^{({\rm{corr}})} = T_i +q_i \mathop \sum \limits_{j(r_j < r_i )} \left( {\frac{{q_j }}{{r_{ij} }}} \right)-q_i \mathop \sum \limits_{k(r_k < r_i )} \left( {\frac{1}{{r_{ik} }}} \right).\end{equation} \tag{ 9 }$$

The energy correction corresponds to the potential energy V_ions and was used in the previous work [21] where electron trapping was mostly insignificant.

Each of the three energy corrections presumes a spherical charge distribution inside the nanoplasma. In practical calculations, smoothed Coulomb potentials for the ion–electron and electron–electron interactions were used in equations (5), (6) and (9), in accordance with the trajectory calculations [9, 16–18, 21, 22].

The Xe₁₀₆₁ cluster is still small enough to check the three energy correction formulae against long-time ion kinetic energies, where T_i^(corr) ≈T_i . For the outer ionization case (A), figure 3(a) shows the maximum corrected and uncorrected ion energy, E_M and T_M, respectively, as a function of time. The three correction formulae lead to nearly identical results and are nearly converged to the long-time T_M value already at t = 200 fs. For the average ion energy (figure 3(b)), the values corrected by the pairwise interaction formula with effective ion charges and by the central charge formula are identical, converged already at t = 20 fs and in close agreement with the long-time 〈T〉 value. In contrast, the pairwise interaction formula without trapping gives too low values. The reason is that using bare ion charges in the energy correction formula (equation (9)) does not take into account the Coulomb repulsion between trapped electrons and other electrons which are closer to the cluster center. As the trapped electrons form a subsystem together with their host ion, the energy correction must involve all interactions of the subsystem with other ions and electrons. When the interparticle distances are large compared to the distances inside the ion-trapped electron subsystem, it suffices to reduce the ion charges rather than considering the repulsion of the subsystem electrons with other nanoplasma electrons explicitly. For the maximum ion energy (figure 3(a)), the pairwise interaction formula with bare ion charges and the full number of nanoplasma electrons (equation (9)) works well, because the maximum energy stems from an ion near the cluster surface, where trapping does not play a role for these laser parameters.

Zoom In Zoom Out Reset image size

In all calculations an ion–electron cutoff distance of 2 Å was used; the ion energies turned out to be insensitive to the particular choice of the cutoff distance. The residence times of the trapped electrons increase with the dilution of the nanoplasma. A cursory analysis revealed that for the Xe₁₀₆₁ cluster irradiated by a pulse of I_M = 10¹⁶ W cm⁻², τ = 10 fs, the average residence time is of the order of 300–400 fs at t = 1 ps.

Figures 3(c) and (d) exhibit the maximum and the average ion energy for outer ionization category (B), the Xe₁₀₆₁ cluster irradiated by a laser pulse of I_M = 10¹⁴ W cm⁻², τ = 10 fs. The pairwise interaction formula with trapping (equation (6)) and the central charge formula (equations (7) and (8)) give nearly identical results for E_M, whereas the pairwise interaction formula without trapping (equation (9)) results in unphysical fluctuations of E_M (figure 3(c)). As equation (9) sums ion–electron interactions r_k < r_i, the number of summands oscillates with the instantaneous coordinates of the trapped electrons orbiting their host ion, indicating that for short weak pulses electron trapping plays a role also for surface ions. For the average ion energies (figure 3(d)), the pairwise interaction formula with trapping and the central charge formula converge to the same final value. The pairwise interaction formula without trapping (equation (9)) fails completely and gives negative average ion energies. Comparing the ion energetics of cases (A) and (B), not only the average and maximum ion energies are much lower in case (B). Rather, the striking difference lies in the slow convergence of ion energies, in particular of the average ion energy, in case (B). Whereas in case (A) the average ion energy is converged already at 20 fs after the temporal peak of the laser envelope function, in case (B) convergence is achieved only after ≈5 ps. This point will be analyzed in section 3.2 utilizing the lychee model.

In figures 3(e) and (f), the time-dependent maximum and average ion energies are presented for the example (C), the Xe₆₀₉₉ cluster irradiated by a 230 fs pulse at I_M = 5 × 10¹⁴ W cm⁻². The energies obtained by the pairwise interaction formula with and without trapping (equations (6) and (9)) and by the central charge formula (equations (7) and (8)) are nearly identical, as trapping is prevented by the high electron kinetic energies attained during resonance heating. The corrected average and maximum ion energies, 〈E〉 and E_M, pass a maximum during the resonance of the nanoplasma electrons with the laser field (t = 140–180 fs). During the resonance, the electron population within the ion framework is lowered (cf figure 1(c)), so that the energy correction formulae overestimate the ion energies. After the resonance situation, the electrons are decelerated by the laser field and partly return from the vicinity of the cluster (<10R) into the ion framework. The ion energies are nearly converged at t ≈ τ, when the laser electric field has abated.

Our calculated single-trajectory results are in reasonably good agreement with the results of Holkundkar et al [23]. Their maximum ion energies for 25 fs infrared laser pulses are 7 keV for Xe₂₁₇₁ at I_M = 10¹⁵ W cm⁻² and 77 keV for Xe₆₁₀₀ at 2 × 10¹⁶ W cm⁻²; our results are 5.2 and 76 keV, respectively. For the Xe₆₁₀₀ cluster irradiated by a 125 fs pulse at I_M = 10¹⁶ W cm⁻², Holkundkar et al [23] obtained 〈E〉 = 58 keV and E_M = 100 keV; our results are 〈E〉 = 56 keV and E_M = 121 keV. For a Xe₆₁₄₀ cluster irradiated by a 100 fs pulse at I_M = 10¹⁶ W cm⁻², Petrov and Davis [13] have simulated an average ion energy of 68 keV, which is considerably higher than our value of 38 keV obtained for the Xe₆₀₉₉ cluster.

Erk et al [50] have recently published experimentally obtained ion kinetic energy distributions (KEDs) for xenon and argon clusters. To make contact with their experimental results, the simulated single-trajectory KEDs must be doubly averaged over the cluster size distribution and over the laser intensity profile [22]. For the cluster sizes we assume a log-normal distribution with an FWHM equal to the modal cluster size [51]. For the laser intensity we take a Gaussian beam intensity profile [52]

$$\begin{equation}I_{{\rm M}} \left( {r,z} \right) = I_{{\rm M}}^{{\rm{max}}} \frac{{z_{\rm{R}}^2 }}{{z_{\rm{R}}^2 +z^2 }}\,{\rm{exp}}\!\left[ {-\left( {\frac{{2r^2 }}{{w_0^2 }}} \right)\left( {\frac{{z_{\rm{R}}^2 }}{{z_{\rm{R}}^2 +z^2 }}} \right)} \right],\end{equation} \tag{ 10 }$$

where I_M^max is the spatially and temporarily overall peak intensity at the center of the focus volume, w₀ the beam waist radius and z_R = πw₀² /λ is the Rayleigh length, with λ = 800 nm being the laser wavelength. The cluster beam, with a radius r_B, is directed along the r-axis, which is perpendicular to the focus offset axis z. The relation between r_B and z_R determines which part of the intensity profile is relevant for intensity averaging. For r_B ≫ z_R, the averaging over the entire three-dimensional (3D) laser focus volume has to be carried out; for z_R < r_B, only a slice of the laser intensity profile around z = 0 is involved, so that equation (10) reduces to a 2D intensity profile

$$\begin{equation}I_{{\rm M}} \left( r \right) = I_{{\rm M}}^{{\rm{max}}} \,{\rm{exp}}\!\left( {-\frac{{2r^2 }}{{w_0^2 }}} \right).\end{equation} \tag{ 11 }$$

Compared to the 3D profile, in the 2D case the relative weight of low intensities is considerably lower. In the paper by Erk et al [50], r_B = 0.75 mm, but z_R is not specified. However, from experience z_R can be assumed to be ≈2.5 mm [53], nearly realizing the 2D scenario of equation (11). The KEDs in the experiment of Erk et al were obtained for $I_{{\rm M}}^{\max } = 5 \times 10^{15} \,{\rm{W}}\,{\rm{cm}}^{-2}$ and τ = 35 fs. For the double averaging in [22], we used ion energies from trajectories in the cluster size range 55 ⩽ n ⩽ 4093 and intensities 10¹⁴ W cm⁻² ⩽ I_M ⩽ 5 × 10¹⁵ W cm⁻². Erk et al give KEDs for two average xenon cluster sizes, 〈n〉 = 200 and 700, which are in our simulated cluster size range. The experimental maximum ion energy for 〈n〉 = 200 is about 22 keV and for 〈n〉 = 700 27 keV. Applying the same definition of the maximum ion energy as in the experiment, that is, when the KED drops to the given fraction of its maximum value (2 × 10^–4 for 〈n〉 = 200 and 4 × 10^–5 for 〈n〉 = 700), the simulated maximum ion energies on the basis of a 2D intensity profile are 11 keV for 〈n〉 = 200 and 18 keV for 〈n〉 = 700. The corresponding maximum ion energies on the basis of the 3D intensity profile are much lower, 7 keV for 〈n〉 = 200 and 12 keV for 〈n〉 = 700. That is to say, all simulated values are considerably lower than the experimental values, but better agreement with the experiment is obtained for the 2D scenario. A more detailed analysis of the experimental KEDs will be given in a forthcoming paper. Trying to assess our simulation results for the ion energetics, one may expect in general that our trajectory simulation model tends to underestimate ion energies, since the neglect of the excitation-autoionization and recombination-autoionization EII channels (see section 2) causes too low ion charges. For 〈n〉 = 700, Erk et al observed a maximum charge state q_M = 19, whereas our simulations give values q_M = 13 for Xe₁₀₆₁ and q_M = 14 for Xe₂₁₇₁ (I_M = 5 × 10¹⁵ W cm⁻², τ = 35 fs). The deviations between simulation and experiment may be partly due to the uncertainties of the experimental cluster size and laser intensity distributions. As pointed out by Schultze et al [54], the assumption of a Gaussian beam intensity profile is a significant oversimplification.

Insight into the ion energetics can be gained by using the electrostatic lychee model [14, 21]. The model presumes that (i) inner and outer ionizations are vertical, (ii) all nanoplasma electrons are confined to the cluster interior and form, together with the ions, a neutral sphere where the sum of ion and electron charges is equal, so that (iii) only ions in the electron-free cluster periphery participate in the Coulomb explosion. Knowing the initial radius R₀ and the homogeneous initial density ρ of a spherical cluster, the lychee model gives analytical expressions for the central charge Q(r) as a function of the distance r from the cluster center, for the radius R_p of the neutral interior cluster sphere, for the ion energies E(r), for the average and maximum ion energies 〈E〉 and E_M and for the ion kinetic energy distribution (KED) P(E) [21]. While in previous work [21] the lychee model was expressed for an average ion charge q_av, an obvious improvement is extension to a charge distribution, in particular since E_M and therefore the width of P(E) are determined by the maximum ion charge q_M at the cluster surface. The extension to a multiple-charge state lychee model is straightforward, for simplicity presuming that the ions of each charge state are homogeneously distributed over the cluster with a density

$$\begin{equation}\rho _q = \rho \frac{{n_q }}{n},\end{equation} \tag{ 12 }$$

where n_q is the number of ions of charge q and n is the total number of ions. The resulting expressions for R_p, Q(r) and 〈E〉 are identical with those of the single-charge state lychee model (all energies in the equations are given in atomic units)

$$\begin{equation}R_{\rm{p}} = R_0 \left( {\frac{{n_{\rm{p}} }}{{q_{{\rm{av}}} }}} \right)^{1/3} ,\end{equation} \tag{ 13 }$$

$$\begin{equation}Q(r) = \frac{{4\pi }}{3}\left( {r^3 -R_{\rm{p}}^3 } \right)\rho q_{{\rm{av}}} \quad (r \geqslant R_{\rm{p}} ),\end{equation} \tag{ 14 }$$

$$\begin{equation}\langle E\rangle = \frac{{4\pi }}{5}\rho q_{{\rm{av}}}^2 R_0^2 \left[ {1-\frac{5}{2}\frac{{n_{\rm{p}} }}{{q_{{\rm{av}}} }}+\frac{3}{2}\left( {\frac{{n_{\rm{p}} }}{{q_{{\rm{av}}} }}} \right)^{5/3} } \right].\end{equation} \tag{ 15 }$$

For the energy E_q(r) of an ion of charge q at a distance r from the cluster center, one obtains

$$\begin{equation}E_q (r) = \frac{{Q(r)q}}{r},\end{equation} \tag{ 16 }$$

and for the maximum ion energy

$$\begin{equation} E_{{\rm M}} = \frac{{Q(R_0 )q_{{\rm M}} }}{{R_0 }}.\end{equation} \tag{ 17 }$$

The energy distribution function P(E) is the superposition of the partial distribution functions P_q(E) of each charge state

$$\begin{equation}P\left( E \right) = \mathop \sum \limits_q \left( {\frac{{\rho _q }}{\rho }P_q \left( E \right)} \right),\end{equation} \tag{ 18 }$$

with

$$\begin{equation}P_q \left( E \right) = \frac{3}{{E_{{\rm M},q} }}\frac{{\frac{{r^2 }}{{R_0^2 }}\left( {1-\frac{{R_{{\rm p}}^3 }}{{R_0^3 }}} \right)}}{{2\frac{r}{{R_0 }}+\frac{{R_{{\rm p}}^3 }}{{R_0^3 }}\frac{{R_0^2 }}{{r^2 }}}},\end{equation}\noindent \tag{ 19 }$$

where E_M,q = E_q(R₀) is the maximum ion energy for charge state q; cf equations (16) and (17). The P_q(E) and, consequently, also P(E) are normalized to 1 − n_p/q_av

$$\begin{equation}\int_0^{E_{{\rm M},q} } {P_q (E)\,} {\rm{d}}E = 1-\frac{{n_{\rm{p}} }}{{q_{{\rm{av}}} }},\end{equation}\noindent \tag{ 20 }$$

excluding the neutral lychee core, which is consistent with the average ion energy 〈E〉, equation (15). Reasonably, one applies equations (12)–(20) to the effective ion charges and the free nanoplasma electron population rather than to the full ion charges and total nanoplasma electron population.

Applying the lychee model to the analysis of the slow convergence of the average ion energy in case (B), we note that the expression for 〈E〉, equation (15), contains powers of n_p/q_av and q_av² in the prefactor. With x = n_p/q_av and the average net charge per atom, q_net = q_av − n_p, equation (15) can be recast:

$$\begin{equation}E = \frac{{4\pi }}{5}\rho R_0^2 q_{{\rm{net}}}^2 \frac{{1-\frac{5}{2}x+\frac{3}{2}x^{5/3} }}{{(1-x)^2 }}.\end{equation}\noindent \tag{ 21 }$$

To find out to what extent 〈E〉 is dependent on x = n_p/q_av, in figure 4 the fractional polynomial in equation (21), which is equivalent to the relative average ion energy 〈E〉(x)/〈E〉(0) for constant q_net,

$$\begin{equation}\frac{{1-\left(\frac52x+{\frac32}x\right)^{{5/3}} }}{{(1-x)^2 }} = \frac{{E(x)}}{{E(0)}},\end{equation} \tag{ 22 }$$

is plotted against n_p/q_av. Note that 〈E〉(0) is the average ion energy for the lower limit n_p/q_av = 0, when all nanoplasma electrons are trapped. The upper limit, n_p/q_av = 1, where the total cluster charge is zero and 〈E〉(n_p/q_av)/〈E〉(0) = 5/6, is approached by outer ionization category (B). Large clusters can come arbitrarily close to this limit. As shown by figure 4, the fractional polynomial, equation (22), varies at most by 17% over the entire n_p/q_av range, so that 〈E〉 depends mainly on q_net. Since q_net is directly related to the total cluster charge Q(R₀)

$$\begin{equation}Q(R_0 ) = \frac{{4\pi }}{3}\rho R_0^3 \left( {q_{{\rm{av}}} -n_{{\rm p}} } \right),\end{equation} \tag{ 23 }$$

we conclude that the average ion energy depends mainly on the total cluster charge rather than on the charge distribution inside the cluster, and Q may be a key parameter for understanding the convergence behavior of the ion energetics and the underlying physical processes.

Zoom In Zoom Out Reset image size

For the Xe₁₀₆₁ cluster, I_M = 10¹⁶ W cm⁻², τ = 10 fs (outer ionization category (A)), in figure 5 the central charge Q(i) (equation (8)) experienced by each ion i is plotted against the ion number i, where the ions are ordered according to their distance from the cluster center. In this way, figure 5 contains information on the charge distribution inside the cluster as well as on the total cluster charge; Q(1061) represents the central charge experienced by a surface ion. Q(i) is given for t = 0.02, 0.5 and 3.5 ps. For all the considered times, Q(i) reflects a typical spatial lychee charge distribution, low Q values in the cluster interior and a steep rise toward the cluster periphery. Q(1061) is practically converged at t = 0.5 ps. At t = 20 fs, when 〈E〉 is already converged but E_M still by 1.5 keV above its final value of 16.5 keV (figure 3(a)), Q(1061) exceeds its final value by 20%. Apparently, the lowering of Q(1061) for t > 20 fs is due to the rapid expansion of the cluster periphery into the surrounding electron cloud.

Zoom In Zoom Out Reset image size

For the outer ionization category (B), the spatial dependence and time evolution of Q(i) are very different from that of category (A). In figure 6(a), the Q(i) function is given for t = 0.02, 0.2, 1, 5, 10 and 15 ps. Q(i) assumes much lower values than in case (A) and is not converged until at least t = 10 ps. To elucidate the different roles of ions and electrons, in figure 6(b) Q(i) is calculated from the reduced ion charges only, i.e. excluding the free nanoplasma electrons. When the free nanoplasma electrons are excluded, not only are the Q(i) values much higher, as expected, but Q(i) also converges much faster. This provides additional evidence that the slow convergence of the total Q(i) function (figure 6(a)) is caused by a gradual loss of nanoplasma electrons, a continuous evaporation of nanoplasma electrons from the ion framework on the picosecond time scale. The consequence of the slow electron evaporation on the ion energetics is that the ions are accelerated by a variable central charge Q, precluding a meaningful estimate of the ion potential energy residues from short trajectories. The electron evaporation causes the outer ionization to be highly non-vertical, despite a negligible cluster expansion during the short laser pulse. A comparison of figures 3(d) and 6(a) shows that 〈E〉 is earlier converged (around t = 5 ps) than Q (not before t = 10–15 ps), reflecting that the effect of the electron evaporation on the ion energetics fades away at a late stage of the cluster expansion.

Zoom In Zoom Out Reset image size

For a quantitative comparison of the KEDs from the trajectory calculations and those of the lychee model, it must be taken into account that in the trajectory calculations, P(E) is approximated by a histogram h_k of ion energies of finite-energy intervals ΔE

$$\begin{equation}P\left( {\left( {k+\frac{1}{2}} \right)\Delta E} \right) \approx \frac{{h_k }}{{\Delta E}},\quad k \cdot \Delta E \leqslant E < (k+1) \cdot \Delta E.\end{equation} \tag{ 24 }$$

Unlike in the lychee model, in the trajectory calculations P(E) is normalized to 1, as well as the ions in the nearly neutral lychee core participate in the cluster expansion and contribute to the h_k of the lowest-energy interval. Therefore, to be comparable with the KED of the trajectory calculation, the constant term n_p^(free) /q_av^(eff) δ _k has to be added to h₀ of the lychee model

$$\begin{equation}\frac{{h_k }}{{\Delta E}} = \frac{1}{{\Delta E}}\left[ {\int_E^{E+\Delta E} {P_{{\rm{lychee}}} \left( {E'} \right)\,} {\rm{d}}E'+\frac{{n_{\rm{p}}^{({\rm{free}})} }}{{q_{{\rm{av}}}^{({\rm{free}})} }}\delta _k +\frac{{n_0 }}{n}\left( {1-\frac{{n_{\rm{p}}^{({\rm{free}})} }}{{q_{{\rm{av}}}^{({\rm{eff}})} }}} \right)\delta _k } \right].\end{equation}\noindent \tag{ 25 }$$

The last term in equation (25), n₀/n(1 − n_p^(free) /q_av^(eff) ) ⋅ δ _k,/ accounts for ions that have reached charge neutrality by trapping the corresponding number of electrons, with n₀ and n being the number of neutral atoms and the total number of atoms in the cluster, respectively.

For a numerical application of equation (25), with P_lychee(E) from equations (18) and (19), n_p^(free) and the ion charge state abundances $\left\{ {n_q^{({\rm{eff}})} } \right\}$ are taken from trajectory calculations. The question arises then at which instant of the trajectory, n_p^(free) and $\left\{ {n_q^{({\rm{eff}})} } \right\}$ should be taken, as the charge build-up in the cluster is not completely vertical even for outer ionization case (A). One possible choice is at a comparatively long time t_L of the trajectory, when the ion energies are converged. A better alternative seems to be the charge distribution, which is dominant for the cluster expansion; we therefore take n_p^(free) and $\left\{ {n_q^{({\rm{eff}})} } \right\}$at the instant t_acc when the average ion acceleration reaches its maximum, t_acc = 20, 200 and 110 fs for examples (A)–(C), respectively. For the example (A) the instant t = 20 fs practically coincides with the termination of outer ionization. In case (B), the delay of t_acc to 200 fs, which is long after the termination of the laser pulse, reflects the slow charge build-up in the cluster by electron evaporation. In case (C), the instant of maximum ion acceleration coincides roughly with the enhanced outer ionization at the nanoplasma electron resonance.

Figure 7 shows the distributions of effective and bare ion charges at the maximum ion acceleration times t_acc. For comparison, the corresponding ion charge distributions at long times t_L have been included. Using the abundances of effective ion charges and free nanoplasma electron populations at t_acc, the KEDs calculated by the lychee model are shown in figure 8, together with the corresponding results of the trajectory calculations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In case (A) (figure 8(a)), the KEDs of the trajectory calculation and of the lychee model are in overall good agreement and the average ion energy is in perfect agreement. The overestimate of E_M by the lychee model by 3 keV is a non-vertical effect, as the cluster expansion sets in before the charge build-up is complete. The low-energy component h₀, which constitutes the maximum of the KED, is reproduced by the lychee model within 35%.

The underlying effective ion charge abundances of case (A) (figure 7(a)), at t_acc = 20 fs and t_L = 3.5 ps, differ mainly for the low ion charges 0–4, whereas the highest ion charge q = 13, which defines the width of the KED, is practically unaffected. Effective ion charges of 0 and 1 are a consequence of the high electron density in the practically unexpanded lychee core at t_acc = 20 fs; under these conditions, electrons are only formally assigned to ions according to their ion–electron cutoff distance (2 Å) rather than actually trapped. (Ions that formally acquire negative charges by close encounters with electrons are counted as neutral atoms; the resulting excess electrons are counted as free nanoplasma electrons.) Using the ion charge abundances at t_L = 3.5 ps, where the effective ion charges of 0 and 1 are missing, the lychee model gives nearly identical results, 〈E〉 = 1.92 keV and E_M = 19.6 keV, and a very similar KED (not shown in figure 8(a) for t = t_L), provided that the number of nanoplasma electrons is chosen such that Q(1061) is adjusted to the same value as that at t = t_acc. Note that Q(1061) at t = t_L is lower by almost 500 charge units than that at t = t_acc (figure 5) due to the rapid expansion of the cluster periphery into the surrounding dilute electron cloud. While the reentry of the electrons into the expanded cluster periphery has a relatively small effect on the ion energetics of the trajectory, the 500 charge units make a considerable difference to n_p^(free) as one of the key input parameters of the lychee model; the inclusion of the 500 electrons lowers 〈E〉 and E_M to 1.25 and 15.9 keV, respectively. In conclusion, the lychee model results for 〈E〉 and E_M are insensitive to the abundances of low ion charges, as long as the total cluster charge and the highest ion charge remain unaffected, reflecting the essence of equations (17) and (21).

The ion charges are not homogeneously distributed over the entire cluster as presumed by the lychee model. Rather, higher ion charges are enriched toward the cluster periphery. The presumption of a spatially homogeneous distribution of ion charges does not affect E_M but is expected to influence the shape of the KED. In particular, the larger number of higher charge states in the cluster periphery will enhance the high-energy tail of the KED.

Proceeding to case (B) (figure 8(b)), the deviations of the lychee model results manifest extreme non-vertical effects. With the abundances of reduced ion charges and the number of free nanoplasma electrons taken at t_acc = 200 fs, the lychee model underestimates 〈E〉 by a factor of 5. To reproduce the trajectory calculation result of 〈E〉 = 53 eV, the lychee model requires a central charge value Q(1061) of ≈420, which is approached only at t ≈ 4–5 ps (cf figure 6(a)), when 〈E〉 of the trajectory calculation also converges. In view of this slow charge build-up, it is therefore at first surprising that the lychee model overestimates E_M. Taking the maximum effective ion charge q_M^(eff) = 4 (figure 7(b)), the lychee model yields E_M = 0.44 keV as compared to E_M = 0.33 keV from the trajectory calculation. However, a closer inspection of the trajectory reveals that although the maximum effective ion charge of 4 is reached already at an early stage of the trajectory, it is not constantly attributed to the same ions. Short-time averaging shows that the maximum ion charge fluctuates between 2 and 4, reaching its final constant value of 4 only at t ≈ 0.5 ps. The shape of the KED calculated by the lychee model is dominated by the low-energy component h₀, which is by one order of magnitude higher than that of the trajectory calculation. In the trajectory calculation, the sharp h₀ peak is spread over a larger number of low-energy intervals, reflecting the large departure of the spatial charge distribution (figure 6(a)) from that of a lychee structure.

In case (C) (figure 8(c)), the E_M value calculated by the lychee model is 201 keV, being too high by almost a factor of 3, as a consequence of the neglect of the large cluster expansion until t_acc, R(t_acc) = 8R₀ (figure 1(c)). Furthermore, the lychee model predicts a large low-energy component h₀ of the KED, since the lychee model presumes all nanoplasma electrons to be located in the lychee core due to the neglect of electron dynamics. A neutral cluster core is not found in the trajectory calculation, as the electrons are nearly equally distributed over the entire ion framework by resonance heating.

For carrying out trajectory calculations, it is of practical interest to define some guidelines on what is the temporal length of a trajectory that is necessary for achieving convergence of the ion energies. From the foregoing discussion one can infer that a sufficient criterion for the convergence of the corrected ion energies is the convergence of the charge distribution inside the cluster and of the total cluster charge. This includes the absence of nuclear overrun subsequent to the truncation of the trajectory, as nuclear overrun can lead to an exchange of ion energies. However, as examples (A) and (B) have shown, the effect of the remaining variations of the charge distribution on ion energies fades out at a late stage of the cluster expansion, so that the criterion of a converged charge distribution is sufficient but not necessary. As an empirical rule of thumb, we found that the corrected ion energies are sufficiently converged when 〈T〉/〈E〉 ⩾ 0.85.