1. Introduction

Atomic clusters are aggregates of approximately 10² to 10⁷ atoms held together by van der Waals forces. Clusters are formed when gas is injected into vacuum at super-sonic velocity and cools. Clusters form when the induced interatomic/molecular dipole-dipole potential (van der Waals interaction) cannot be overcome by the thermal energy. The result is bonding of many atoms or molecules to form clusters. Each cluster has a density close to that of a solid, while the distance between two clusters is large giving the gas of clusters a mean mass density comparable to that of the original gas. When a strong laser pulse irradiates the clusters, its electric field couples with the clusters as a solid material, while it interacts with a large number of clusters because the average material density is low due to the space between clusters. The interaction uniquely combines features of laser-solid and laser gas interactions.

Laser-cluster interactions have several applications, for example, as a source of radiation such as EUV or X-rays that is applicable to lithography [1] or X-ray microscopy [2], and as a source of energetic electrons or ions [3–6]. The dynamics of the exploding clusters also gives rise to interesting nonlinear optical effects such as harmonic generation [7] and self-focusing [8].

In order to analyze the laser-cluster interaction by computer simulation, a range of physical processes must be modeled including ionization, collisions and charged particle motion [9–19]. Since clusters formed in a gas jet are concentrations of neutral atoms, the strong laser field first ionizes the atoms in the cluster, which is followed by heating and acceleration of the generated charged particles. In some cases the generation of the charged particles and their heating and acceleration take place simultaneously, so that the simulation must include all these effects. Even after atoms in a cluster are ionized, there is a period of time before the cluster expands. During this time the density is sufficiently high that collisions between charged particles are important.

In our previous work [16–18], the dynamical behavior of ions and electrons in an argon (Ar) cluster excited by a strong laser field was analyzed using a particle-in-cell (PIC) code. We found strong heating caused by a resonant interaction between the oscillating field of the laser and the transiting of electrons through the cluster. In the simulation, ions in the cluster were assumed to be in a fixed charge state and only the elastic scattering of electrons by ions was included. This assumption is valid when the laser field is high enough and ionization occurs early in the laser pulse before significant heating occurs. If the pulse duration of the laser is sufficiently short, ionization processes must be solved simultaneously with the plasma dynamics.

In order to apply our PIC code to a wider parameter range of cluster plasma, we have developed a new simulation code. Our new code includes the following processes: field ionization, impact ionization, electron-ion scattering and electron-electron collisions. Since the maximum electron velocity is not relativistic in the considered intensity range (I<1017 W/cm2), we treat these processes non-relativistically. Ion-ion collisions are neglected because the collisional time between ions is so long that it is not important during the ultra-short laser-plasma interaction.

Using the new code, we have analyzed two phenomena. One is the large energy absorption by resonant heating, which we described in the previous papers. By simulating the cluster dynamics from the neutral state, we confirm the existence of the resonant heating and estimate its threshold intensity more accurately.

The other is the angular dependence of the ion energy spectrum produced by cluster explosion. Recently, Symes et al. [20] and Chen et al. [21] found an anisotropy in the ion energy distributions. We have calculated the ion energy distribution using our new code and have compared the simulation results with their experimental data.

In this paper, Sec. 2 is devoted to explain the basic features of the PIC code and how collisional processes are implemented. The implementation of ionization processes in the PIC code is described in Sec. 3. Resonant laser heating accompanied by ionization processes is described in Sec. 4. The analysis of the anisotropic ion energy spectrum is described in Sec. 5. Conclusions and discussions are in Sec. 6.

2. Collision processes in a PIC code

Particle-in-cell (PIC) codes are usually used in the simulation of microscopic plasma physics. The major difference between a PIC code and other particle methods such as molecular dynamics is the method for calculating the force acting on each particle. In molecular dynamics simulations [9,10,18] all binary forces acting on all particle pairs are calculated directly. Hence the calculation of forces requires on the order of N ² numerical operations, where N is the number of particles. On the other hand, a PIC code introduces a spatial grid, on which electric field and charge density values are accumulated. The force on a particle located in a cell is determined by interpolating values of the electric field on neighboring grid points. Similarly, the charge density on each grid point is determined by the number and position of charges in neighboring cells. The electric field on the grid points is determined as the solution of Poisson’s equation. The required number of operations for the calculation of the density and for the interpolations of field values at the particle positions are both proportional to N. Computation of the field on the grid is independent of the particle number, and does not significantly increase the total number of operations as long as the total number of particles is much larger than the total number of grids. As a consequence, the total number of numerical operations to advance one time step in a PIC code is proportional to N.

