A concept to generate ultrashort ion pulses for pump-probe experiments in the keV energy range

The impact of an energetic particle onto a solid surface generates a strongly perturbed and extremely localized non-equilibrium state, which relaxes on extremely fast time scales. In order to facilitate a time-resolved observation of the relaxation dynamics using established ultrafast pump-probe techniques, it is necessary to pinpoint the projectile impact in time with sufficient accuracy. In this paper, we propose a concept to generate ultrashort ion pulses via femtosecond photoionization of rare gas atoms entrained in a supersonic jet, combined with ion optical bunching of the resulting ion package. We calculate the photoion cloud generated by an intense focused laser pulse and show that Arq+ ions with q = 1–5 can be generated with a standard table-top laser system, which are then accelerated to energies in the keV range over a very short distance and bunched to impinge onto the target surface in a time-focused manner. Detailed ion trajectory simulations show that single ion pulses of sub-picosecond duration can be generated this way. The influence of space charge broadening is included in the simulations, which reveal that flight time broadening is insignificant for pulses containing up to 10–20 ions and starts to increase the pulse width above ∼50 ions/pulse.


Introduction
The interaction of ions with solids is of considerable interest for many fields such as atomic physics, solid state physics and chemistry as well as materials research. Ions can be used to modify either the surface or the bulk of the irradiated material by changing the electronic properties (e.g. by implanting ions for doping or defect engineering) [1][2][3][4], cleaning, etching or patterning the surface by sputtering (i.e. the removal of surface material) [5][6][7][8], fabrication of thin films (by collecting the sputtered material on a substrate) [9,10], imaging [11], or chemical surface analysis via mass spectrometry of the sputtered material [12][13][14].
If an energetic ion hits a solid surface, it deposits its kinetic impact energy within the target material. The choice of mass, composition, charge state and kinetic energy of the ions affects the energy density which is deposited in a sample, thereby offering a way to choose between different energy transfer mechanisms [15,16]. At high impact energy (>1 MeV/u), the primary energy loss mainly occurs to the electronic system of the solid, thereby generating a strongly localized electronic excitation which may then transfer energy to the target nuclei via electron-phonon coupling. In the limit of low impact energies (<10 keV/u), on the other hand, the primary energy dissipation occurs via mostly elastic collisions with the target atoms, leading to fast, billard-like collision dynamics which may then couple energy into the electronic system via electronic stopping of moving particles. In any case, the projectile impact generates a strongly perturbed non-equilibrium state, which is extremely localized both in space (nm) and time (sub-picoseconds). Up to date, the relaxation dynamics following such a sudden excitation are experimentally inaccessible.
In view of the rapid advancement of techniques investigating ultrafast processes via pump-probe experiments, the question arises whether such methods could be used in order to gain an experimental handle on ion impact-induced relaxation dynamics as well. Most of the currently used pump-probe experiments utilize

Methods
The principal setup of the proposed ion source is depicted in figure 1. A supersonic jet of rare gas particles is crossed with an intense femtosecond laser pulse in order to produce positive ions of the desired species and charge state. The resulting ion packet is extracted towards the target surface via a two-stage electric field generated in a three-electrode bunching configuration, with the target surface itself acting as the third electrode. The two remaining electrodes are set to suitable potentials to ensure first-order flight time focusing conditions, where the generated ions impinge onto the target surface at nearly the same time. Depending on the operation conditions of the supersonic expansion, single atoms as well as clusters of several atoms can in principle be formed and ionized. Moreover, by varying the intensity of the ionization laser, ions of different charge states can be produced. The targeted impact energy onto the sample surface is chosen as 5 keV for singly charged ions. In order to keep the total ion flight time short and minimize the temporal spread, the entire ion source is strongly miniaturized, with the distance between two electrodes being as small as possible.

