Fully microscopic analysis of laser-driven finite plasmas using the example of clusters

We discuss a microscopic particle-in-cell (MicPIC) approach that allows bridging of the microscopic and macroscopic realms of laser-driven plasma physics. The simultaneous resolution of collisions and electromagnetic field propagation in MicPIC enables the investigation of processes that have been inaccessible to rigorous numerical scrutiny so far. This is illustrated by the two main findings of our analysis of pre-ionized, resonantly laser-driven clusters, which can be realized experimentally in pump–probe experiments. In the linear response regime, MicPIC data are used to extract the individual microscopic contributions to the dielectric cluster response function, such as surface and bulk collision frequencies. We demonstrate that the competition between surface collisions and radiation damping is responsible for the maximum in the size-dependent lifetime of the Mie surface plasmon. The capacity to determine the microscopic underpinning of optical material parameters opens new avenues for modeling nano-plasmonics and nano-photonics systems. In the non-perturbative regime, we analyze the formation and evolution of recollision-induced plasma waves in laser-driven clusters. The resulting dynamics of the electron density and local field hot spots opens a new research direction for the field of attosecond science.


4
The purpose of this paper is to test the validity and applicability of MicPIC and to identify some new physical processes that have become accessible through MicPIC. In the first part of the paper, the code is applied to investigate the regimes of linear and nonlinear laser cluster interactions by using pre-ionized clusters in resonance with few-cycle infrared laser fields. We consider singly ionized spherical nanoplasmas with a density chosen such that the Mie plasmon is in resonance with the driving 800 nm laser field. In a real pump-probe experiment performed at optimal delay for resonant Mie plasmon excitation, this simplified scenario corresponds to the situation prior to the arrival of the probe pulse. Although the pivotal role of resonant collective electron excitation is well accepted for laser-cluster interactions [39][40][41], little is known about how the collective enhancement proceeds in large clusters for which electromagnetic propagation effects cannot be neglected. In the second part of the paper, we model a smaller Ar cluster starting from its neutral ground state and demonstrate the capacity of MicPIC to describe the full evolution process of laser-cluster interaction including ionization. The neutral cluster goes over into a warm, highly overdense and strongly coupled plasma that expands, driven by Coulomb and hydrodynamic forces. The results are found to be in very good agreement with MD simulations. The following conclusions can be drawn from our analysis.
First, considering resonant nanoplasma excitations in the linear response regime, we extend a previous systematic analysis of the radius-dependent absorption and scattering of [34] up to larger cluster sizes. The comparison of MicPIC and MD simulation results shows that electrostatic MD loses its validity for cluster radii of R 20 nm, where the propagation regime is entered and effects such as radiation damping, field attenuation (skin effect) and polaritonic plasmon red-shifts begin to substantially influence the optical response. Peak absorption crosssections of up to five times the geometric cross-section indicate the occurrence of a pronounced absorption paradox in resonant clusters in the size regime of some tens of nanometers. By fitting the MicPIC and MD data to the corresponding continuum approaches, i.e. to Mie theory and the so-called nanoplasma model, we extract the surface contribution to the collision frequency of the nanoplasma. The MicPIC analysis enables us for the first time to determine quantitatively the contributions of bulk and surface collisions to the collision frequency as well as their detailed impact on the electromagnetic nanoplasma response. It is shown that the contribution of surface collisions remains profound up to the propagation regime. Based on our results, the maximal plasmon lifetime can be related quantitatively to the competition between surface collisions and radiation damping. It should be noted that the capacity to calculate the collision frequency has important implications in the fields of nanophotonics and nanoplasmonics.
Secondly, on resonant nanoplasmas excited at higher laser intensity, we provide a detailed analysis of the generation and evolution of recollision-induced nonlinear plasma waves. The existence of such waves has recently been reported for noble gas clusters [34]. Here we perform a more detailed study of the formation and dynamical evolution of these plasma waves in a simpler system, namely a hydrogen-like cluster. The MicPIC analysis reveals that the plasma waves are created by electrons that are first expelled from the cluster, then accelerated outside the cluster and are finally driven back and recollide with the cluster surface. Because of the high nanoplasmonic field enhancement of the resonant cluster, this recollision process is most pronounced at the two cluster poles, creating plasma waves that propagate from there into the cluster. Due to the spherical cluster geometry, waves are focused toward the cluster center, where they eventually collide and break. The plasma wave dynamics results in strong density fluctuations and electric field hot spots on the nanometer space and attosecond time scales. The intensity of the transient field hot spots exceeds the laser field by more than two orders of magnitude, demonstrating the feasibility of strong nanoplasmonic field enhancement via nonlinear plasma wave breaking. Compared to previous results on resonant Xe clusters [34], the absence of atomic ionization leads to much clearer signatures of the plasma waves and substantially higher field enhancement. Because of the microscopic resolution of the plasma processes, MicPIC opens a route toward the fully ab initio description of turbulent plasma wave phenomena in classical plasmas. Moreover, as the electron density oscillations are occurring on an attosecond time scale, nonlinear waves in nanoplasmonic systems open a new area for the application of attosecond metrology.
Thirdly, the modeling of a small Ar cluster (R = 5 nm), starting from its neutral ground state configuration, demonstrates the capability of MicPIC to model overdense, strongly coupled plasmas including atomic-scale ionization processes. It shows the applicability of the method to typical dynamical aspects of laser-cluster interactions, such as microfield effects in the ionization dynamics, hydrodynamic nanoplasma expansion and Coulomb explosion.
The paper is structured as follows. Section 1 provides a brief description of the MicPIC and MD methods and the applied continuum models. In section 2, we discuss the radius-dependent absorption, scattering and Mie plasmon lifetime of the clusters in the linear regime, the formation and propagation of plasma waves for strong nonlinear excitation and the ionization dynamics in an initially neutral rare-gas cluster. Conclusions are drawn in section 3.