The efficiency of the force calculations is one of the great advantages of a PIC code, and hence PIC codes are widely used for modeling plasmas. However, the use of a grid alters the description of binary processes, i.e. collisions. One can think of the force of interaction between two charges as being filtered by the grid. This can be an advantage for PIC codes in situations where collisions are unimportant in the physical system under study. The PIC method allows for the simulation of such cases with fewer simulation particles than real particles. A gridless simulation accurately accounting for binary forces would require a number of simulation particles equivalent to the actual number of particles, so as to correctly model collisions. Using a smaller number of particles, but each with a bigger weight results in an overestimate of the effect of collisions. Further the integration of the trajectories for two particles passing close to each other requires a much smaller time step if the full interparticle force is calculated.

When the plasma density is as high as solid density, binary collisions become important in determining physical processes such as heat conduction and energy relaxation. These processes are not treated in a conventional PIC code, and must be restored by modifying the PIC algorithm. Takizuka and Abe (TA) proposed one such method [22], which has been further refined by Ruhl [23] and Mikhailova [24]. They treat collisional processes as random momentum exchanges, which are assumed to occur many times during one time step of the simulation. In order to calculate the processes efficiently, they divide the whole simulation box into small cells (which may encompass several PIC cells to guarantee enough particles to give good statistics) pick pairs of particles in each cell, and give random momentum increments to each pair, based on the Coulomb collision cross-section. When the size of a cell is small enough, the number of particles in a cell is also small, so that the calculation of pairs of collisions is not as large as N ². However, a sufficiently large number of pairs must be selected so as to achieve good statistics.

The method we propose is based on the Fokker-Planck description of Coulomb collisions [25]. In our implementation each particle will receive a random change in momentum during each time step. The properties of the momentum change will be selected based on the average plasma properties in a cell so as to mimic the effects of binary collisions. Further, the momentum changes occurring in an individual cell will be adjusted so as to conserve momentum and energy in each cell. Effectively, our calculation models collisions with all particles in a given cell

The Fokker-Planck equation for the evolution of the distribution function f(r,v,t) of electrons in the presence of random velocity increments is given by

Here the velocity drift and diffusion coefficients are given by

In our previous work we treated only electron-ion scattering. When electrons scatter from ions the direction of their momentum changes, but not its magnitude. Within the Fokker-Planck treatment, this is achieved by specifying $\Delta \text{v}$ such that D and F have the form

To treat electron-electron collisions, we introduce the model collision operator

We note that Eq. (4) is just a model collision operator and has a number of features that need to be explained and modified before implementation in our code. First, Eq. (4) differs from the collision operator for scattering of test electrons by a Maxwellian distribution of field electrons [26] in that the diffusion in velocity is isotropic. For the case of test electrons scattering from field electrons [26] the scattering in velocity space is anisotropic; pitch angle scattering occurs at a different rate from scattering in magnitude of velocity. The rates are similar in order of magnitude, but not equal. Since electron pitch angle scattering is dominated by scattering from ions with Z > 1, we have simplified the situation by making the electron pitch angle scattering rate equal to the scattering rate for the magnitude of the velocity. Second, Eq. (4) describes collisions of test electrons with a zero mean velocity, thermal ensemble of scatterers. Thus, to implement it in the code we compute in each cell the mean electron velocity and the mean kinetic energy of electrons in a frame moving with the mean velocity. Velocity changes are then computed in the moving frame based on a temperature equal to 2/3 of the mean kinetic energy in that frame. Finally, Eq. (4) does not conserve momentum or energy as electron-electron collisions must do. To correct this shortcoming we add velocity increments in addition to the random ones that guarantee conservation of momentum and energy.