Supersonic jet
The supersonic jet is formed via adiabatic expansion from a pulsed nozzle. In order to permit a high repetition rate, a commercially available piezo-driven valve will be used which can be operated at pulse lengths down to ∼10 μs and repetition rates up to 5 kHz [26]. The nozzle will be operated with argon or neon gas at backing pressures p 0 between 1 and 10 bar and placed at a distance of ∼10 cm before a first skimmer of 1.5 mm entrance diameter. In the context of the present work, important parameters of the resulting jet are the density and temperature of the gas particles (atoms or clusters) in the beam center at a given distance from the nozzle. According to theory, the density n should scale as a function of the distance d as [27,28]  ) [27], yielding d ref = 6×10 −3 cm with a nozzle diameter of 2R nozzle = 150 μm. Placing the ionization region at a distance of d= 100 cm from the nozzle, one would therefore expect a beam density of the order of the order of several 10 11 cm −3 at the point of interaction with the ionizing laser pulse.
Experimentally observed densities, on the other hand, are usually found to be about one order of magnitude smaller than the theoretical values due to clogging of the skimmers [29]. This phenomenon becomes more significant with decreasing skimmer aperture size and increasing gas density in the skimmer. The resulting transmission was calculated by Luria et al [29]. For a 1.5 mm diameter skimmer, they find a transmission of 15% at a gas density of 3×10 15 cm −3 for argon [30]. According to equation (1), we expect a beam density of about 1.5×10 15 cm −3 at the position of the first skimmer (100 mm away from the nozzle), so that this value should represent an upper limit to the beam attenuation induced here. There is also a heating effect induced by the skimmer walls, which, however, should be small (<0.5 K) provided a sufficiently sharp skimmer with 3 μm wall thickness is used [29]. In order to set up an additional differential pumping stage and define the final beam size, further skimmers will be placed at larger distances downstream from the nozzle, where the density should already be small enough to prevent significant skimmer interference. More specifically, a second skimmer with 1 mm entrance aperture will be located at d = 200 mm before a second differential pumping stage. A third skimmer of several 100 μm entrance aperture will then be used to define the final beam size upon entering the ultrahigh vacuum chamber containing the ion source.
From the above considerations, we expect a gas density of the order of 10 11 atoms per cm 3 in the ionization region. Using xenon gas at a backing pressure p 0 =2 bar and a conical nozzle of 500 μm diameter, which was pulsed with a 300 μs opening time, Schofield et al [31] have measured the density profile of the resulting xenon beam via laser ionization tomography and found a center density of 4×10 12 cm −3 at a distance of d = 301 mm from the nozzle. In their experiments, they also used two skimmers located at d = 50 mm (1 mm aperture) and at d = 280 mm (100 μm diameter), so that the beam conditions should be comparable to the ones projected here. Also in this case, the measured density is by about a factor 10 lower than that expected by theoretical considerations. These results lend credibility to the expected beam density estimate given above.
An important aspect in the context of the present work is the fact that the gas is effectively cooled during a supersonic expansion. Using the sudden freeze model, it is assumed that the particles undergo collisions, until at a certain 'freezing distance' d F from the nozzle their density has decreased enough to effectively isolate their motion. Beyond d , F the particles basically follow straight stream lines, which can be traced back to a virtual source of radius R . source Up to the freezing distance, a temperature can be assigned to the particles, which is the same for their relative motion both along (T  ) and perpendicular (T^) to the beam propagation direction. The theoretical description of T  and T^has been discussed in detail by Beijerinck and Verster [32]. Both temperatures are identical up to the freezing distance and remain constant afterwards. On the other hand, the velocity components v^associated with T^lead to diverging particle trajectories, where particles with larger vd iverge more strongly. Therefore, if only the center part of the beam is probed at a distance d downstream, the effective 'temperature' T eff of the remaining beam particles is reduced with increasing flight length. If the probed volume is restricted to a diameter R 2 in the direction perpendicular to the beam, the radial velocities of the ionized particles are geometrically limited to v v R d.


Here, v k T m 5 B 0 =  is the terminal beam velocity, which is calculated as 568 m s −1 for argon. In our concept, the beam will be collimated to a diameter of less than a millimeter by a third skimmer located at a distance of about 50 cm from the nozzle. Assuming R = 0.5 mm in combination with d = 50 cm, this restricts the possible starting velocity component of ions along the extraction direction perpendicular to the beam axis to values smaller than 0.6 m s −1 . If converted to energy, this corresponds to 7×10 −8 eV or ∼1 mK. This so-called 'geometric cooling' effect represents the key to the short pulse ion source concept presented here. Using a more sophisticated approach to describe the (in principle non-Maxwellian) velocity distribution yields [32] Calculation of d F requires the potential describing the interaction between the particles [32]. Measurements of the axial beam temperature yield values of T T d F = =  ( ) 2.3 K and 0.6 K for Xe [31] and Ar [30,31], respectively. With d F » 8 mm [31] and d= 50 cm, this results in T eff < 1 mK.
It should be noted that the above estimate neglects the skimmer influence, which will inevitably heat the beam and increase the beam divergence. In that context, the most critical device is the third skimmer, since possible increases in beam divergence generated at the first two skimmers will be remedied by the geometric cooling effect. As a worst possible case, the virtual source may be placed at the position of the third skimmer. Due to the fact that the laser beam is tightly focused (see below), the relevant ionization volume is restricted to a dimension below 100 μm in the direction perpendicular to the supersonic beam propagation, thereby again leading to a geometric cooling effect which limits the maximum radial velocity to < 0.3 m s −1 . In summary, the radial starting temperature of the ionized atoms will therefore be restricted to values of the order of 1 mK. In the remainder of this paper, we will nevertheless assume a value of 10 mK in order to arrive at an upper estimate for the role played by the initial motion of the generated ions.