Methods
We begin with a brief description of the key idea and the implementation of the MicPIC method; for further details see [34].

The microscopic particle-in-cell approach
In MicPIC, each plasma particle (electron or ion) is described by a charge density ρ i (r) = q i g(|r − r i |, w 0 ), where q i and r i are the charge and position of the ith particle and g(r, w) = exp(−r 2 /w 2 )/π 3/2 w 3 is a normalized Gaussian-shaped function. The shielding parameter w 0 softens the Coulomb forces on a length scale of the order of one atomic unit and has to be used to ensure numerical stability. It emulates the effective shielding of Coulomb singularities by quantum uncertainty and by the finite width of particle wavefunctions. For example, it prevents electrons recombining classically below the quantum energy levels. This is a generally accepted procedure that is also used in MD simulations. The effective width w 0 can be adjusted, so that the minimum potential energy for an electron-ion pair is equal to a specified binding energy.
The classically exact dynamics of the ith plasma particle is driven by the force obtained by averaging the electric (E) and magnetic (B) fields over the particle charge density. The self-consistent evolution of the fields is given by Maxwell's equations as where the current density j = iṙ i ρ i depends on the velocitiesṙ i and effective charge densities ρ i (r) of plasma particles.
Departing from the above exact description, we introduce a PIC approximation by representing each particle by a wider Gaussian charge density ρ pic i (r) = q i g(|r − r i |, w pic ) with w pic w 0 . The larger PIC particle width is key to an efficient solution of Maxwell's equations, allowing the use of a coarse numerical grid and large time steps. Note that PIC underestimates and softens the fields of charged particles close to their origin. Therefore, PIC electric and magnetic fields are smoother than the actual fields. They evolve as The PIC fields are driven by the current density j pic = iṙ i ρ pic i . At the PIC level, radiation fields are fully accounted for; however, the microscopic nature of the particles is lost due to the large PIC particle size. The PIC force on the ith particle is To identify and approximate the missing short-range forces, the actual force on particle i can be formally split into a microscopic part f mic i and a long-range PIC part f pic i by subtracting/adding the PIC forces from/to the full expression in (1). This splitting has motivated the acronym MicPIC and yields This expression is identical to the force in (1). The short-range character of the microscopic contribution f mic i becomes evident after decomposing PIC and actual electric and magnetic fields into their individual particle contributions The above sum describes the force on the ith particle, created by the fields of all other ( j) particles. This decomposition is justified because of the linearity of Maxwell's equations. Note that fields and charge densities with superscript 'pic' and without superscript refer to the PIC and actual fields and charge densities, respectively. For all particles j far enough away from the ith particle (r i j w pic ), the actual and PIC fields produced in the region r ≈ r i are identical. Further, the variation of the (actual and PIC) fields over the PIC particle extent can be approximated by a linear Taylor expansion around r i . Due to the even symmetry of the charge density, integration over the linear field terms gives zero. As a result only the constant field terms remain and can be pulled out of the integral. The remaining integral over the actual and PIC particle density gives zero for each index j, as their total charge is equal. This proves the short-range nature of the microscopic correction. It should be stressed that the dynamics of the plasma particles is described exactly and independent of the width of the particles on the PIC level. The value of the PIC particle width only determines the softness of the force on the PIC level and, in turn, the radius within which the microscopic forces contribute.
Up to this point, everything has been derived in full generality. Now, as the only formal approximation in MicPIC, field retardation is neglected locally, i.e. within the microscopic correction, by taking the non-relativistic, electrostatic limit of (8). This greatly facilitates its numerical evaluation. Dropping magnetic fields and expressing electric fields by the respective Coulomb interaction yields For Gaussian-shape functions, the double integral can be evaluated analytically and yields the difference of the particle interaction energies for actual and PIC particles. The interaction energy of two Gaussian particles with width parameter w reads and allows rewriting (9) as Here r i j = |r j − r i | is the inter-particle distance. By combining the electrostatically approximated microscopic correction in (11) and the PIC force in (6) we arrive at the final MicPIC force The MicPIC dynamics is now determined by the self-consistent integration of Newton's equations of motion with the force specified in (13), together with (4) and (5).

