Propagation effects in multipass high harmonic generation from plasma surfaces

Multipass high harmonic generation from plasma surfaces is a promising technique to enhance the efficiency of the generation process. In this paper it is shown that there is an optimal distance between two targets where the efficiency is maximized, depending on the laser and plasma parameters. This can be explained by the Gouy phase shift, which leads to the relative phase between the colours being changed with propagation in free space. A simple model is used to mimic the propagation of light from one target to another and to observe this effect in 1D particle-in-cell (PIC) simulations. The results are also verified using 2D PIC simulations.


Introduction
Surface high harmonic generation (SHHG) is a promising route towards intense ultrashort XUV pulses generated in laser-plasma interactions [1][2][3][4][5][6][7][8][9]. Several different techniques have been proposed since to enhance the efficiency of the harmonic generation process. One technique is the mixing of several frequencies in the driving laser pulse, either by including the second harmonic [10] or by using multiple passes [11][12][13][14]. While the former has been demonstrated experimentally [15], the latter still awaits experimental investigation. An important factor, which has not been considered so far, is the distance between two targets in a multipass setup. It is especially important for possible experiments to find the optimal distance, while still being practical from a technical point of view. As already outlined in previous work [12], the waveform of the driving laser pulse is important for the efficiency of the harmonic generation process. But due to the Gouy phase [16][17][18][19], the waveform of a Gaussian laser pulse with multiple colours will change [see also equation (4)], even when propagating in free space [17,20]. Furthermore, the Gouy phase shift is different, depending on the number of spatial dimensions. While this does not seem to be important for experiments in the three-dimensional real world, often two-dimensional particle-in-cell (PIC) simulation are used to support experimental data. The paper is structured as follows: in section 2, the differences between a Gaussian beam in paraxial approximation in two and three dimensions are explained. In section 3, the consequences of the Gouy phase shift for the waveform of a multicoloured Gaussian beam are described. In section 4 the numerical methods are outlined in detail. Finally, in section 5 the results of the numerical study are presented and conclusions for possible future experiments drawn.

Differences of Gaussian pulses propagating in 2D and 3D geometry
The behaviour of a Gaussian pulse is usually described with the slowly varying envelope approximation. For a linearly polarized beam in 3D space propagating in z direction it reads with v(r) the slowly varying part in paraxial approximation v(x, y, z) where ϕ = arctan −z/z r is the Gouy phase and z r = πw 2 0 /λ 0 the Rayleigh length. If not stated explicitly otherwise, in this paper z r always refers to the Rayleigh length associated with the fundamental wavelength of the incident laser pulse. The focus position is in this coordinate system at z = 0. Naturally, this equation is derived for a 3D geometry, where the pulse propagates in z and has transverse components in x and y. Plasma physics simulations, especially PIC simulations of SHHG are often run in 2D for computational efficiency. It is therefore important to note that the behaviour of Gaussian pulses is different in 2D geometry compared to 3D. In particular, the intensity falloff is slower, the Gouy phase is only half the value and the focus is a line focus rather than a point focus. The paraxial equation in 2D for a Gaussian pulse propagating along the x axis reads v(x, y) = A 0 The important consequences are that running a 2D simulation with two targets, the intensity and the Gouy phase on the second target will not match the 3D case and will lead to different results. However, if one takes this effect into account, 2D simulations will still be a valuable tool to demonstrate the principle.