Photoionization
Ionization of the gas particles (atoms and possibly clusters) entrained in the supersonic beam will be accomplished via strong-field photoionization in an intense ultrashort laser pulse. For that purpose, we intend to use a chirped-pulse-amplification titanium-sapphire laser operated at a wavelength of 800 nm, which delivers pulses of up to 2 mJ energy and<50 fs duration at a repetition rate of 1 kHz. Since the photon energy of 1.5 eV is small compared to the ionization energy of argon (15.6 eV), efficient ionization can only be achieved via multiphoton or field ionization processes which require relatively high peak power densities in excess of 10 14 W cm −2 . As a consequence, the laser beam will be focused to a spot size of<10 μm within the extraction volume probed by the ion buncher.
The resulting ionization probability p i q ( ) for production of a q-fold positively charged ion can be estimated using, for instance, Ammosov, Delaney and Krainov (ADK) tunnel ionization theory [33], which has been shown to describe the strong field ionization of rare gas atoms quite well [34]. Applying the ADK model to argon, we obtain the laser intensity dependence of single and multiple ionization probabilities as shown in figure 2. For reference, the result of a calculation based on perturbative multiphoton ionization has been included as a dashed line, revealing that the exact form of the laser intensity dependence is not critical. It is seen that singly charged Ar + ions are efficiently produced at peak power densities around several 10 14 W cm −2 . If the laser intensity is increased further, higher charge state ions are predominantly generated, so that Ar q+ ions with q=1-6 can in principle be produced with the available laser system.
In order to estimate the number and spatial distribution of the generated ions, the photoionization probability resulting from the laser intensity distribution must be convoluted with the atom density distribution of the supersonic neutral beam. Since the laser beam will be focused to a small spot size of a few μm diameter within the ionization region, the extension of the effective ionization volume is small in the direction perpendicular to the laser beam. In the direction along the laser beam, on the other hand, the ionization volume is in principle determined by the Rayleigh range of about 100 μm. It is further ion-optically restricted by a small (80 μm f) aperture in the extraction electrode (see below). For that reason, it would be sufficient to limit the diameter of the supersonic jet to about 100 μm as well, using the third skimmer at a distance of about 20 cm before the ionization region. In practice, we will use a larger skimmer aperture in order to homogenize the atom density distribution, which will in the following be assumed constant throughout the entire effective ion extraction volume.
The ionization laser features a nearly Gaussian beam profile of about 14 mm diameter (full width half maximum, FWHM) prior to focusing. The beam will be deflected by 90°and focused into the ionization chamber by an off-axis parabolic mirror of 100 mm focal length. Standard calculation of the beam parameters results in a Gaussian intensity profile with a minimum waist of w 0 = 4.25 μm (5 μm FWHM) and a Rayleigh length of l R  70 μm. The resulting peak power densities which are reachable in the focal region range up to ∼10 17 W cm −2 . Due to the steep laser intensity dependence of the ionization probability, ions will be efficiently produced within a few μm around the laser beam axis and within the confocal parameter l 2 R = 140 μm along the laser beam. For the subsequent calculations, we assume a laser intensity profile with ρ and l being the radial and axial distance from the laser beam center, respectively. The number of ions generated at a position l , r ( ) is then given by Using the density of n = 10 11 atoms cm −3 as given above, we obtain a total number of ions generated per laser pulse as depicted in figure 3. Two points are immediately evident, namely (i) the ionization conditions can in principle be tuned from delivering single ion pulses to pulses containing up to>1000 ions; and (ii) charge states up to q = 5 can be obtained with the available laser power (see figure 2).
At low laser intensity, singly charged ions are predominantly formed with a spatial distribution depicted in figure 4. The data were calculated for a central laser intensity I 0 = 2×10 14 W cm −2 , i.e. before the ionization probability starts to saturate in the center of the laser beam. The upper panel shows that the ion distribution is narrower than the laser beam profile due to the strongly nonlinear intensity dependence of the photoionization probability. Using a threshold value of p i  1%, one obtains an elliptic effective ionization volume with At larger laser intensity, the single ionization probability becomes saturated in the center of the laser beam. At the same time, multiple charge states start to be produced in this region. It should be noted that ions of different charge states are not spatially distributed in the same way. As an example, figure 5 shows charge state selected cross sections of the ionization probability distribution generated at I 0 = 1.2×10 15 W cm −2 . For singly charged ions, the effective ionization volume increases with increasing laser intensity and starts to become depleted in the center of the laser beam once multiple ionization becomes efficient. Higher charge state ions are  produced in this region, while more singly charged ions are produced in the wings of the laser beam profile. These features become important when discussing the space charge broadening of the generated ion pulses (see below).
An important point concerns the role of the photoelectrons generated during the photoionization process. On one hand, the leaving electron will transfer rebound momentum to the ion, which may alter the ion's starting velocity. For a photoelectron excess energy of 1 eV, this 'recoil kick' would impart a velocity of about 8 m s −1 to an Argon ion, which, converted to energy or temperature, would correspond to 1.3×10 −5 eV or 0.15 K, respectively. If this velocity was directed along the ion extraction field, it would significantly interfere with the geometric cooling effect. Therefore, the ionization laser will be linearly polarized in the direction parallel to the propagation of the supersonic beam. In that case, the momentum imparted to the ion will be directed perpendicular to the ion extraction, so that the starting velocity distribution along the extraction direction remains unaffected.
On the other hand, the photoelectron will be accelerated in the ion extraction field and may in principle generate further ions via electron impact ionization of gas atoms. These ions would be produced at different times and larger distances from the beam center and could therefore broaden the temporal width of the generated ion pulse. Using a typical electron impact ionization cross section [35] of the order of 10 −16 cm 2 along with a gas density of 10 11 cm −3 , the probability for the production of such a 'secondary' ion during an electron flight path of 1 mm towards the nearest electrode amounts to 10 −6 , rendering this effect negligible. If the photoelectron hits an electrode surface, electron stimulated desorption may in principle also produce additional ions of the respective electrode material, which will also be accelerated towards the target in the ion extraction field. These secondary ions would, however, arrive at the target surface at significantly different flight times, so that any effects induced by their impact can safely be separated in time from those induced by the ultrashort argon ion pulses generated by photoionization.