Numerical implementation
For the efficient numerical solution of this coupled set of equations, we apply the P 3 M scheme, leading to linear scaling of the numerical effort with particle number N . For the particle-mesh part, Maxwell's equations for the electric and magnetic fields on the PIC level are calculated using the finite-difference time-domain method (FDTD) on a coarse equidistant three-dimensional grid using standard Yee-cell staggering of the field components [42]. The time step is chosen according to the Courant criterion as t = x/ √ 3c, with c being the vacuum speed of light. The width of the Gaussian PIC particle w pic can be chosen similar to the grid spacing. For current and force calculation we use the current and charge density of PIC particles at the grid points. Only grid points closer than 3.5w pic to the PIC particle are evaluated, which conserves 99.998% of the total particle charge. The laser field with a flat beam profile is treated via the total-field-scattered-field scheme [42]. Hence, only the scattered fields have to be propagated on the PIC grid. On the surface of the simulation box, we apply uniaxial perfectly matched layer (UPML) absorbing boundary conditions [42] to remove scattered fields.
The microscopic correction represents the particle-particle part of the P 3 M decomposition and is implemented as a local MD. Since the microscopic correction forces are short range, the corresponding binary forces need to be evaluated only within a finite cut-off sphere with radius r cut . A cut-off radius of r cut = 3w pic turned out to be sufficient for the present study. To 8 evaluate the microscopic correction within the respective cut-off spheres, neighboring particles within the cut-off radius r cut have to be identified. For an efficient identification we apply the cell index method [43], resulting in linear scaling with particle number N per step. Note that a direct treatment (comparison of all particle pairs) results in an N 2 scaling of the identification procedure and would destroy the linear scaling of MicPIC. The system is advanced after the determination of the Mic and PIC forces for each time step.
PIC particle size and PIC grid spacing have to be chosen fine enough to resolve the length scale of field propagation processes. Apart from that constraint, x and w pic can be chosen freely. This freedom can be utilized for load balancing between Mic and PIC parts, which can be seen from the following scaling relation. It should be noted that there is no limitation on the number of PIC particles per grid cell. By exploiting the fact that the cut-off radius, PIC particle width and grid spacing are chosen in a fixed relation to each other, the total numerical effort for the propagation of a fixed time interval can be expressed as [34] where the parameters α, β and γ are constant prefactors for the microscopic correction, the current injection and force evaluation on the PIC level and the solution of Maxwell's equation on the PIC level, respectively. For very small/large values of r cut , the PIC or microscopic parts produce the dominant numerical load. Hence, the cut-off radius (or equivalently the PIC particle width) determines the load balancing between the microscopic correction and the PIC part. It can thus be chosen to minimize the total numerical load. This typically leads to values for w pic that are smaller than the ones needed to resolve the physical processes. In our study, we used a 300 3 -cell PIC grid with mesh size x = 5 Å, UPML boundary conditions [42], a PIC particle width w pic = 1.15 x and a cut-off radius r cut = 3w pic . This results in almost equal numerical load between Mic and PIC parts of the code for the largest investigated systems. Note that the numerical effort is independent of actual particle widths, which is of the order of w 0 = 1 Å for electrons and ions in our study. The particular values are given in the discussion of the considered scenarios.
For the response experiments, scattered and absorbed energies have to be determined. The total scattered radiation energy is measured using the Larmor formula [44] wherep(t) denotes the time-dependent dipole acceleration. The energy absorption by the system contained in the box follows from the total energy difference before and after excitation. The total box energy reads where the individual terms on the right-hand side describe the kinetic energy, the electromagnetic energy on the PIC level, the energy resulting from the microscopic correction and the energy renormalization to remove the spurious self-energy of the particles on the PIC grid. 9

Molecular dynamics
As an electrostatic reference in the small cluster limit and for the identification of propagation effects via comparison to MicPIC, we use an MD code [32]. The representation of the actual plasma particles in MD and MicPIC is the same; we use Gaussian charge distributions with width w 0 . The difference is that in MD the electrostatic fields obtained from these charge distributions are used throughout the total simulation volume. MD takes into account the electrostatic interactions between all particles and includes the external laser field in dipole approximation, resulting in the usual N 2 scaling of the numerical effort. The MD force on the ith plasma particle is given by and is evaluated using massively parallel computation techniques. In MD, the energy absorption is given by the total energy difference before and after the laser excitation. As in the MicPIC case, the scattered energy is determined by the Larmor formula. It should be stressed that the determination of absorption and scattering cross-sections with MD is only valid in the limit of weak light scattering, as the energy loss due to the re-emission of light is not taken into account.