The final prescription for the random velocity increment of the i-th electron takes the form

To demonstrate how our numerical procedure works, we calculate the temporal evolution of an ensemble of electrons, whose distribution function is initially non-Maxwellian. Our simulation consists of a spatially uniform plasma consisting of 32² cells, each containing approximately 160 simulation electrons, and the electron density is 1.742 × 10²¹ cm⁻³. The black curve in the Fig. 1 represents the initial distribution as a function of energy. It is non-Maxwellian with average electron energy of 15eV. As time passes, the distribution gradually approaches to a Maxwellian. The electron-electron collision time for these parameters is 3.4 fs. The final light-colored curve well matches to a Maxwellian with temperature of 10 eV, which is 2/3 of the initial average energy. The result indicates that our scheme provides the electron-electron collisions appropriately to a PIC platform.

 

Fig. 1 Temporal evolution of electron energy distribution. The black curve is the initial non-Maxwellian distribution and the light-colored curve is the asymptotic distribution simulated by the code. The dashed curve is the Maxwellian distribution.

Download Full Size | PDF

3. Ionization processes

Our code includes both field ionization and impact ionization. Since field ionization is determined by the strength of the electric field at the position of each atom or ion, the number of operations scales linearly with N. Field ionization is implemented as a stochastic process. First, the ionization rate appropriate to the instantaneous electric field for each atom is calculated using the ADK formula [27]:

Here, U_i is the binding energy of the i^th electron measured in $U_0=e^2{(4\pi {\epsilon }_0{a}_B)}^{-1}$ (one hartree), E is the magnitude of the electric field strength normalized by $E_0=e{(4\pi {\epsilon }_0{a}_B^2)}^{-1}$, where a_B is the Bohr radius, and ω_A is the atomic frequency, ω_A = $U_0/\hslash$. $C_{n^{*}{l}^{*}}^2$ and f (l, m) are defined as follows:

After the rate of ionization W(E) for each atom is calculated, one random number a (uniformly distributed, 0<a<1) for each atom is generated and it is compared with 
R = 1-exp[-W(E)Δt], where Δt is a time step for the calculation of the ionization process. When a<R, a new electron is added and the charge state Z of the atom is increased by one.

Field ionization occurs for each atom independent of other charged particles. Impact ionization is caused by an interaction between an atom and electrons, similar to a collision. This means that direct calculation of impact ionization would require a number of operations proportional to the product of the number of electrons and the number of ions, N_iN_e. In order to avoid this situation, we treat impact ionization as follows. According to Lotz [28], the cross section for impact ionization by an electron with normalized kinetic energy U is represented by the following empirical formula:

Here, U_i is the lowest binding energy of electrons remaining in the ion, n_i is the number of equivalent electrons with binding energy U_i, and a_i, b_i and c_i are individual constants based on in Ref [28]. Using the cross section, we first prepare the impact ionization rate,

In the case of impact ionization, the new electrons are created with a kinetic energy $U^{\prime }$. Since $\Delta t{W}_i(U_i)$is the probability for producing a free electron with any energy $U^{\prime }>0$ in the time interval $\Delta t$, we consider the probability of generating and electron with energy $U^{\prime }>\Delta U$ in the same time interval to be $\Delta t{W}_i(U_i+\Delta U)$. The impact must be such so as to raise the target electron’s energy by an amount $\Delta U$ in excess of the binding energy. We can then use the definition of conditional probability ($P(U^{\prime }>\Delta U|U^{\prime }>0)P(U^{\prime }>0)=P(U^{\prime }>\Delta U\cap U^{\prime }>0)=P(U^{\prime }>\Delta U)$) to find the probability, given an ionization, that the new free electron has energy $U^{\prime }>\Delta U$,

The probability density distribution for creating an electron with energy $U^{\prime }$in the range $d(\Delta U)$centered at $\Delta U$is thus given by

In the code, rather than compute this probability distribution function and initialize new electrons with a random energy consistent with it, we initialize new electrons with a Maxwellian distribution of velocities with a temperature determined by the expected value of $\Delta U$. Specifically, $T_{new}=(2/3)〈\Delta U〉$, where