Ion extraction
The generated photoions are extracted from the ionization region using an ion-optical bunching configuration sketched in figure 6(a). The setup consists of two (nearly) homogenous electric fields, which are tuned to provide first order flight time focusing conditions at the position of the target surface. A straightforward calculation yields that these conditions are achieved at [36]  denote the corresponding electrical field strengths. The quantity 0 f denotes the starting potential of the ions within the first gap. Since the ion extraction axis is perpendicular to the laser propagation direction, the extension of the effective ionization volume is restricted to<10 μm in that direction by the tight laser focus (see figure 4). Therefore, first order time focusing is sufficient as will be shown below. The intermediate electrode E 2 contains a tapered hole of 80 μm diameter which acts as an aperture to determine the extension of the probed ion extraction volume in the plane perpendicular to the ion extraction axis. This geometrical restriction is important in the direction along the laser beam propagation, where ions are being generated along a range of about 100 μm (see figure 4). The target itself acts as a third electrode E 3 and must therefore be extremely flat and electrically conducting. The implication of these requirements will be discussed below. In order to keep the total ion flight time as short as possible, the entire bunching geometry must be miniaturized as much as possible, while still sustaining the necessary electric fields. The targeted impact energy of singly charged ions impinging onto the sample surface is E 0 = 5 keV. Since the laser beam must enter the buncher from the side, sufficient space must be provided in the first gap in order to prevent intensive laser light hitting the electrode surfaces. Therefore, the laser focus must be placed in the center of the first electrode gap, yielding U 2. 0 1 f = Moreover, flight time focusing according to equation (5) requires that most of the ion acceleration is done in the first gap. As a consequence of these considerations, the distance d 1 was chosen as 2 mm, and d d 2 1 = was chosen in order to provide symmetric time focusing conditions. For this geometry, application of equation (5) yields Together with the requirement U U E 2 , 1 2 0 + = this is fulfilled for U 0, 2 = so that practically the entire ion acceleration is performed in the first gap.
From our practical experience, we know that electric field strengths up to about 10 kV mm −1 can be maintained under ultrahigh vacuum conditions. We therefore chose to conservatively restrict the field strength to values5 kV mm −1 . The resulting principal layout of the buncher is shown in figure 6 (bottom panel). The intermediate and target electrodes E 2 and E 3 form a unit and are spaced by a machined ceramics. Since the two electrodes generating the first acceleration field must be exactly parallel, the top electrode E 1 is mounted on a separate nanopositioning stage equipped with three piezo motors. This way, it is possible to finely adjust this electrode in order to tune the time-focusing conditions while the ion source is running. At the same time, the electrical insulation of this electrode is well separated from the ceramics, thereby reducing the risk of arcing. All three electrodes must be straight and flat, so that we envision to use polished and highly doped silicon wafers for that purpose. The 80 μm diameter hole can be drilled to high precision using a focused laser beam.