Continuum models
To determine the absorption and scattering of the pre-ionized metal-like clusters with continuum theory, the nanoplasma is treated as a dielectric sphere with radius R and frequency-dependent relative dielectric constant ε(R, ω). The dielectric response of the nanoplasma electrons is described by a Drude-like dielectric function with χ 0 a real-valued background susceptibility, ν the collision frequency and ω p = e 2 n e /m e ε 0 the plasma frequency. Note that the frequency of the Mie plasmon of a small, perfectly metallic sphere is ω mie = ω p / √ 3. To include both electron-ion and electron-surface collisions in the cluster response, we use a cluster radius-dependent collision frequency ν(R) = ν 0 + ν 1 /R. Here ν 0 is the regular bulk collision frequency and ν 1 /R describes the additional surface collisions [45], where ν 1 is an effective electron velocity.
The microscopic effects contained in the dielectric function are accounted for by collision frequency and susceptibility. Once these parameters are known, the dielectric function can be used together with macroscopic models, such as Mie and nanoplasma theory, to model laser absorption and scattering by a cluster.
In Mie theory, the frequency-dependent absorption and scattering cross-sections are defined as [44,46] and Here the a n and b n are the well-known Mie coefficients; see [44,46] for details of their calculation. The nanoplasma model is the electrostatic small-sphere limit of Mie theory (see p 136 of [46]) and neglects propagation effects. In this electrostatic limit (indicated by superscript 'stat') the single-frequency absorption cross-section becomes Note that this result is equivalent to the heating rate used in the well-known nanoplasma model [47]; therefore we will denote the small-sphere electrostatic limit as nanoplasma theory henceforth. The corresponding small-sphere expression for the scattering cross-section is The MD analysis and the nanoplasma model do not account for propagation effects and are only applicable to the small cluster limit. However, as MD and the nanoplasma model are based on the same (electrostatic) approximation, they describe the same physical scenario and must fit to each other even outside their range of validity. The Mie theory and MicPIC model include propagation and are valid for all cluster radii. However, the validity of Mie theory is confined to the linear response regime. Finally, the microscopic parameters of the dielectric function in continuum theory, such as the collision frequency, have to be determined. The dielectric function for both continuum approaches must be the same.
In order to compare the Mie and nanoplasma models with MicPIC and MD for excitation scenarios with short laser pulses, the cross-sections need to be averaged over the spectrum of the laser pulse. We consider Gaussian pulses with carrier frequency ω 0 , duration τ , field envelope exp(−t 2 /τ 2 ) and the corresponding normalized intensity spectrum I (ω) = τ exp(−[ω − ω 0 ] 2 τ 2 /2)/ √ 2π . The effective cross-sections then follow from Normalizing the cross-sections to the geometrical cross-section of the sphere yields the single frequency and spectrally averaged efficiency factors for absorption and scattering.

Results
In the following, we present the results for (i) the linear and nonlinear excitations of resonant pre-ionized clusters with τ = 7 fs few-cycle laser pulses with Gaussian field envelope exp(−t 2 /τ 2 ) and (ii) the nanoplasma formation in an initially neutral Ar cluster. Beginning with scenario (i), the model cluster for the MicPIC and MD analysis is a homogeneous spherical droplet with one conduction electron per ion, initialized with an electron temperature of T e = 5 eV. The width of the particles is w 0 = 1 Å. Further ionization of bound atomic electrons is disregarded. We have chosen face-centered cubic (fcc) ion geometry with an ionic Wigner-Seitz radius of r s = 3.6 Å. The resulting cluster has an electron density of n e = 5.1 × 10 21 cm −3 and is in resonance (Mie plasmon) with a near-infrared laser wavelength of 800 nm. In terms of the critical plasma density n crit = ω 2 las m e 0 /e 2 , this scenario describes a spherical plasma that is three times overcritical and has a skin depth of 73.5 nm. Although our cluster is a model system, it also has experimental relevance. It describes a fully ionized (metallized) and pre-expanded hydrogen nanodroplet, which could be prepared, e.g., via the pump pulse in a pump-probe experiment.