To keep electron energy in balance we tabulate the energy required to make the new electrons, $U_{new}={\displaystyle \sum _{new}\left(U^{\prime }+U_i\right)}$. We then decrease the energy of each of the impacting electrons by multiplication by a factor

The numerical procedures of the impact ionization explained in this section have been demonstrated for uniform plasmas and the results are in a good agreement with the analytical prediction.

4. Resonant heating of an argon cluster

In previous papers [16–18] we reported the dynamical properties of an argon cluster irradiated by a strong laser as calculated by a PIC code. The results show that strong heating takes place when the laser intensity exceeds a threshold value, which is proportional to the square of cluster size. If the atoms in the cluster are ionized, the laser field is shielded from the interior of the cluster by the response of free electrons. The strong laser field can extract electrons from the cluster in the first half of a laser cycle and accelerate them back to the cluster in the second half of the cycle [29]. If the electrons gain sufficient energy, such that they transit through the cluster and emerge on the other side in the accelerating phase of the electric field they will continue to gain energy [16–18]. We note that a number of different resonance heating mechanisms [16–19,30–33] have been proposed. The mechanism we describe, in which the electron oscillation period depends on the laser intensity, is distinct from the originally proposed mechanism in which resonance occurs when the laser frequency matches the linear electrostatic mode frequency of the expanding plasma sphere [30,31]. It is related to, but distinct from resonance heating in which the period of an energetic electron in the confining electrostatic potential matches the laser frequency [32]. A partial review of these mechanisms can be found in Ref [34].

In order to evaluate the threshold intensity and to study the mechanism of resonant heating, a large number of simulation runs were done for a wide range of laser intensities and cluster sizes. However, our previous code assumed that the charge state of argon atoms is fixed initially and the ionization process was not studied. The heating was found to be relatively insensitive to the initial charge state. In this section, we describe the resonant heating simulated by our new code including collisional and ionization processes.

Figure 2 shows the temporal evolution of the average charge state of argon ions for different laser intensities. In these runs, the laser field is imposed on the cluster as a uniform oscillating electric field whose frequency corresponds to a laser with a wavelength of 800 nm. The pulse shape of the laser is a Gaussian, which is shown by the dashed line in the figure, and its pulse duration is set to be 100 fs (FWHM). The cluster consists initially of neutral argon atoms of density is 1.7x10²² cm⁻³. The simulation is two dimensional (x-y), and the cluster is a cylinder with the diameter of 38 nm. These conditions are the same as in our previous papers.

 

Fig. 2 Temporal evolution of the average charge state of argon ions in a cluster for different laser intensities.

Download Full Size | PDF

The atoms are neutral until the laser intensity reaches the threshold value for field ionization. Thus, ionization appears earlier for higher laser intensity. Once ionization starts, electrons are heated by the laser field and contribute to impact ionization. As intensity is increased from I = 1x10¹⁵ W/cm² to I = 1x10¹⁶ W/cm² there is a large increase in the average charge state due to the resonant heating mechanism mentioned above.

In order to evaluate the threshold intensity more clearly, the absorbed energy as a function of the laser intensity is shown in Fig. 3 for this set of cluster parameters as well as others to be discussed later. The figure shows the total energy in electrons, ions and electric field at the time of 306 fs, divided by the number of electrons. Figure 3 shows that the absorbed energy abruptly increases at the intensity of 7x10¹⁵ W/cm². The results are almost the same as those obtained by the calculation with the fixed charge state Z = 8 in the previous papers. The main difference from the previous results is the threshold intensity, I = 7x10¹⁵ W/cm², is a little higher than the result of the fixed charge state, I = 5x10¹⁵ W/cm². The difference of the threshold intensity mainly results from collisional effects. The stopping range of electrons inside the cluster is shorter than that calculated by the previous code, and hence a higher intensity is necessary to satisfy the condition that the electron transit time back and forth through the cluster match the laser period. On the other hand, we note that the amount of heating is not affected by the variation of the charge state. Additionally, we have examined the time dependent cluster polarizability that is responsible for the self-focusing observed in a gas of clusters [8]. This too is essentially unchanged from that observed in our previous simulations with fixed charge state.

 

