3D simulations of surface harmonic generation with few-cycle laser pulses

The mechanism of harmonic generation in the interaction of short laser pulses with solid targets holds promise for the production of intense attosecond pulses. Using the three-dimensional code , we have performed simulations pertaining to an experimentally realizable parameter range by high-power laser systems to become available in the near future. The emphasis of the investigation is on the coherent nature of the emission. We studied the influence of the plasma scale length on the harmonic efficiency, angular distribution and the focusability using a post processing scheme in which the far-field of the emission is calculated. It is found that the presence of an extended density profile reduces significantly the transverse coherence length of the emitted extreme ultraviolet (XUV) light. The different stages of the interaction for two particular cases can be followed with the help of movies.


Introduction
The early studies of harmonic generation from solid targets using tabletop laser systems [1]- [3] revealed some fundamental properties regarding the spectrum and its dependence on laser intensity and plasma conditions. More recently there has been renewed interest in this harmonic generation process with the goal to scrutinize its origin and study in detail the spectral but also temporal characteristics of the emitted radiation. The reason is the envisaged potential as a method for the efficient generation of phase-locked harmonics leading to attosecond light bunching [4]- [9]. The motivation is the rapid technological advancements in laser technology that promises in the near future tabletop lasers delivering several tens of TW (10 12 W) power with kHz repetition rate [10]- [14] and thus the availability of small scale easily accessible high intensity set-ups. The functional production of intense attosecond pulses in the XUV spectral region would constitute a breakthrough in the quest towards real time studies with atomic resolution of fast dynamical processes. In contrast to the attosecond sources currently available using the process of harmonic generation in rarefied gases, the new mechanism is expected to provide attosecond pulses with enough intensity to enable XUV-pump XUV-probe types of investigations.
One of the most significant and recent findings is associated with the roll-off of the harmonic spectrum produced by the laser-plasma interaction at high intensities and its cutoff. It has been shown theoretically [5] and confirmed experimentally [15] that it exhibits a power-law dependence of the form I XUV ∝ ω −q with q ≈ 2.5-2.7. Moreover its maximum frequency (cutoff) increases with the laser intensity I L according to ω co ∝ I 3/2 L [16], but unlike the case of harmonics from an atomic medium there is no limitation on the maximum intensity that can be used (no saturation intensity).
The implications of this finding are substantial as they predict a high laser to XUV light conversion efficiency. For a laser intensity in the range of I L 10 20 W cm −2 , the estimated conversion efficiency η XUV in the 20-70 eV range is a few percent, while η XUV ≈ 10 −3 and ≈ 10 −4 in the 80-200 eV and 400-1000 eV range correspondingly [7]. Furthermore due to the coherent nature of the generation process the individual harmonics are phase-locked, giving the possibility to generate diffraction limited attosecond (1 as = 10 −18 s) pulses extending to the keV spectral range, which can be focused to higher intensities than the incident beam. This, in combination with the availability of compact laser systems delivering exceptionally high intensities, accentuates the advantage of the new harmonic generation mechanism as a method to produce intense attosecond pulses and explains the interest associated with it.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
To date, there are a number of reports on experimental investigations dealing mainly with the origin of the generation mechanism based on the frequency spectrum produced by it. The most recent experimental results have been obtained using high contrast laser pulses [15], [17]- [20]. Two important findings have resulted from these investigations. (i) In the relativistic limit I L 10 18 Wcm −2 , high efficiencies (η XUV 10 −6 ) at wavelengths <4 nm are indeed possible and in addition, the spectrum exhibits a slow decay of harmonic efficiency to high orders that follows very closely the theoretically predicted power law. (ii) The long standing dispute over a density dependent cutoff of the power spectrum was settled. In independent experiments [17]- [19] it was demonstrated that at weakly relativistic intensities I L 10 19 W cm −2 there is indeed a density dependent cutoff because in this intensity range the dominant harmonic generation mechanism is due to plasma oscillations within a density gradient that lead to the emission of harmonics up to the plasma frequency.
Due to the coherence of the generation process, the harmonics are phase locked, giving the possibility to generate attosecond pulses when proper filters are applied. No measurement in the time domain has been performed so far although efforts towards characterization of the temporal behaviour of a filtered part of the spectrum are currently in progress. Theoretically, the predictions concerning the suitability of the harmonic spectrum generated at the plasma vacuum interface to produce a train of or single attosecond pulses are based on one-dimensional (1D) analytical models or 1D particle-in-cell (PIC) simulations. A simplified but very helpful generation model is that of the oscillating mirror [7,21] where the electrons at the plasma vacuum interface execute forced oscillations near the edge of an immobile step-like ion background driven by the ponderomotive force of the incident laser pulse.
The spatial coherence of the emitted harmonics cannot be verified in 1D modelling. Especially the transverse variation of the laser intensity and the subsequent deformation of the vacuum-plasma interface may result in dephasing of the oscillating dipoles. As no experimental results exist to confirm the spatial coherence necessary for the formation of attosecond pulses, the only possibility to investigate this issue is to use 3D-PIC simulations.
Here, we present 3D simulations of high harmonic generation on a plasma-vacuum interface using the 3D-PIC code ILLUMINATION [22]. The three-dimensionality of the code allows for the finite transverse extension of the interaction region and for the resulting deformation. This enables us to address not only the issue of coherence under these circumstances, but also the effect of the surface deformation on the angular divergence of the emission. The radiation emanating from the simulation box is post-processed, and the far field thus obtained furnishes a clear demonstration of the dependence of the angular spread on the plasma scale length in front of the plasma surface, thereby explaining some older observations [23].