Linear-response dynamics for small and large clusters
The linear-response regime of our resonant nanoplasma system is studied with a laser intensity of I = 6 × 10 11 W cm −2 . The time evolution of selected observables, as predicted by MicPIC and MD simulations, is displayed in figure 1 for a small (R = 10 nm) and a larger resonant cluster (R = 30 nm).
Both the small and large cluster results show the excitation of pronounced dipole oscillations of the cluster electrons. The resonant nature of the interaction is reflected by the increase of the dipole amplitude even after the pulse peak. The decay of the oscillations after the laser pulse reflects damping of the collective motion via electron-ion and electron-surface collisions, as well as radiation damping. While collisional damping converts energy from the collective motion into thermal nanoplasma energy, part of the collective excitation energy is removed from the system via the emission of light. Signatures of both effects can be found in figure 1.
Firstly, the conversion of collective excitation energy into heat can be directly extracted from the decay of the mutual transformation of potential and kinetic energies; see figures 1(c) and (d). The collective center of mass motion of the electron cloud periodically converts potential into kinetic energy and vice versa. Once thermalization is completed, the energy stored in the collective motion has been fully transformed to heat.
Secondly, the importance of radiation effects as a function of cluster size is revealed by the 10 and 30 nm simulations. Whereas the dipole signals and energies calculated by MD and MicPIC for the 10 nm cluster are in almost perfect agreement, they substantially differ for the 30 nm system. The total cluster energy after the laser pulse is higher for MD than for MicPIC, as MicPIC accounts for the attenuation of the laser field in the cluster due to the skin effect. Note that the cluster diameter is comparable to the skin depth. Further, the MicPIC dipole signal amplitude reveals a faster decay. Both differences are due to energy lost from the cluster by radiation. Finally, the frequency of the dipole oscillation after the pulse is slightly lowered in the MicPIC data, which can be attributed to the polaritonic red-shift of the Mie plasmon via field retardation [45]-another propagation effect. The difference between MicPIC and MD final total energies is of the order of the total absorption and impressively demonstrates the large error incurred by a pure electrostatic treatment of the 30 nm system. The neglect of propagation effects limits MD analysis to small cluster sizes.
The excellent agreement of MD and MicPIC for 10 nm shows that microscopic effects are resolved in the same way in both codes. As MD is a well-established numerical method, it proves that collisional effects in MicPIC are correctly accounted for.
Finally, the comparison of the two cluster sizes further demonstrates the radius-dependent effect of the cluster surface on the optical response. Both the decay time of the dipole amplitude In the bottom panels, E kin , E pot and E tot denote the total kinetic energy, the potential energy and the total energy of the cluster, respectively. E kin,col shows the kinetic energy of the center of mass motion of the electron cloud. For better comparability all data are normalized to the number of electrons N . The given plasmon energies and lifetimes are the results from fitting a damped harmonic oscillator to the decaying dipole moments (fit regions as indicated). and the total energy absorption per electron are clearly enhanced in the larger cluster system, as can be most clearly seen by comparing the MD results. This shows the reduced damping effect of electron-surface collisions due to the smaller surface to volume ratio in larger clusters.

Radius-dependent absorption and scattering: comparison with continuum models
For a systematic analysis of the propagation and surface effects, we performed a radiusdependent comparison of the MicPIC and MD results with the predictions of the Mie theory Absorption (Q abs ) and scattering (Q sca ) efficiencies of metallic clusters excited at resonance (800 nm) by 7 fs laser pulses with a peak intensity of 6 × 10 11 W cm −2 . The efficiency is defined as the absorbed and scattered energies normalized by the pulse fluence and geometrical cluster cross-section. Part of the data (MD absorption; MicPIC absorption and scattering up to 40 nm) are taken from [34]. and nanoplasma model. In order to make this comparison possible, the parameters in the dielectric function (18) used in the continuum models have to be determined from the absorption and scattering efficiencies obtained from MicPIC and MD simulations. The absorption and scattering efficiencies, as obtained from the MD, nanoplasma model, MicPIC and Mie theory, are plotted in figure 2. MD data points are given only up to 30 nm, as this was the maximum size we could simulate with the available computation resources. Part of the data (MD absorption; MicPIC absorption and scattering up to 40 nm) has been reported previously in [34]. Figure 2 reflects our ongoing effort to model increasing cluster sizes with MicPIC. The goal is to test the validity of continuum models over a wide range of cluster sizes. Up to R = 50 nm cluster radius, the quality of the Mie theory based on the radius-dependent dielectric function is found to be outstanding and justifies its applicability in section 2.3.
As the excitation takes place in the linear response regime, MicPIC results for any cluster size can be compared with Mie theory [44], which fully accounts for propagation effects (including skin effect, field attenuation and radiation damping). In the small cluster limit, propagation effects are negligible and MD describes the classical cluster dynamics exactly. In this electrostatic limit, Mie theory simplifies to the nanoplasma model. The nanoplasma model can be compared to MD, as it uses the same approximations. Matching the nanoplasma model to MD results yields the coefficients for the real-valued background susceptibility χ 0 = 0.15, the bulk collision frequency ν 0 = 0.102 fs −1 and the effective electron velocity ν 1 = 2.4 nm fs −1 as 14 a measure of the surface collision contribution to the optical response. Using these parameters in the Mie model yields excellent agreement between MicPIC and Mie theory for the entire size range considered and validates MicPIC for large systems. The agreement of the microscopic simulations with the respective continuum model based on the same dielectric function proves the consistency and, in particular, the validity of MicPIC.
The impact of field propagation can be extracted directly from the comparison of the electrostatic and electromagnetic descriptions. For larger clusters, i.e. in the size regime R 20 nm, the MicPIC absorption is substantially reduced when compared to the prediction of MD. This reflects the onset of field attenuation and radiation damping as clear signatures of field propagation effects. Values of Q abs beyond unity for R 8 nm imply absorption cross-sections higher than the geometric cross-section of the sphere. In analogy to the extinction paradox [48], we refer to this phenomenon as the absorption paradox. The maximum absorption efficiency for R ≈ 37.5 nm indicates an optimal cluster size for efficient resonance enhanced absorption for the investigated pulse parameters. For substantially larger clusters, the strong increase of field attenuation and radiation damping leads to a rapid decrease of the absorption efficiency.
The main result obtained from the comparison in figure 2 is the set of macroscopic material parameters in the dielectric function. These parameters quantify the effect of microscopic processes in continuum models. Without knowledge of their magnitude, continuum models are confined to qualitative predictions. So far, experimental measurements or theoretical determination of optical material parameters for nanosystems has proven difficult. Here, MicPIC opens up new capabilities. In particular, when two or more different processes contribute to a parameter, such as bulk and surface collisions in a cluster or nanoparticle, it is very difficult to unravel the various processes and determine their relative importance. For example, the quantitative description of surface versus bulk collision effects in metallic nanoparticles and nanoplasmas based on first principles and including the full electromagnetic response has not been possible so far. The new insights that can be gained from the direct determination of the optical parameters will be demonstrated below in terms of the plasmon lifetime.