Fig. 3 Average absorbed energy per electron versus laser intensity for several simulations.

Download Full Size | PDF

5. Anisotropic expansion of a cluster plasma

Generation of energetic ions is one important application of laser-cluster interactions. In our intensity regime, fast ions are accelerated by the electrostatic field produced by the charge separation created by energetic electrons that attempt to leave the cluster. Accordingly, ion acceleration depends upon the distribution of the expanding electrons, and the strong heating of electrons by the laser should affect ion acceleration.

Recently, Symes et al. [20] and Chen et al. [21] measured the angular dependence of the energy distribution of ions in laser-cluster interaction experiments using time-of-flight (TOF) detectors. Ion spectra were measured for the laser polarization rotated along and perpendicular to the ion TOF tube axis. The results show, in some cases, anisotropic energy distributions that depend upon the size of clusters. In this section, we describe a simulation study of ion expansion using our new code.

For these simulations we first consider hydrogen clusters, in which case any resulting explosion anisotropy must not be due to the anisotropic distribution of the ion charge state [20]. We choose parameters corresponding to the experiments of Chen et al [21]. Specifically, we consider irradiation of 22 nm clusters (ionic density $4.1\times {10}^{22}$ cm⁻³) with various intensity laser pulses of 60 fs FWHM duration. Our simulations in this case are both 2D, with reflection symmetry in the direction transverse to the laser electric field, and 3D, with reflection symmetry in the two directions transverse to the laser electric field . For these simulations we use open boundary conditions, meaning the Poisson equation is solved in an infinite domain assuming no charge outside the simulation volume, and particles that leave the simulation domain only interact with the laser field, and their positions are not tracked. This is intended to model the case of very low cluster density in which little inter-cluster interaction and heating can be expected.

The first thing we note in Fig. 3 is that the energy absorbed as a function of intensity curve for the 3D case is more gradual and with a lower threshold for strong heating than in the 2D case. The lower threshold for 3D simulations was observed in our early work [16,17] on argon clusters and was attributed to the fact that for a given intensity the electric field at the surface of a spherical cluster is greater than that for a cylindrical one. The magnitude of the electric field at the surface determines the threshold because it determines the energy of electrons passing through the cluster. The more gradual transition in 3D is at this time unexplained. We note that in our previous 3D simulations of 38 nm Argon clusters [16] without electron-electron collisions and ionization, a clear threshold was observed. In the new simulations there appears to be more heating at low intensities in the 3D simulations than in the 2D simulations. This extra heating obscures the transition and its cause is under investigation. The ion spectra at the end of the simulations strongly depend on the laser intensity.

This is illustrated in Fig. 4 , where spectra for three different intensities, (2 × 10¹⁴ W/cm², 7 × 10¹⁴ W/cm², and 2.5 × 10¹⁵ W/cm²) at 308 fs (well after the laser pulse, FWHM = 60 fs) from 3D simulations are displayed. Two curves are displayed for each intensity showing the spectra for ions with a range of angles (π /16) close to being parallel (solid symbols) or perpendicular (open symbols) to the laser field. For the lowest intensity case (2 × 10¹⁴ W/cm²) the spectra decrease monotonically with energy, out to a maximum energy of about 2.5 KeV. There is anisotropy in this case with more ions emerging in the parallel than in the perpendicular direction. As the intensity is increased the spectra become flat and then inverted, and extend to higher energy (~10 keV). In the high intensity cases there are beam like features at high energy that have been noted previously [17].

 

Fig. 4 Ion energy spectra for a 22 nm hydrogen cluster and intensities I = 2x10¹⁴ W/cm², I = 7x10¹⁴ W/cm² and, I = 2.5x10¹⁵ W/cm², as labeled, and a weighted average of spectra for intensities 0.1-3x10¹⁵ W/cm². In each case, the solid symbols show the distribution parallel to the polarization, while the open symbols show the distribution perpendicular to the polarization.

Download Full Size | PDF