Gouy phase effect on multicolour Gaussian beams
Multicolour Gaussian beams are a superposition of Gaussian beams with different central wavelength λ 0 , such as those generated during SHHG on overdense targets [1-6, 8, 9]. Here, a beam with wavelength λ 0 and its second harmonic λ 0 /2 is modelled, as the contribution of the second harmonic is the dominant one. Of course, higher harmonics contribute to the waveform as well and can be treated analogously. The impact of the Gouy phase on the waveform of the multicolour Gaussian pulse is explained here in this example to show the principle. The length of one cycle (2π in radians) is defined by the wavelength, but the different colours in the Gaussian beam are viewed at the same absolute distance z. This is important when calculating the relative phase between different colours. The complete phase of a Gaussian beam on axis is The term k z − ωt associated with the propagation of the beam increases linear with z with a slope proportional to ω. Therefore, the difference k 1ω z − ωt − (k 2ω z − 2ωt) defines a point of constant phase relation between the 1ω and 2ω beam. As we are interested in the relative phase at an absolute z position for every time t, both propagation terms have to cancel each other. Therefore, to be comparable, one has to decide to measure the phase difference either in radians of the 1ω or 2ω frequency.
Here, it is given in terms of the 2ω beam, hence there has to be a factor two in front of ϕ 1ω . The relative phase between the 1ω and 2ω component is thus Δϕ = 2ϕ 1ω − ϕ 2ω or where z R,1ω = πw 2 0 λ −1 L and z R,2ω = πw 2 0 λ L /2 −1 , assuming a common source size w 0 for both colours. The case of deviating source sizes is discussed in the supplementary. For a fixed absolute propagation distance, the Gouy phase shifts for the individual wavelengths will be different and thus the relative phase between them will change while propagating. This effect is already known from broadband few-cycle single colour pulses [21,22]. This will change the waveform dramatically. They are shown in figure 1 for three different cases: in focus (figure 1(b)) and 1z r (measured in terms of the wavelength of the fundamental) before and after the focus (figures 1(a) and (c)). As a result, a symmetrical waveform in focus changes to an asymmetric form out of focus. It is known from experiments with 2-colour laser pulses, that the waveform and relative phase is crucial to the generation efficiency of surface high harmonics [15]. It is therefore a good assumption that this will also play a role in the enhanced harmonics generated in a multipass geometry [12,13]. Here, the relative phase between harmonics generated on the first target is determined by the distance they need to travel until they hit the second target. In addition to the waveform, the second important factor for the generation efficiency is the intensity reduction when propagating further away from the focus. However, as it is shown in section 5, the optimal waveform has a larger impact and outweighs the reduced intensity. We tested this assumption by conducting a wide range of 1D and 2D PIC simulations.

PIC simulations
The simulation setup (figure 2) consists of two targets A and B of an overdense plasma (n = 200n c , where n c is the critical density for the fundamental laser pulse frequency for a quasi-neutral plasma) with a step-like density distribution (no scale length). Target A is placed at the origin (0, 0) of the coordinate system at an angle of θ = 45 • . The second target is placed at position (0, −d) under an angle of θ = −45 • . This alternating target orientation is important to be able to achieve an enhancement at all. If the targets would both be aligned at the same orientation (in this example, both at θ = 45 • ), the harmonic intensity after the second pass would have been about the same as after the first pass and no enhancement would be possible. Only for alternating target orientations the second pass will enhance the harmonics. The reason for this is that one not only has to consider the intensity, but the orientation of the electric field. After reflection from target A, the shape of the pulse is asymmetric with the sharpest peaks in amplitude in the negative half of the cycle. Only if the second target is oriented alternating to the first target, these peaks of the electric field will pull the electrons across the plasma-vacuum boundary and generate the enhanced harmonics. A Gaussian laser pulse with the central wavelength λ 0 is injected at the −x boundary. The laser pulse has a maximum amplitude of a 0 = 10, a full width half maximum (FWHM) duration of τ = 5 cycles and is p-polarized. On target B, the laser pulse is still p-polarized, as the polarization remains the same after the reflection from target A. These laser and plasma parameters are varied in section 5.1. Where not noted otherwise, the laser pulse is focussed at the origin (0, 0), i.e. at the surface of target A. For focussing on target B, the focus of the laser pulse is placed at (d, 0) with target A acting as a mirror. During the interaction with target A, high harmonics are generated in reflection due to the relativistic oscillating mirror [1,2]. After propagating the distance d, the laser pulse, now containing high order harmonics next to its fundamental frequency, interacts with target B creating even more intense XUV pulses [13]. In 1D PIC simulations, the 2D geometry can be emulated by a Lorentz transform to the boosted frame, in which the laser pulse is incident normally on the target [1,2,23]. In this frame, the plasma moves with a velocity of v y = sinθ and the electric field is transformed as E y = E y cosθ. As there is no focus in 1D, the laser pulse is injected at the boundary with the intensity it would have in focus in 2D or 3D. Between the targets, the intensity and the Gouy phase is modelled according to the paraxial approximation (see section 4.2). All 1D PIC simulations, conducted with the code PICWIG1D [24], had a length of 4λ 0 with a resolution of 500 grid cells per wavelength. The plasma slab had a thickness of 0.8λ 0 and there were initially 500 macro-particles per grid cell. All 2D PIC simulations have been performed with the code EPOCH [25]. The 2D PIC simulation for the first pass had a box size of 50μm × 50μm with a resolution of 400 cells per μm.
The laser pulse had a central wavelength of 800 nm and a FWHM spotsize on target of 5μm. The target thickness was 1.0μm. Initially there were 20 macro-particles per grid cell. The pulse was retrieved after the first target and propagated in free space to the distance needed using the method described in [26]. It was then inserted in another EPOCH simulation with a box size of 50μm × 30μm and 400 cells per μm resolution. The target had the same properties as in the first pass, but was rotated by θ = −45 • . The results are shown in terms of the attosecond pulse (AP) generation efficiency η = I AP /I L , where I AP and I L are the maximum intensities of the AP and the incident laser respectively. The AP is obtained by applying a highpass filter which cuts out everything below the 20th harmonic. For better comparison to the ordinary SHHG process with one pass, the axis are further normalized to the efficiency η S of the single pass harmonics.