Competition of bulk and surface effects with radiation damping in resonant clusters
Using the parameters of the dielectric function determined above, the relative importance of surface and bulk effects can be studied quantitatively as a function of the cluster radius. Here we concentrate on the size dependence of the absorption of the dipolar Mie plasmon. Therefore, we compare the absorption efficiencies from Mie theory using the dielectric function with and without the surface collision term. The second case (ν 1 = 0) corresponds to the bulk limit. Figure 3 shows the radius-dependent spectral profiles of the absorption efficiencies. We also indicate the intensity spectrum of the laser pulse used in the above analysis for the calculation of spectrally averaged efficiencies. Inspection of the absorption efficiencies in figures 3(a) and (b) shows that the peak efficiency (solid line) exhibits a pronounced polaritonic redshift with increasing particle size, irrespective of the surface damping term. However, the radius dependence of the peak value and width of the absorption efficiency profiles shows a clear sensitivity to surface effects up to about a radius of 80 nm. In the absence of the surface term, the efficiency profiles increase more rapidly with radius for small particles, are narrower, and reach a substantially higher maximum peak value near 40 nm when compared to the prediction including the surface contribution; compare figures 3(a) and (b). This obvious deviation of the absorption properties highlights the strong influence of surface collisions. An experimentally relevant and accessible parameter is the width of Mie plasmon resonance. To analyze its size dependence, the linewidth has been determined from a Lorentzian fit of the absorption profiles. The resulting plasmon lifetime (inverse linewidth) is displayed in figure 3(c) for both versions of the dielectric function. The following conclusions can be extracted from the comparison.
Including the bulk collision term only, the lifetime is maximal for small particles and decreases monotonically with increasing cluster radius due to radiation damping. In contrast to that, the lifetime becomes peaked with the inclusion of the surface terms. The reduction for small particles directly reflects the surface damping, which is increasingly important here because of the high surface to volume ratio. For large particles, the full result converges asymptotically to the bulk behavior. The maximum lifetime thus results from the competition between the surface and bulk collisions with radiation damping.
Maximal lifetimes of the Mie plasmon in the few tens of nanometers size range have been observed experimentally for spherical noble-metal nanoparticles [49][50][51]. Both the measured optimal particle size and the peak lifetime are of the same order of magnitude as our theoretical result. Deviations in the range of 50% are attributed to our model system. A more detailed comparison will require the modeling of the particles at solid density, which will be a next step in the application of MicPIC.