3D-PIC simulations with the code ILLUMINATION
The simulations are performed with the 3D-PIC code ILLUMINATION. The typical simulation box is 500 × 500 × 1500 cells in x-, y-and z-direction respectively and has a volume of 12 × 12 × 18 µm 3 . The target is placed at z ≈10 µm down the propagation direction of the laser pulse (see figure 1). Each target cell contains one macro-particle corresponding to 4 million particles. With the typical grid size of ≈10 nm the highest resolvable wavelength is the 20th harmonic of an 800 nm laser pulse. The target is assumed to be fully ionized with immobile ions as a neutralizing background. We have checked that allowing for ion motion gives almost the same results. Throughout the paper the laser pulse is linearly polarized in the x-direction, propagates in the z-direction and has a central wavelength of λ 0 = 800 nm. The normalized vector potential of the incident laser pulse, a = eA/mc 2 , depends on radius and time in the form Here e and m are charge and rest mass of electrons, and c is the velocity of light. For the reference pulse, we choose a focal radius of w 0 = 5 µm and a full-width-at-halfmaximum (FWHM) duration of the intensity of τ 0 = 5 fs. This implies a Rayleigh length of z R = πw 2 0 /λ 0 = 98 µm. Further the intensity, the peak power, and the total pulse energy can be calculated as I = 2.2 × 10 18 × a 2 0 W cm −2 , P = πw 2 0 I/2 and ε = 1.67Iτ 0 w 2 0 , respectively. The vacuum-plasma interface is modelled with an exponentially decaying density profile n = n 0 exp(z/L p ), where L p is the plasma scale length and n 0 the plasma electron density. It is chosen as n 0 = 1 × 10 22 cm −3 , corresponding to ≈6 times the critical density for λ 0 = 800 nm 5 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT light. In the simulations presented here, all the parameters are held fixed with the intensity at I = 5.5 × 10 19 W cm −2 corresponding to a 0 = 5. Only the plasma scale length L p in front of the target is varied between 0 and 5 µm.