Simulations
In order to determine the flight times of ions generated at different starting positions in the photoion cloud, detailed trajectory simulations were performed. The procedure to calculate the flight time distribution of ions impinging onto the target surface was as follows. First, the electric field configuration was calculated with sufficiently high lateral resolution using two different software packages, namely SIMION  in the SIMION code or, more importantly, the Generalized Particle Tracer (GPT, version 3.10) [39]. The GPT code was needed in order to include the influence of space charge broadening, which is not tractable in SIMION with sufficient accuracy. The starting conditions for all trajectory simulations were defined by statistically setting up clouds of different ions, which were generated according to the photoionization profiles discussed in section 2.2. Single ion pulses were calculated first, where only one singly charged ion was positioned at a location r x y z , , 0 =  ( )around the nominal beam center located in the middle of the first gap (z=0) and in the center of the intermediate aperture (x=y=0). The coordinate system used throughout the trajectory simulations is set up as follows. The ions are accelerated along the z-axis towards negative values of z. In the plane perpendicular to the ion extraction, the coordinate r x y 2 2 = + describes the distance from the ion optical axis, where the xand y-coordinates refer to the directions along and perpendicular to the supersonic gas jet, respectively. The ydirection is therefore identical with the propagation direction of the ionization laser beam. For each ion, the vector r 0  was statistically chosen according to the ionization probability distribution displayed in figure 4, and all ions were assumed to start with a velocity v 0  which was statistically chosen according to [40] f where m is the particle mass (40 amu for Ar), k B is the Boltzmann constant, n 0 is the gas atom density, u is the The SIMION code solves the Laplace equation using a finite differences method with the electrodes acting as boundary conditions. For most of the present calculations, a 2D configuration using cylindrical symmetry around the ion extraction axis was used, and the space between the electrodes was discretized in steps of 1 μm in directions of r and z, respectively. In some of the calculations, the planar version was used with 1 μm steps in yand z-directions. In all cases, it was verified that a further reduction of the step size did not change the calculated results. The iterative refinement of the potential field was performed with a goal accuracy of 0.5% ('convergence objective' of 5 mV for an electrode potential of 1 V), and all remaining parameters regarding the field calculation were set to the SIMION default values. In some of the simulations, the ion trajectory calculation was performed using the particle tracker included in the SIMION code. In these calculations, the 'Tquality' parameter was increased until a further increase did not change the calculated ion flight times by more than a femtosecond, resulting in a chosen value of 50. In all calculations, the electrodes were extended all the way to the boundary of the simulation box. In that case, SIMION assumes the electrodes to extend to infinity. The accuracy of the flight time calculation was examined by simulating the field and ion trajectories within an ideal parallel plate capacitor, where the hole in the central electrode was covered with a virtual grid. Comparison of the results with those calculated analytically for this configuration reveals that the SIMION flight time calculation is indeed accurate within less than 10 fs.

CPO field calculation
An inherent problem in the SIMION code is the fact that the boundary conditions for the electric field outside the buncher cannot be properly accounted for. As described above, the SIMION simulations were performed by assuming an infinite lateral extension of the buncher electrodes. In reality, however, the electrodes have a finite diameter and are moreover surrounded by a number of grounded mounts and walls. In order to appropriately account for the role of these influences, one is forced to extend the simulated volume such as to include all outside electrodes. If such electrodes are placed far away from the volume of interest, the SIMION code must therefore be run with full precision throughout a very large space, which is impractical if femtosecond flight time resolution is required. The CPO code, on the other hand, invokes a different strategy to solve the Laplace equation. In the boundary element method applied here, the field is calculated by placing an electrical charge density distribution onto all electrodes. For that purpose, each electrode surface is divided into segments containing a suitable charge, and the Coulomb field of all segments is superimposed to calculate the field at any desired position r .
 The charge q i on each segment (i) is determined by solving the set of linear equations q r , where the values of j f are known from the respective electrode potential. Once the surface charge distribution is fixed, the electric field can be calculated only for the volume of interest, and the calculation can therefore be performed with the desired high precision without overloading the computer memory. This way, the influence of boundary conditions arising from distant electrodes can be studied. In the present case, this was done by enclosing the buncher electrode setup in a grounded cylinder of varying radius and height. The CPO calculations were therefore performed using the same 2D cylindrical symmetry around the ion extraction axis as used in SIMION.
Simulations were performed for different electrode diameter (4 versus 8 mm) and distance between the electrodes (2 versus 4 mm). The electrodes were split into segments using the standard CPO procedures, where the parameters governing the segment density were selected such as to provide a finer density with decreasing distance to the volume of interest, the latter being defined by a radius of 500 μm around the ion extraction axis and therefore sufficient to enclose all relevant ion trajectories. More details of the resulting segment distribution are provided in the supporting material, available online at stacks.iop.org/NJP/21/053017/mmedia. The potential distribution was calculated for a regular grid with 1 μm step size both in r-and z-directions within the volume of interest, with the 'inaccuracy' parameter of the CPO code set to 10 −5 .
When performing ion trajectory simulations using the resulting field, slight discontinuities of the ion flight time were observed which were traced to corresponding discontinuities in the potential energy surface calculated by CPO. These artifacts arise from the rather complicated switching routines between different ways in which the field generated by a specific segment is calculated, which are based on relative distances and orientations between the point of interest and the respective segment. The resulting potential jumps at the switching points were found to be of the order of 1 mV. Since these discontinuities influence the calculated ion flight times on the sub-picosecond level, they were smoothed by iteratively averaging the potential in all cells (i, j) within the volume of interest according to Note that this describes a discrete solution of the Laplace equation in cylindrical symmetry, where the cell indices i and j code the radial distance r and the axial coordinate z, respectively. The equation system (8) was iteratively solved using the CPO-calculated potential surface around and the electrode surfaces within the volume of interest as boundary conditions, until the deviation between subsequent iteration steps fell below 10 −8 V for an electrode potential of 1 V. With the highest electrode potential set to 10 keV, this ensures the potential of all relevant grid points to be accurate within less than 0.1 mV.
The accuracy of the resulting field configuration was investigated by performing calculations both with SIMION and CPO for the same electrode configuration. The results of such a comparison are presented in the supporting material. In the relevant region of interest, both calculations are found to agree with sufficient accuracy to deliver practically the same flight time dispersion.