The ion spectra in the experiment are the accumulated results of the explosion of a number of clusters, which have sampled different intensities due their different locations in the plane transverse to the axis of the laser pulse. We can simulate this effect by forming a weighted average of the spectra calculated for different intensities. This weighted average is also shown in Fig. 4. To form the average we weight spectra computed for intensities between 1 × 10¹⁴ W/cm² and 3 × 10¹⁵ W/cm² for the data points indicated on Fig. 3. The weighting is inversely proportional to intensity and proportional to the interval in intensity surrounding a data point. This corresponds to equal area weighting assuming that the transverse profile of the laser intensity is Gaussian. It is clear that the spectra produced by this weighting process are substantially different than the spectra produced by single intensity simulations. Thus, any quantitative attempt to understand the spectra in experiments must include the effect of intensity variations. Since the threshold for strong heating is cluster size dependent ($I\propto d^2$ [16,17], ) variations in cluster size will play a role similar to variations in intensity in determining the spectra.

The weighted spectra produced by the 3D simulations can be compared with those in Ref [20]. The simulated and measured spectra have the following features in common: spectra are monotonically decreasing with energy, the maximum energy is in the range 10 – 20 keV, there is an anisotropy with the number of high energy ions in the parallel direction exceeding that in the perpendicular direction, and the maximum energy in the parallel direction is only slightly higher than the maximum in the perpendicular direction. To examine the origin of the anisotropy in the simulations we plot the laser cycle averaged radial electric field in the parallel and perpendicular directions at different times during the simulation of the cluster explosion in the 1 × 10¹⁵ W/cm² case. Specifically, Fig. 5 shows $〈E_x(x,y=0)〉$ and $〈E_y(x=0,y)〉$ where the laser electric field is directed parallel to x, at t = 90, 125, and 143 fs (The laser pulse peaks at 104 fs with a FWHM of 60 fs.) and the angular brackets indicate a time average over one period of the laser field. Since the plotted field is averages over a laser cycle, what is plotted reflects the charge density surrounding the cluster. We note that early in the pulse (t = 90 fs) the core of the cluster is still intact, and the electric field is small for r < 10 nm due to shielding by the dense core. The profiles of the radial electric field in the two directions are similar, but with a significant difference. Specifically, the electric field parallel to the laser electric field is larger at smaller radius and where there are ions to be accelerated in qualitative agreement with the model of Ref [33]. The overlap of the field and the ions causes the number of ions accelerated in the parallel direction to be greater than the number in the perpendicular direction. The phase space of electrons in the (x, v_x) plane at y = z = 0 is shown in the inset of Fig. 5 for t = 90 fs. It shows the presence of energetic electrons that have been pulled from the cluster core and then reaccelerated back into it [16,17]. These electrons have sufficient energy to travel through the cluster and interact with the field on the other side creating a halo of dynamic electrons surrounding the cluster. Thus, in this regime the electric field cannot be calculated assuming the electrons are a cold fluid in force balance as considered in Ref [33].

 

Fig. 5 Radial profiles of radial electric fields as a function of distance from the center of the cluster at several different times during the cluster explosion for I = 1x10¹⁵ W/cm². For the solid (dashed) curves distance is measured along the direction parallel (perpendicular) to the laser polarization.

Download Full Size | PDF

Profiles of the radial electric field at later times are also shown in Fig. 5. By t = 125 fs the cluster has expanded to the point where the core is no longer shielded. The parallel and perpendicular radial electric field profiles are still quite similar at this time indicating that the electron cloud surrounding the cluster is isotropic. However, by 143 fs there is a clear anisotropy in the electric field profiles. Figure 6 shows the cycle averaged ion density profiles at the times for which the electric field profiles are shown in Fig. 5. The solid curves correspond to the density along the x-axis, i.e. parallel t the laser field, and the dashed curve correspond to the density along the y-axis, i.e. perpendicular to the laser field. At t = 90 fs the ion density profiles in the two directions are similar indicating the core of the cluster is still spherical. However, the ion density profile in the parallel direction has developed a pedestal extending from 13 to 19 nm. This pedestal represents ions that are beginning to be accelerated. A similar, but smaller pedestal appears in the profile corresponding to the perpendicular direction. At t = 125 fs the ion density profile has expanded significantly in the parallel direction and appears to have contracted in the perpendicular direction. Actually, this apparent contraction is the result of ion motion away from the plane x = 0, the ion density being ellipsoidal at this time. In summary, the observed anisotropy in the ion spectra in these simulations is a consequence not so much of large differences in the electric fields in the parallel and perpendicular directions, but rather a consequence of differences in the number of ions experiencing the average accelerating electric fields. More ions find themselves in the accelerating field in the parallel direction than in the perpendicular direction. A consequence of this is that the maximum energy attained by ions accelerated in the two directions is not as different as the number of ions in each direction accelerated to high energy.

 