Influence of the plasma scale length
An important experimental issue in the generation of harmonics from solid targets is the influence of the plasma scale length resulting from the ionization and heating of the target by the presence of the inevitable pedestal and/or prepulse in the laser pulse. Even early studies showed that the existence of anything else but an abrupt density profile in front of the target surface has detrimental consequences to the efficiency of the harmonic production [3].
A successful but cumbersome way to circumvent this problem is with the help of plasma mirrors as described in [15], but here it was achieved by losing ≈70% of the energy at the mirrors. It is therefore rather important to know what density scale length can be tolerated before the efficiency and the coherence reach unacceptable levels.
Using the 3D capabilities of the ILLUMINATION code, we studied the influence of the scale length of the expanding plasma L p on the overall development of the 5 fs pulse interaction with the overdense plasma. The emphasis is primarily on the efficiency of the high harmonic generation, but more importantly on the spatial coherence of the back reflected XUV radiation. Figure 1 shows a snapshot after the interaction of the laser pulse with a step-like density profile (L p = 0). At this time, the pulse has been reflected and propagates to the left. The upper two panels present the electron density in the target and the intensity of the reflected laser pulse. To obtain short XUV pulses, we have applied a virtual pass-band filter to the back-reflected light, allowing only the 5th to 10th harmonic to pass. Although higher harmonics are present in the simulations, due to the limited spatial resolution, they propagate with substantial numerical dispersion and their inclusion would give unphysical results. The lower two frames of figure 1 show the intensity of the filtered light with spectral range comprising the 5th to 10th harmonic in the x-z-plane (right) and the y-z-plane (left). As it can be seen in the time-domain, the result of the filtering is a pulse train with the pulses separated by half a period of the fundamental laser wave and having a duration of a fraction of this. For an 800 nm laser pulse, this corresponds to individual pulse durations of few 100 as. The extent of the train and the wavefront curvature appear to be slightly different in the two planes. A line-out along the laser axis is shown in figure 2. Beside the instantaneous intensity of the laser field, the pulse train of the XUV field is shown having a carrier frequency close to the 5th harmonic. The existence of an attosecond pulse train itself provides strong evidence that the emission in the near field and within the selected spectral range is coherent.
The various stages of the interaction can be followed in movie 1: it shows a sequence of frames of the same quantities as in figure 1. The time t = 0 is when the laser pulse is located at about 9 µm in front of the plasma vacuum interface, and it is followed for about 70 fs afterwards. After 20 fs the pulse reaches the vacuum-plasma boundary and starts interacting with the electrons inside the target. The combination of the ponderomotive force of the laser pulse and the restoring electrostatic force of the immobile ions leads to the characteristic anharmonic but periodic oscillation of the surface. This can be seen in the movie during the next 17 fs in which the main part of the laser pulse is reflected. Throughout the entire interaction time the regularly oscillating surface and the steep plasma gradient are maintained. This observation provides a sound justification for the use of the simple oscillating mirror model to gain insight into the mechanism of high harmonic generation at high intensities [7,21]. Because of the short pulse duration the interaction surface stays intact despite the high intensity. Longer pulses would deform the oscillating surface by pushing it into the target. Thus, the reflecting surface will go through a denting phase which is initially likely to affect divergence but not focusability. The timescale for this deformation depends on the target density and laser intensity and evolves on the ion motion timescale (typically in the 100 fs range). At even later times rippling develops, which destroys the integrity of the reflecting surface and therefore the focusability of the reflected light [24] (see also the L p = 5 µm case below). The reflected pulse is ≈10 fs long, it contains ≈90% of the incoming pulse energy and has an intensity in excess of 10 18 W cm −2 . The lower frames in the movie show propagation of the filtered high harmonic pulse in the x-z-plane (right frame) and in the y-z-plane (left frame). During the whole period of reflection high harmonic spikes are generated, separated by half the laser wavelength. According to the oscillating mirror model, the peak of the high harmonic emission coincides with the peak of the laser pulse. Following the pulses as they propagate out of the simulation box suggests that the harmonics are coherent. This will be studied in the subsequent far-field section.
The situation changes dramatically when the plasma scale length is increased to L p = 5 µm. Figure 3 shows the same quantities as figure 1 for this case. At the time of 37 fs in movie 2, the reflection of the main pulse is almost over and the pulse propagates to the left. The incoming laser pulse propagates up to the critical density in the density ramp, i.e. until it reaches n crit = 1.75 × 10 21 cm −3 . In contrast to the step-like profile, where the reflection takes place at a well defined surface, the laser pulse is now partially reflected from different parts of the interaction surface, while still propagating towards higher density. The total thickness of the reflecting volume is ≈5 µm. Further the reflecting surface is no longer a mirror-like plane, but it is strongly deformed. Only 50% of the incoming pulse energy is reflected. This highly complex reflection process can be followed in movie 2. The formation of density bubbles reminiscent to earlier 2D simulation results [24] are clearly seen in the density evolution during the interaction. Due to the curved and dented surface the reflected high harmonic pulse consists of several pulses propagating in different directions. This leads to a strong reduction of the peak intensity, by one order of magnitude compared to the steep gradient. Here the timescale for the growth of this Rayleigh-Taylor-like electron instability is in the fs-range and thus comparable to the laser pulse duration.
An interesting feature appears in the movie after 15 fs: behind the laser pulse at z = 5 µm the first peak of a wake field pattern appears with a peak density of 5 × 10 21 cm −3 . The leading part of the incoming laser pulse is already reflected at ≈7 µm and propagates in the negative z-direction towards this overcritical density spike. This leads to a re-reflection and doppler upshift of the laser pulse. This can by seen by the short high harmonic pulse propagation in the +z-direction (towards the target). This intricate influence of the plasma scale length in  the harmonic emission is revealed only by a 3D simulation and was not observed in previous 1D studies [3]. The differences between 1D and 3D simulations will be discussed in section 5.