Trajectory simulations
Once the field configuration was determined, ion trajectories were calculated using the GPT software. The electric field distribution calculated by SIMION or CPO along with the starting position and velocity of the ions was introduced into the GPT code via the respective input files. In addition, the coordinates of the intermediate electrode surface were entered via the 'Rmax' command in order to ensure the stopping of ions hitting the intermediate electrode and exclude them from the calculated flight time distribution. The trajectory calculation was performed for all ions simultaneously, and snapshots of the momentary particle positions were recorded in steps of 1 ns.
The influence of space charge was included by superimposing the Coulomb field of each individual ion to the external electric field. For that purpose, the 'Spacecharge3D' function of GPT was employed, where the coordinates of all particles at a specific time step are used to calculate the force on each particle for the next time step. This way, it is possible to treat the propagation of particle clouds consisting of ions in different charge states, such as shown in figure 5. This possibility is important due to the charge state dependent inhomogeneity of the initial ion distribution. For the situation depicted in figure 5, for instance, the repulsive interaction between the central doubly charged ions is superimposed to the space charge generated by the singly charged ions located in the outer region of the cloud. It should be noted, however, that doubly and singly charged ions quickly separate during the extraction/acceleration process, so that this effect is only operative during the first few picoseconds of their flight.
The trajectories of all ions were followed until the ion hits an electrode surface. Only those ions reaching the target electrode E 3 were counted as part of the generated ion pulse, and their flight time as well as their impact location and energy were recorded using the 'Screen' command of the GPT code. The flight time of each individual ion was calculated to femtosecond accuracy and histogrammed for all simulated ion trajectories. The FWHM of the resulting flight time distribution is then interpreted as the temporal width of the generated ion pulse.

Flight time optimization
All calculations presented in this paper were performed with the potential 1 f of the first electrode E 1 set to +10 kV and the target electrode E 3 kept at 3 f =0 (ground potential). The potential 2 f of the intermediate electrode E 2 was adjusted to deliver optimum flight time focusing for ions starting around the beam center position r 0, 0, 0 0 =  ( ) located in the center of the first gap and the extraction aperture. The procedure used to find these conditions was as follows. First, the electric field distribution was calculated for each electrode separately, with this electrode being set to a potential of 1 V and all other electrodes set to zero potential. Utilizing the linearity of the Laplace equation, the field was then calculated by superimposing the electrode fields multiplied by the actual electrode potential. This way, the high resolution field calculation only needs to be performed once, allowing a fast switching of electrode potentials. For a particular setting of , 2 f a set of trajectories was then run for ions starting with zero velocity on the z-axis (i.e. the extraction axis) at equidistant values of z. Figure 7 shows the resulting z-dependence of the flight time calculated for different values of . Calculating the same flight times using the CPO-generated field, we find an offset by 0.35 ps as shown in the supporting information. Apart from this shift, however, the two curves look identical, indicating that the offset is irrelevant for the shape of the calculated flight time distributions. If the buncher electrode setup is enclosed in a grounded cylinder, the optimum value of 2 f is found to increase, with the shift becoming stronger with decreasing cylinder radius. For the example shown in the supporting material, an optimum value of 2 f = 1080 V is found. If the electrode potential is set to this value, the flight time dependence on the starting coordinate z is identical to that displayed by the blue curve in figure 7. Again, it is in principle possible to bring the optimum value of 2 f closer to zero by either shifting the target surface towards the intermediate electrode or by increasing the gap between the first and intermediate electrodes.
As a consequence, we conclude that the flight time focusing conditions described by the curves depicted in figure 7 can be obtained even in the presence of realistic outside electrodes.

Results
The expected performance of the short-pulse ion source can be deduced from detailed trajectory simulations of the laser-generated photoions. The results of such simulations are presented in this chapter, which is organized as follows. First, we present the arrival time and impact energy distribution which is calculated for single ion pulses, where the flight time distribution of the ions is not influenced by space charge. As shown in section 2.2, pulses of this kind can be generated at low ionization laser intensity up to a peak power density of about 2×10 14 W cm −2 . They represent the ultimate limit of the time resolution which should be achievable with the ion source concept presented here. The following subsection then deals with possible deviations from the assumed ideal conditions regarding the geometry and potential of the buncher electrodes and their influence on the achievable time resolution. If more than one ion is generated in a single pulse, the flight time distribution will be influenced by the Coulomb interaction between the ions during their trajectory. The magnitude of this space charge broadening is estimated in the next subsection. In order arrive at a realistic estimate, we generate ion clouds containing all ions which are produced at a specific laser intensity, which is chosen to deliver the desired number of ions in a specific charge state for a given neutral density.