Fig. 6 Radial profiles of ion density (here n₀ = 4.1 × 10²² cm⁻³) as a function of distance from the center of the cluster at several different times during the cluster explosion for I = 1x10¹⁵ W/cm². For the solid (dashed) curves distance is measured along the direction parallel (perpendicular) to the laser polarization. In each case, the solid symbols show the distribution parallel to the polarization, while the open symbols show the distribution perpendicular to the polarization.

Download Full Size | PDF

A different picture emerges when we consider larger argon clusters. Figure 7 shows ion energy distributions for exploded argon clusters for three different peak intensities, 5x10¹⁵ W/cm², 1.5x10¹⁶ W/cm² and 3x10¹⁶ W/cm². The pulse duration of the laser is again set to be 60 fs FWHM, the diameter of the cluster is 35 nm and the ionic density is $1.7\times {10}^{22}$cm⁻³. These simulations are two dimensional because of the added computational burden associated with treating collisions and ionization in a higher Z material. In the argon runs, open boundary conditions were used to calculate the electric field, and particles that leave the simulation domain are tracked in momentum, but not position. The curves in Fig. 7 are snapshots of the distribution functions at the time 440 fs, at which point the laser pulse has already terminated. In the figure, the solid curves indicate the energy distributions of ions observed in the direction parallel to the laser polarization (θ₀ = 0), and the dashed curves indicate them observed in the direction perpendicular to the polarization (θ₀ = π/2). The distributions in each angle, θ₀, are again evaluated by counting the number of ions within 
π /16 of θ₀.

 

Fig. 7 Ion energy spectra for 35 nm argon clusters and intensities I = 5x10¹⁵ W/cm², I = 1.5x10¹⁶ W/cm² and I = 3x10¹⁶ W/cm². In the graph, the solid curves show the distributions parallel to the polarization, while the dashed curves show the distributions perpendicular to the polarization.

Download Full Size | PDF

As shown in Fig. 7, when the laser intensity is low, the ion energy distribution is almost isotropic. The distributions parallel and perpendicular to the laser field are equal, except for the maximum energy, which is higher in the parallel direction. This is in contrast to the hydrogen case where anisotropy was seen even at lower intensity, and low energy. The difference in isotropy is due to the increased importance of electron collisions in the argon case. This is primarily due to the higher charge state and higher electron density of the argon clusters. As a result the electron distribution is more isotropic, and ions are mainly accelerated in the radial direction, and the ion distribution function becomes isotropic as well. Although for the most part of the parallel and the perpendicular distributions in I = 5x10¹⁵ W/cm² are the same, we note the maximum energy of ions in the direction parallel to the laser field is slightly larger than that perpendicular to the field. The most energetic ions tend to be the ones that were first accelerated, and the presence of these very high-energy ions indicates that at early times the acceleration process was not isotropic. When the laser intensity is high, it is found that the ion energy distribution becomes slightly anisotropic in the sense that at moderate energies there are more ions in the parallel direction than in the perpendicular direction. However, the maximum energies in the parallel and perpendicular directions are now more equal than in the low intensity case. This change correlates with the intensity threshold discussed earlier. Figure 3 shows the threshold intensity for resonant heating in this case is estimated to be about 8 × 10¹⁵ W/cm², which is intermediate to the intensities 5x10¹⁵ W/cm² and 1.5x10¹⁶ W/cm². Thus, the anisotropy observed in the higher intensity cases occurs for intensities above the threshold, where electron energies are higher and the effects of collisions are reduced.