The far-field of the reflected light
In order to investigate coherence, focusability, and divergence of the high harmonic pulse one has to know the far-field distribution of the reflected light. Since the fields inside the simulation box are totally known, the far-field can be calculated using standard Kirchhoff diffraction theory [25]. If we denote as V(Q, t) an arbitrary field component at time t and at point Q on a surface S inside the simulation box then the corresponding field component V(P, t) at point P on a plane at distance s is given by: where ∂/∂n means differentiation along the inner normal of the surface S, and the terms enclosed in square brackets are calculated at the retarded time t − s/c. Let us first consider the case of figure 1 with sharp plasma gradient. The intensity spectrum after propagation in vacuum 200 µm off from the target is given in figure 4. It exhibits the characteristic for high harmonic generation from surfaces according to the power law scaling I(ω) ∝ ω −p . Theoretical analysis of the surface harmonic generation mechanism by Baeva et al [16] predicted a value of p = 8/3 for the exponent. This is in excellent agreement with our simulations as can be seen by the fitted curve in figure 4. The cutoff at the 20th harmonic is an artifact due to the limited temporal resolution of the code.  In figure 5(a) we have plotted the 3D iso-intensity contour of the infrared (red) and XUV (blue) pulse resulting from the filtering of the 5th to 10th harmonic. The infrared pulse is given in terms of its instantaneous intensity [E 2 IR (t)], and the XUV pulse in terms of [ E 2 XUV (t) ]. The localization in space of the attosecond pulse train is clearly seen. In figure 5(b) the energy flux is plotted for the fundamental and the high harmonic pulse. The far-field exhibits excellent symmetry of both pulses with considerably narrower angular distribution for the XUV pulse. This clearly indicates the coherent nature of the high harmonic generation and the focusability of the high harmonic pulse. One can compare the calculated width with that of a Gaussian pulse after propagation over a distance z = 200 µm obtained from the expression w(z) = w 0 × 1 + (z/z R ) 2  coherently generated almost over the whole spot size illuminated by the laser pulse on the target with a corresponding transverse coherence length of l XUV ≈ 2w 0 .
The situation is quite in contrast to the case where an exponential density profile with a scale length of L p = 5µm is present. A similar 200 µm propagation of the fields shown in figure 3 and movie 2 yields IR and XUV beams that are strongly divergent. This is depicted in figure 7 where a line-out of both fields along the x-axis is given. The spot size in this case has increased by a factor of eight compared to the step-like density profile. Accordingly, the transverse coherence length in this case is reduced by the same factor, i.e. l XUV (L p = 5λ) ≈ l XUV (L p = 0)/8 ≈ w 0 /4. This is attributed mostly to the strong deformation of the 'soft' interaction layer and the reflection of different parts of the pulse in different directions due to indenting of the plasma surface (see movie 2), resulting not only in increasing the spot size but also leads to a much broader and higher pedestal for the XUV-pulse compared to the sharp gradient case. As can be witnessed in movie 2, the plasma surface at the beginning of the interaction stays intact for a short time but later as the instabilities develop it breaks up into irregular reflection regions. It is reasonable to presume that the sharp peak at the centre of the profile in figure 7 is coherent and stems from the first phase of the interaction, while the broad pedestal is produced during the second phase.