Single ion pulses
The flight time distributions calculated for single Ar + ions generated by an ultrashort laser pulse at t=0 are shown in figure 8. The calculations were performed with SIMION, but those calculated using CPO and GPT exhibit the same behavior. Technically, the trajectories of 10000 single ions were followed, with the starting positions r 0  of the individual ions being statistically chosen according to the ionization probability function depicted in figure 4. The starting velocity v 0  was statistically chosen according to equation (7) using u = 568 m s −1 and T =  0.54 K. The parallel beam velocity u results in a slight shift of the impinging ion beam with respect to the 80 μm hole in the intermediate aperture, but does not influence the relative timing of the individual ion arrivals. According to the considerations in section 2.1, two different values were inserted for the radial temperature T , which represents the particles' starting temperature with respect to their motion along the extraction direction relevant here.
The first and most important observation is that ion pulses of sub-picosecond duration can in principle be obtained with the concept presented here, provided the radial temperature is restricted to 10 mK or below. As shown in section 2.1, this is achievable via the geometric cooling effect. The arrival time distribution calculated under these conditions exhibits a most probable flight time of 26.6148 ns with a FWHM of 0.45 ps. The comparison in figure 8 reveals that even the reduced temperature of several hundred mK achievable in a supersonic expansion would still limit the achievable time resolution to about 2 ps, and the symmetric distribution calculated for T = 540 mK is apparently determined by the thermal spread of the starting ion velocity. At T = 10 mK, on the other hand, the asymmetric shape of the distribution indicates that the influence of the starting velocity is small. In this case, the distribution is more strongly influenced by the local spread of the ions' starting coordinates in combination with the flight time focusing conditions of the applied buncher setup. This is demonstrated in figure 9, which shows the influence of the starting conditions in more detail. The left panel visualizes the relative flight time of the ions versus their impact location on the target surface. Since the starting conditions of all ions were centered around the point x y , 0,0, 0 0 = ( ) ( ) the shift of the impact zone by about 20 μm in x-direction results from the collective particle velocity u in the supersonic beam. The image shows that the major part of the total flight time dispersion arises from the spread of ion starting coordinates in y-direction i.e. the direction along the laser beam, where the ionization volume extends over ∼100 μm. In the upper right panel, the influence of different starting coordinates is switched off by letting all ions start at x y , 0,0 0 0 = ( ) ( )in the center of the ion beam axis. The resulting distribution reveals the role of the assumed starting temperature of 10 mK, which apparently limits the time resolution to about 300 fs (FWHM). In the lower right panel, the influence of the starting temperature is switched off by setting the starting velocity of all ions to zero. The resulting distribution now reveals the influence of the time focusing conditions achieved in the applied bunching geometry. Figure 10 shows the impact energy distribution of the ions. It is seen that the energy spread E D around the nominal impact energy of about 5 keV is small (∼10 eV). This is different from 'normal' bunching conditions applied in commercial ion sources, where the entire buncher gap is usually filled with ions, thereby producing time focusing at the expense of rather large energy spreads. The reason for the quasi-monoenergetic ion pulses generated here is given by the small spatial extension of the ionization region particularly in the direction of ion extraction.

Experimental influences
An experimental realization of the ion source concept presented here will inevitably introduce imperfections leading to deviations from the ideal situation considered in the preceding subsection, which assumes ideally flat and exactly aligned electrodes at exact potentials. In this section, we estimate the influence of such imperfections on the achievable time resolution. Possible sources of flight time broadening are either the roughness or a relative tilt of the electrode surfaces limiting the first acceleration gap. In order to investigate their influence on the time resolution, simulations were performed with either steps or a tilt introduced to the surface of the first electrode. In that context, the most critical direction is the laser beam propagation direction, since the ionization volume is extended to about 100 μm in that direction. Therefore, the steps as well as the tilt are assumed to occur in that direction. The result is depicted in figure 11(a), which shows the flight time of an ion as a function of its starting coordinate along the laser beam (y-) axis. Different non-ideal surface conditions of the first buncher electrode E 1 are considered, namely (i) a single step of variable height in the center of the electrode and (ii) random steps of 10 μm height statistically distributed over the extension of the ion extraction volume. While the first is intended to mimic a macroscopic tilt of the  electrode surface, the latter illustrates a worst case scenario for the influence of surface roughness (which is usually smaller in amplitude). It is seen that the flight time dispersion calculated for the second case is practically identical to that calculated for the ideally flat surface, indicating that electrode surface roughness, as mimicked by the random step distribution, does not seem to influence the time focusing conditions very much. In contrast, the calculated flight time dispersion reacts very sensitively to a single step distorting the electric field exactly where it is relevant, i.e. in the center of the ionization region. In principle, such a step appears similar to a macroscopic tilt of the entire electrode. The flight time dispersion resulting from such a scenario is displayed in figure 11(b). It is seen that a tilt influences the flight time dispersion in quite a dramatic way and may lead to a destruction of the sub-picosecond time resolution if the tilt angle exceeds 0.5°. As a consequence, we conclude that the electrode angle must be controllable from outside the vacuum chamber, and the electrode must be aligned on-line in order to achieve optimum time resolution while the ion source is running.
Apart from the geometric effects, another possible source of increased time dispersion may arise from instabilities of the electrode potentials. The dependence of the relative ion flight time on small deviations of the electrode potentials is shown in figure 12, which reveals that deviations of the order of ±0.1 V can be tolerated, before the voltage stability starts to limit the achievable time resolution. In connection with a nominal electrode potential of 10 kV, this requires a stability of the order of 10 −5 for the high voltage power supply which appears feasible.