A simple estimate of the collision rates in the two cases confirms the difference between the argon and hydrogen cases. We note that the electron-electron scattering rate for electrons with speed v is given by Eq. (5). The scattering rate for electrons colliding with ions is similar, but multiplied by the charge state Z. The importance of collisions for intensities at the threshold can be estimated as follows. For an electron to be resonant its speed must be high enough to carry it through the cluster in half a laser period, $\text{v}=2cd/\lambda$, where d is the cluster diameter and λ the laser wavelength. Collisions will be important if the collision rate approaches the transit time of electrons through the cluster ${\nu }_{s}d/\text{v=}{\nu }_s\lambda /\text{(2}c\text{)}\approx 1$. For the parameters of our simulations we find ${\nu }_{s}d/\text{v}=0.4$ for the argon case, and ${\nu }_{s}d/\text{v}=6.0\times {10}^{-3}$ for the hydrogen case. Thus, at threshold, the electron distribution will be isotropic due to collisions in the argon case, but not in the hydrogen case. This explains the qualitative difference between our hydrogen and argon simulations. We note the above arguments imply an interesting scaling with cluster size. Normally, one would expect collisions to become more important as the diameter of a cluster is increased. This will be the case if intensity is held fixed. However, it requires a higher intensity to reach the resonance threshold for a larger cluster. The speed of electrons at resonance is twice the diameter divided by the laser period. Because of this higher speed, collisions are less important at threshold for larger clusters. This has been verified by comparison of 2D simulation results for large (22 nm) and small (5 nm) hydrogen clusters. A second issue to consider is that we have been considering short pulses (< 100fs), and the threshold for resonant heating can be lower for longer pulses [18]. This is due to the fact that resonance can occur when the transit time is multiples of the laser period, which occurs for lower electron energies and correspondingly lower intensities. However, these resonances are only observed in simulations when the pulse duration is increased to 400 fs and beyond. The role of collisions on cluster explosions with longer pulses remains to be explored.

6. Conclusion and discussion

We have developed a new particle-in-cell code, which includes collision and ionization processes in order to analyze laser-cluster interaction. The schemes used for the implementation of these processes are selected to keep the number of numerical operations proportional to the number of particles. In order to calculate the effect of electron-electron collisions, we have developed a method based on the Langevin equation, which not only maintains numerical efficiency but also conserves the total energy and momentum of electrons. The ionization processes included in the new code consist of field ionization and impact ionization. Both of these are calculated via Monte Carlo techniques.

We have applied our new code to investigate the resonant heating that was found in our previous work, and that was analyzed with a PIC code assuming fixed charge-state ions. The new results show that the resonant heating also occurs when the dynamics of ionization are included and the averaged charge state of ions increases during the pulse. The difference between the result of the previous code and that of the new one is that collisions can raise slightly the threshold intensity. This effect will be dependent on the cluster material.

As one of the applications of our new code, we have investigated the anisotropic ion acceleration observed in experiments [20,21]. The results show that the occurrence of anisotropy depends upon the intensity of the incident laser, the size of the cluster, and the composition of the cluster. In the case of size and constitution of the clusters. High Z clusters are more collisional and tend to exhibit reduced anisotropy. In this case when intensity is raised the anisotropy increases due to the decreased importance of collisions. However, the anisotropy is less in the argon case than in the hydrogen case.

We expect our new code to be applicable to the wide range of high-density plasma, and strong laser pulse interactions, especially laser-solid or laser-liquid interactions. The applicable intensity range of the current code is limited to the non-relativistic regime (I<10¹⁷ W/cm²) since it is electrostatic and collision rates are calculated based on nonrelativistic formulas. However, since the probability of collisions becomes small in the relativistic regime, our model can be applied to problems with ultra high intensities. However, at these intensities the magnetic fields and the Lorenz force must be taken into account.

One of the advantages of using PIC codes is their ability to simulate large numbers of particles. Recently, simulations using molecular dynamics in which all particle interactions are binary have been considered for large clusters [34]. Comparison of these results with PIC results is a topic for future investigation.