Nonlinear plasma wave dynamics in large resonant clusters
Next, we investigate resonant cluster excitation at a higher intensity (6 × 10 13 W cm −2 ), where nonlinear processes become important. Our model system describes a fully ionized hydrogen plasma droplet that has been expanded to fulfill the resonance condition, which could be realized in a pump-probe experiment. Except for the higher intensity, the same laser parameters are used as in figure 1. In a recent work [34], the creation and evolution of plasma waves have been identified in the interaction of noble gas clusters with intense laser pulses. Electric field enhancements caused by the plasma wave dynamics result in the increased ionization of deeper bound electrons, which in turn influences the plasma wave dynamics. Here we study the same process, but in a purer system, where the plasma wave dynamics is not tainted by ionization processes. Although qualitative results are similar, the signatures of the plasma wave dynamics come out more clearly, making a hydrogen droplet the ideal system for a first experimental demonstration of plasma wave dynamics in clusters.
Selected charge density snapshots of the dynamics of an R = 30 nm cluster are shown in figures 4(a)-(f); the corresponding time instants are put in relation to the laser field and the evolving dipole moment (see the inset). The analysis of the plots shows that plasma waves are excited at the cluster surface and propagate into the cluster (figures 4(a) and (b)). The waves are mainly produced at the cluster poles, i.e. where the x-axis intersects the cluster surface. With increasing dipole amplitude, stronger plasma waves are created and penetrate deeper into the cluster ( figure 4(c)). The waves are a mixture of spherical and plane waves, as they evolve in the mixed geometry of a sphere and a linear cluster polarization field. When the laser intensity is high enough that the generated plasma waves reach the cluster center and collide (figure 4(d)), wave breaking and turbulent electron dynamics are observed (figures 4(e) and (f)). The resulting fluctuations of the charge density evolve on the attosecond time scale with peak amplitudes of the order of the ion density and spatial structures on the nm scale.
The origin of the plasma wave dynamics can be understood from the electron phase-space distributions f (x, v x ) (see figures 4(g)-(l)). The snapshots correspond to the same times as the density plots above. Figure 4(g) shows a tail on the left side of the cluster that was created in the previous half-cycle, when the electron cloud was pushed over the left boundary of the ion sphere. At the moment of the snapshot, the main electron cloud has already moved to the right, leaving behind unscreened ion charges at the left border of the cluster (see also the yellow surface region in figure 4(a)). The resulting polarization creates a strong force that accelerates a part of the escaping electrons back into the cluster. When they hit the surface, a plasma wave is created (figure 4(h)) that propagates into the cluster together with the recolliding electron bunch (figure 4(i)). The excited plasma wave appears as the dip in the (yellow/green) cluster electron distribution in figure 4(i). At higher plasmon amplitude, the generated particle-wave units copropagate deeper and faster into the cluster; see also the ring structure in figures 4(a)-(d). In plasma physics the above-described process is known as wake field acceleration; in clusters, it has not been observed so far. Wave breaking occurs when the plasma waves from the left and right cluster surfaces reach the center and collide (figure 4(j)), leading to strong fragmentation of the electron phase space distribution and high local electric fields (figures 4(k) and (l)). Note that both the density snapshots and phase-space plots show much clearer signatures of the plasma waves and recolliding electron bunches when compared to the previous results in [34] on resonant Xe clusters, where the inclusion of further atomic ionization in the cluster leads to additional damping. Figure 5 provides more details of the dynamics inside the 30 nm cluster. In particular, it shows plasma waves being created at the left and right boundaries (see figure 5(a)). In the rising edge of the laser pulse (t < −2 fs), the weak plasma waves die early without propagating far into the cluster. The corresponding field evolves close to that of a polarized metallic sphere, showing harmonic oscillations with a nearly homogeneous profile inside the cluster ( figure 5(b)). However, the propagation and collision of the increasingly strong plasma waves being created near and after the pulse peak completely change this behavior and induce strong spatial and temporal density and field fluctuations; see figures 5(a) and (b). At the cluster center, the electric field intensity (∝ E 2 ) exceeds the laser peak intensity by more than two orders of magnitude and changes on the attosecond time scale.
Finally, it should be noted that the generation of plasma waves in small clusters could, in principle, be observed also using MD techniques. However, plasma wave dynamics becomes less prominent with decreasing cluster radius, which could be the reason why it has not been observed in MD simulations so far. As MicPIC can be more easily scaled to larger cluster sizes and accounts for wave propagation, it is better suited to the investigation of this phenomenon.