Modelling the Gouy phase in 1D PIC simulations
Performing full scale 3D PIC simulations, especially of SHHG in double reflection, would require a huge number of central processing unit-hours of computing time even for a single simulation, not to speak of simulating different parameters. Even the required resources for an extensive parameter scan in 2D would have been substantial. As this was not feasible for this study, it was possible to reproduce the effect of the Gouy phase by artificially imposing the appropriate Gouy phase on the pulse after the first pass. For the 1D simulations shown in this paper it is assumed that the source size of the generated harmonics matches the spotsize of the fundamental. While this is an approximation, it does not change the conclusions drawn from the simulations. A detailed discussion of the source size is included in the supplementary material (https://stacks.iop.org/NJP/22/093048/mmedia). The model was verified by imposing a 2D Gouy phase and comparing the results to full-scale 2D PIC simulations, which do not make any assumptions about the harmonic source size. In 1D geometry, propagating pulses in free space show no diffraction and no phase shift. To model the free space propagation in higher dimensions between the two targets in a 1D PIC simulation, both the intensity and the phase had to be adjusted accordingly. The reflected pulse was obtained after the first target and, after a Fourier transform, the Gouy phase was artificially imposed on each spectral component of the reflected electromagnetic fields. Likewise the intensity was reduced according to the paraxial approximation. Afterwards the pulse was injected into another simulation for the interaction with the second target. This method was cross-checked with the model detailed in section 3.

Results
The resulting AP generation efficiencies for different distances are shown in figure 3. The contributions of the intensity falloff due to diffraction and the Gouy phase are shown separately in figures 3(a) and (b). As one would expect, the intensity falloff (a) leads to a decreasing efficiency with larger distances. If one is taking into account the Gouy phase only (b), the efficiency reaches a maximum at about 0.2z r (3D) and 0.5z r (2D). The actual behaviour is the combination of both effects shown in (c). This can be explained as follows: it is already known from 2-colour experiments [15], that the harmonic generation efficiency is only enhanced for a certain phase difference between the 1ω and 2ω component of the laser pulse. It is therefore sensible to assume that such a favourable phase relation also exists for a beam containing high order harmonics of the fundamental frequency. As the Gouy phase evolves differently for every frequency as outlined in section 3, the relative phase shift between the frequencies changes constantly with the propagation of the beam. The maximum efficiency is achieved at the distance, where the waveform is favourable for harmonic generation. This optimal distance between the two targets is not as close as possible, but in the range between 0.2z r and 0.5z r . An enhancement, albeit smaller, can be expected for distances up to 1.0z r . For longer distances z > z r , the effect of the decreased intensity becomes dominant and the efficiency is approaching 0. As the beam propagates differently in 2D and 3D geometry, the harmonic generation efficiencies differ as well. Because the intensity falloff is slower in 2D, the efficiencies remain higher for longer distances. Furthermore, in the 2D case the Gouy phase lies in the range of [−π/2, π/2], which is half the range of the 3D case and shifts the maximum efficiency slightly further out. However, the general shape with a maximum at a distance < z r and an asymptotic decay after the maximum remains the same for both 2D and 3D models. To support this model, it is compared to full 2D PIC simulations [blue stars in figure 3(c)]. While they qualitatively agree with the model, there are quantitative differences. This is expected [27], as the model only covers the Gouy phase and intensity falloff, but no other effects like source size shrinking, different divergence of the individual harmonics or coherent harmonic focussing [28]. When taking source size shrinking into account, the model predicts slightly lower efficiencies without changing the shape of the curve. This is backed up by the 2D simulations, where no assumptions about the source size of the harmonics have been made. A detailed discussion and additional figures showing the effect of source size shrinking can be found in the supplemental material. The optimal distance of 0.2z r predicted by the model can be challenging to achieve experimentally, as it requires two targets at very close distance to each other. There are two methods, how this length can be extended. By increasing the focal spot radius w 0 , the Rayleigh length z r will increase with ∝ w 2 0 . While the fraction of the Rayleigh length for the optimal position stays the same, it would allow for a larger actual distance between the targets. This comes with other challenges however, as the required intensity for SHHG may not be achievable with a large focus spot. In this case, refocusing of the harmonic beam between the targets is an alternative. Using an ellipsoidal or two parabolic mirrors, the beam after target A could be reflected and refocused on target B. After reflection from the mirror the Gouy phase evolves in the opposite direction, reducing the relative phase shift between the colours. This would allow for arbitrary large physical distances between the targets, only limited by experimental constraints. Positioning target B behind the focus of the refocused beam would be equivalent to negative distances between the targets, as also shown in figure 3(c). We note that the model predicts in that case an even faster decreasing efficiency without a maximum, which suggests that this case should be avoided in experiments.

Parameter scan
Until now, only one set of parameters (a 0 = 10, n = 200n c , Θ = 45 • , no scale length) have been used to show the principle. To investigate the influence of the various laser-plasma parameters on the optimal distance, more than 200 additional 1D PIC simulations have been conducted. For every set of parameters simulated there was a clear optimal distance which maximized the efficiency. Depending on the specific interaction parameters, the optimal distance varies slightly. As long as the similarity parameter S = a 0 n c /n [3] remains the same, the position of the optimal distance will not change. If S is changed, e.g. for different a 0 at constant n, the optimal distance will vary between 0.1z r and 0.36z r for a range of a 0 = 8, . . . , 20 (corresponding to a S = 0.04, . . . , 0.1), as shown in figure 4. Introducing a preplasma scale length L p , the absolute efficiency even for a single target is expected to increase [4]. A series of simulations for the experimentally relevant cases of L p = 0.05λ 0 and L p = 0.1λ 0 [29] have been conducted and the results are shown in figure 5. Compared to the case without preplasma, the optimal distance is shifted closer to the Figure 3. Dependence of the harmonic generation efficiency η = I AP /I L (normalized to the single pass efficiency η S ) on the distance between the two targets. The AP is obtained by applying a highpass filter which cuts out everything below the 20th harmonic. Each data point corresponds to two 1D PIC simulation, one for the interaction with the first target, and one for the interaction with the second target. The pulse between the two targets has been modelled to reflect the intensity falloff only (a), the Gouy phase shift only (b) and the combined effect of both (c). The combined effect is also compared to full 2D simulations. first target. The magnitudes differ as the single target efficiency used for normalization also increases with the addition of preplasma. The angle of incidence did not have a significant influence on the position of the optimal distance.
In this paper, the focus of the incident laser pulse was always on the first target. Another possibility would be to focus the beam on the second target, where the first target would act as a plasma mirror also generating high order harmonics. In this setup, there is a further constraint on the distance between the two targets. The intensity on the out-of-focus first target must still be high enough to generate at least 10%, better 20%-30% of second harmonic intensity [30]. In case of a peak intensity of a 0 = 10 of the incident laser pulse, for the best second harmonic conversion efficiency the first target should not be more than 0.5z r out of focus. Full 2D simulations were performed for this case and otherwise same parameters as the other simulations in this paper. The highest efficiency gain of 8.5 was achieved for the first target at 0.3z r in front of the focus, which is comparable to the efficiency where the focus was on the first target and the second target the same distance behind focus.
In summary the Gouy phase shift changes the relative phase between multiple colours in a Gaussian beam while propagating. This also changes the resulting waveform, which has implications for multipass high order harmonics generation. As the waveform changes with the propagated distance, there exists an optimal distance between the two targets, where the harmonic generation efficiency on the second target is maximized. This optimal distance depends on the exact geometry, laser and plasma parameters, but can be found for a broad range of parameters in the region between 0.1z r and 0.5z r . To find this distance, a series of 1D PIC simulations has been conducted, where the Gouy phase and intensity are modelled for propagation in 2D and 3D. These simulations are supported by full 2D PIC simulations. For the specific interaction parameters of a 0 = 10, n = 200n c and L p = 0 used in this paper, it is found that the harmonic generation efficiency peaks at 30x for 0.2z r or 0.5z r for 3D and 2D propagation geometry, respectively. It is not optimal to just minimize the distance between the two targets, as one would assume purely from the intensity falloff alone. This highlights the importance to consider the difference in the evolution of the Gouy phase when comparing 2D PIC simulations and experiments.