Comparison between 3D and 1D simulation results
Finally we address the question to what extent a 3D description of the high harmonic generation is really essential. For this we have compared the 3D results presented here for the L p = 5µm case with results from a 1D-PIC code. It is found that, for early interaction times, both codes deliver comparable results. However, at later times the 1D description starts to deviate from the 3D one (see figure 8). In 1D the laser pulse pushes the electrons inside the target to form a sharp density peak at a well defined position inside the density ramp, where the reflection and the high harmonic generation (HHG) takes place. During this stage of the interaction the 3D simulation shows a completely different scenario. Because of the transverse forces of the laser pulse, the electrons are pushed out of the focal region on a transverse plane reducing the peak density in front of the laser pulse (at 7 µm in figure 8). Similar to the propagation of an intense laser pulse in an undercritical plasma, the electrons return to their initial position creating a high density peak just behind the laser pulse (at 5 µm in figure 8). While the laser pulse pushes the first peak inside the plasma, it is at the same time also reflected. Therefore the reflection does not take place in a well defined position, but is rather smeared out over several wavelengths which reduces the coherence and efficiency of HHG.

Conclusions and outlook
We have investigated the characteristics of the surface harmonic generation mechanism using a 3D-PIC code. The investigation is focused on few-cycle TW-to-PW (10 15 W) laser pulses in view of the envisaged high intensity tabletop laser systems to become operational in the near future. Based on 3D-PIC simulations, we have emphasized issues of immediate interest to experimentalists. We present two movies of a reference case, thus providing much more detailed information than has been published so far mainly using 1D codes. The simulations presented here show that the plasma scale length of the expanding plasma plays a major role even for 5 fs laser pulses. While for a step-like density profile the reflected XUV pulse is well behaved with excellent characteristics to warrant focusing of the radiation into a diffraction limited spot, the existence of a pre-plasma leads to severe 12 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT deterioration, manifested by rapid increase of the divergence due to the reduction of the transverse coherence length. The simulations give an upper limit for efficient high harmonic generation of L p < 5 µm for sub-10 fs pulses. The demands on the pulse contrast and pre-pulse duration to fulfill this limit are quite high. Assuming for an Al (Z = 13) plasma a typical electron temperature of kT e ≈ 100 eV from a 10 11 -10 12 W cm −2 pre-pulse [26] results in a plasma scale length [27] L p = c s × t = √ kT e /M · (Z + 1) × t < 5 µm only for pre-pulse durations t < 100 ps. Here, c s is the sound velocity and M and Z are the ion mass and charge respectively. From this one can estimate that the pre-pulse must remain below the plasma formation threshold (≈10 12 W cm −2 for metals, >10 13 W cm −2 for dielectrica) at t > 100 ps before the peak of the pulse.
High harmonic generation has been demonstrated in [15] using long (500 fs) high-contrast PW-range pulses. Few-cycle pulses of corresponding contrast and power will improve the quality of HHG and lead to to ultrabright, single attosecond pulses.