Nanoplasma formation in a small rare-gas cluster
So far we have considered pre-ionized, expanded nanoplasmas and have disregarded atomic ionization. In the last part of the paper, we focus on describing the full plasma dynamics of an initially neutral rare-gas cluster exposed to an intense laser field. The main purpose of this section is to demonstrate that MicPIC is capable of describing the full plasma dynamics including ionization. For this benchmark scenario, we use a spherical R = 5 nm Argon cluster (N = 11100) with fcc structure and atomic Wigner-Seitz radius r s = 2.21 Å. A width parameter of w 0 = 0.81 Å is used to prevent classical recombination below the quantum mechanical energy levels. For cluster excitation a 25 fs laser pulse (800 nm) with a peak intensity of 1 × 10 15 W cm −2 is considered. The atomic ionization dynamics including tunnel and impact ionization is implemented according to Fennel et al [32] and includes the effects of ionization threshold lowering due to local plasma fields. It should be noted that the ionization is described with full atomic scale resolution, i.e. each atom/ion is tested for ionization separately. Since propagation effects are negligible for this cluster size, the results can be compared with MD to validate MicPIC. Figure 6 shows the cluster dynamics calculated with both methods and demonstrates excellent agreement.
Focusing on the physics of this benchmark scenario, the following conclusions can be extracted from figure 6. First, plasma formation starts with tunnel ionization in the leading edge of the laser pulse, see figure 6(c). The threshold behavior reflects the high nonlinearity of the tunnel ionization rate as a function of intensity. The ignition due to tunnel ionization produces only a few electrons, which are then efficiently heated and launch an impact ionization avalanche, rapidly producing a highly overdense nanoplasma (about 60 times overcritical), cf figures 6(b) and (c). Note that electron impact ionization is by far the dominating ionization channel, contributing about 90% of the ionization events. Because of the overcritical nature (the cluster efficiently screens the external laser field), heating due to inverse bremsstrahlung is expected to proceed mainly at the cluster surface (surface heating). During the pulse, the nanoplasma electrons are heated moderately, leading to an electron temperature of about 20 eV at the end of the pulse. The high electron density combined with the moderate temperature result in a strongly coupled nanoplasma. The evolution of the coupling is indicated in figure 6(d) in terms of the Debye number N D = 4π 3 n e λ 3 d , where n e and λ d are the electron density and the Debye screening length extracted from the simulations. The Debye number quantifies the average number of electrons within the Debye sphere and indicates strong coupling for N D 1. The calculated evolution demonstrates that a strongly coupled plasma is formed shortly after the initial ionization. The fact that a strong coupling persists for the rest of the dynamical evolution  In panel (a), E tot , E kin,tot , E kin,ion and E kin,el denote the total absorption of the cluster, the total kinetic energy, the ion kinetic energy and the electron kinetic energy, respectively. requires a fully correlated treatment of the microscopic collisions to describe the dynamics correctly. Note that the concept of Debye screening as a linear screening theory breaks down for strong coupling and is used here for the identification of the coupling regime only.
Focusing on the ion dynamics, cluster expansion begins near the pulse peak and leads to energy conversion of electron thermal energy into ion kinetic energy via adiabatic expansion cooling (hydrodynamic expansion); see figures 6(a) and (b). However, the increase of the total kinetic energy (electrons and ions) shows that also Coulomb repulsion due to outer cluster ionization (Coulomb explosion) and charge spill-out at the surface contributes to the expansion. The cluster expansion begins to efficiently lower the nanoplasma density, which is about 20 times overcritical at the end of the simulation. After substantial further expansion, the cluster is expected to reach the conditions for Mie resonance, see the indicated density (3n crit ) for Mie plasmon resonance.
The excellent agreement with the MD results demonstrates that MicPIC correctly accounts for microscopic processes (ionization and collisions) and is capable of efficiently modeling strongly coupled nanoplasmas. The comparison of the run times for both methods (1 day on 512 cores for MD; 14 days on a single core for MicPIC) highlights the striking numerical performance of MicPIC over MD already for this moderate system size. For larger systems, which will become accessible after parallelization of the MicPIC code, the numerical advantage will be even more pronounced. Most importantly, because of the inclusion of field propagation effects, only MicPIC will be applicable.

Conclusions
In this study, we utilized MicPIC to investigate the resonant interaction of few-cycle laser fields with cluster nanoplasmas in the linear and nonlinear regimes. We demonstrate the potential of MicPIC to bridge the realms of microscopic and macroscopic processes by using the example of laser-driven clusters. Our analysis sheds light on both microscopic collisional processes and macroscopic field propagation effects.
In the linear regime, we have shown the potential to determine material parameters for systems that are strongly influenced by finite-size effects. Their detailed account is key to modeling the optical properties of nanostructures with continuum theory. As the treatment of surface effects requires atomic-scale resolution, microscopic collisions and macroscopic propagation effects need to be treated simultaneously to correctly model laserdriven nanostructures. In this work, this capability of MicPIC has been demonstrated for the plasmon lifetime. Our results reveal a pronounced maximum lifetime of the Mie plasmon of the order of 10 fs for cluster radii of a few tens of nm, which is in reasonable agreement with experiments for spherical noble-metal nanoparticles. We find that the maximum in the lifetime is created by the competition of bulk and surface collision effects with radiation damping. Our analysis even allows us to separate the two collisional contributions to the plasmon lifetime. In a next step, more refined calculations on noble-metal clusters will be performed for the direct quantitative comparison with experiments.
The capacity to determine material parameters including finite-size effects will be of relevance for analyzing processes in nano-photonics and nano-plasmonic structures [52], which can have a thickness of a few atomic layers only. There, resolution on an atomic length scale is likely to be important, and currently used macroscopic models, treating materials on a continuum level, might fall short.
In the nonlinear regime, we have studied the excitation of plasma waves in resonant hydrogen-like clusters with R = 30 nm. Compared to previous results on pre-ionized, resonant Xe clusters of the same size, the plasma wave dynamics comes out more clearly, making a small hydrogen droplet an ideal system for a first experimental demonstration of plasma waves in clusters. The resulting plasma wave dynamics takes place on a nanometer length and an attosecond time scale, opening up a new area for the investigation of attosecond dynamics. Furthermore, the regular wave to turbulent dynamics transition is of fundamental interest for fluid dynamics and plasma physics. MicPIC enables a rigorous analysis of turbulence with microscopic/atomic resolution, which is not feasible with existing methods. Combined with ultrafast x-ray imaging [53,54], we expect new insights into the making and breaking of waves in finite plasmas.