Space charge influence
The influence of space charge broadening must be investigated under realistic conditions regarding the generated ion cloud. At low laser intensity, only singly charged ions are produced, and the resulting ion cloud can be generated by statistically distributing a given number of ions over the ionization volume according to the ionization probability distribution depicted in figure 4. The flight time distributions calculated for such a situation are displayed in figure 13 for pulses containing different numbers of ions. More specifically, clouds with average numbers of 10, 20 and 100 ions are used, corresponding to central laser intensities of I 0 = 2.2, 2.5 and 3.6×10 14 W cm −2 , respectively. The resulting average number of ions actually reaching the target surface is indicated in the figure. For reference, the single ion distribution corresponding to I 0 = 1.5×10 14 W cm −2 is included in panel (a). It is seen that space charge broadening is negligible for a pulse containing 10 ions, but starts to become significant for pulses containing 100 and more ions. Detailed inspection shows that sub-ps time resolution can be maintained up to about 25 Ar + ions per pulse.
At larger laser intensities, doubly charged ions start to be produced in the center of the ion cloud. For I 0 = 4.8×10 14 W cm −2 , one Ar 2+ is on average produced per pulse along with ∼200 singly charged ions. If I 0 is increased to 6.3×10 14 W cm −2 , about 10 Ar 2+ ions are produced along with ∼400 Ar + . It is of note that the singly charged ions are located in the outer regions of the cloud, so that not all of the generated ions actually reach the target surface. While all of the doubly charged ions contribute to the target ion pulse, only about 50% of the generated singly charged ions are transmitted through the intermediate aperture. The resulting flight time distribution calculated for the doubly charged ions is shown in figure 14. For comparison, hypothetical distributions which are calculated by omitting the surrounding singly charged ion cloud are included in the figure. It is obvious that the influence of the space charge generated by the accompanying Ar + ions is small, a finding which becomes understandable when considering the rapid separation of singly and doubly charged ions during the acceleration process.
Another point worth noting is the fact that multiply charged ions impinge with larger kinetic energy and are registered at an accordingly smaller total flight time than the singly charged ions contained in the same pulse. In our system, the impact energy of doubly charged ions is about 10 keV, and the arrival time difference between the Ar 2+ and Ar + sub-pulses amounts to about 8 ns. This difference is sufficiently large to completely separate the dynamics induced by the impact of doubly and singly charged ions, thereby in principle allowing to study both processes at the same time. While the time resolution expected for multiply charged ions is comparable to that achievable by pulses containing the same number of singly charged ions, the flight time distribution of the accompanying Ar + sub-pulse will of course be significantly broadened by space charge. As a consequence, we conclude that pump-probe experiments using singly and multiply charged ions with energies in the keV range appear feasible with sub-picosecond time resolution.

Conclusion
The simulations performed in this work clearly demonstrate that the generation of ultrashort ion pulses with sub-picosecond duration is in principle possible even for heavy ions at relatively low (keV) kinetic energies. The described concept allows to generate pulses containing up to several 10 ions of 5 keV energy with a FWHM temporal width of a few hundred femtoseconds at an energy width of about 10 eV. An experimental realization of the described ion source concept appears principally feasible and is currently under development in our lab. The generation of femtosecond keV ion pulses will open new possibilities to investigate ultrafast relaxation processes which are triggered by an energetic particle hitting a solid surface. Since the ion pulses are generated by means of a laser-based technique, pump-probe experiments using a tuned delay between ion pump and laserbased probe techniques such as those described in the introduction via optical delay lines become possible. These techniques are well established for experiments where the sample is excited by a laser pulse. Using the concept presented here, similar experiments now appear feasible also for the investigation of ion impactinduced non equilibrium states of matter, thereby for the first time providing an experimental access to the ultrafast lattice and electron dynamics following a particle impact, which so far are only accessible via model calculations. The limitation to pulses containing only a few ions is not a problem, since many experiments investigating, for instance, ion induced sputtering processes, are routinely being performed on an event-byevent basis, where the effect of a single ion impinging onto a solid surface is studied [41][42][43][44]. The results can then provide a deepened understanding of the extremely localized non-equilibrium dynamics generated by a particle impact and allow an experimental assessment of the fundamental assumptions underlying the respective model calculations.