Wigner function analysis of high harmonic generation in atoms

The Wigner function provides the expression of phase space dynamics in quantum mechanics. In an application to non-linear optics, we describe its construction from time-dependent wavefunctions generated by numerical simulations, performed in both one and three dimensions, of the interaction of atoms with strong laser fields. From this analysis, the elements of the semi-classical three-step model of high-harmonic generation are extracted directly from the solution of the time-dependent Schrödinger equation. It is demonstrated that information contained in the Wigner function that is not relevant to this analysis may be removed by spatial filtering and state-projection operations.


Introduction
High harmonic generation (HHG) in atoms results from the motion of a single active electron in the the strong electromagnetic field produced by a laser. A precise theoretical understanding of this motion is essential for the fine-tuning and optimisation of HHG sources. Semi-classical methods and quantum mechanical models based on the strong-field approximation have been used extensively to gain insight into the nature this electronic motion and to optimise the interaction parameters. In order to obtain quantitative agreement with experiment in systems for which the single active electron approximation is valid, numerical solutions of the timedependent Schrödinger equation (TDSE) have been the subject of significant development [1][2][3][4][5][6][7][8][9]. Such methods can, in principle, reproduce any experimental observable, but their technical complexity obfuscates any intuitive picture of the underlying dynamics. In contrast, the semi-classical three-step model [10] describes the electronic motion by using Newtonian trajectories, beginning and ending at the nucleus at well-defined times. The use of trajectories is central to the qualitative interpretation of HHG and plays a significant role in phase matching schemes as well as various enhancement mechanisms [11][12][13]. Numerous theoretical studies have been performed examining the time profile of HHG emission and its dependence on laser parameters [14][15][16][17][18][19].
The quantum mechanical description of phase space trajectories is most naturally expressed within the Wigner function formalism [20]. The Wigner function is often useful in identifying how the behaviour of a quantum system deviates from that of a classical system. Recent applications of the Wigner function to quantum physics include systems as diverse as the harmonic oscillator, superconductivity, plasmas and electron transport in materials [21].
Wigner functions and Wigner-like functions, have been applied to intense laser interactions in a variety of contexts. The groundwork for this application was laid by Czirkáj et al through the study of tunnelling in the presence of a static uniform field [22]. Subsequently, it was used for a time-frequency analysis of HHG [23] and a phase space analysis of double ionisation [24]. More recently, Wigner analysis has been utilised to investigate the time-energy dependence of strong-field ionisation [25,26] and the interference between trajectories born at different times [27]. Wigner functions have also been used to compare the effects of long-range and short-range atomic potentials on strong-field ionisation [28] and has been demonstrated to be a useful tool for determining bound state populations [29].
The advantage of a Wigner analysis of HHG is that it facilitates a comparison between a 'classical system', which is the trajectory model, and the quantum mechanical solution that comes from the TDSE. The Wigner function allows phase space features of semi-classical analysis, such as the mapping between the momentum of a returning electron and its time of return, to be brought naturally into a quantum mechanical framework. In this way, predictions of the semi-classical model can be tested directly.
In this article, we describe methods for using Wigner functions to analyse time-dependent wavefunctions obtained by integration of the TDSE. The TDSE integration is performed both for a one-dimensional model system and a three-dimensional system that utilises a spherical harmonic basis. In the latter case, a reduced two-dimensional Wigner function is computed rather than the full six-dimensional function. The reduction in dimensionality results in a function that is easier to visualise and substantially reduces computational requirements.
By applying the analysis to a three-dimensional wavefunction we avoid approximations that are inherent in a one-dimensional potential. It has been shown previously that one-dimensional simulations are insufficient for producing quantitatively accurate results [30]. Although not considered here, a further benefit of the threedimensional simulations is the ability to extend the analysis to fields with elliptical polarisation, or ultra-strong fields where magnetic effects result in electron motion along the direction of laser propagation.
The Wigner function provides a phase space depiction of the electronic evolution by identifying correlations between position and momentum. From this description, direct comparisons can be made between the simulated data and both the qualitative and quantitative predictions of semi-classical approaches.
This analysis reveals that the Wigner function can be complicated by unhelpful interference between the bound and continuum parts of the wave function. To reduce this effect the Wigner function is modified by applying a spatial filtering window that reduces these oscillations. It is also shown that the removal of the bound state contributions to the wavefunction prior to computation of the Wigner function can be used to highlight the semi-classical trajectories of the continuum electrons.

Wigner function
In quantum mechanics, the phase space description of a state is expressed through its Wigner function. For a three-dimensional wavefunction, ψ(r), this is a six-dimensional function of the position and momentum coordinates, defined by Originally formulated by Eugene Wigner [20] to describe deviations from classical behaviour in quantum systems, this function can be thought of as describing a particle's presence in phase space. Strictly, the Wigner function cannot be interpreted as a probability density in the same way as the wavefunction because it can take negative values. Bounds on W can be obtained from the requirement that the wavefunction has unit norm, which for a three dimensional system requires that ( ) The usual position and momentum probability distributions are recovered by integrating the Wigner function over the relevant coordinates, so that The expectation value of any linear combination of operators r and p can be obtained by integrating the Wigner function along a line in phase space. In this work, a time-dependent wavefunction is obtained by the numerical integration of the TDSE. The phase space evolution is then determined by computing the Wigner function directly from the wavefunction at regular time intervals. These computational methods are described in 3.1 and 3.2.

Computational methods
3.1. Application to a one-dimensional wavefunction One-dimensional TDSE simulations are performed in the length gauge using the Crank-Nicolson method. For a time-dependent electric field E(t) and a potential V(x), the one-dimensional TDSE is For these simulations the wavefunction is represented on a linear grid with spacing Δ=0.05 au. The soft-core potential, = - , is used, which has the same ground-state energy, in atomic units, as a hydrogen atom.
In one dimension, the Wigner function is which is the Fourier transform of the Hermitian function For ψ represented by discretisation on a regular grid, this can be computed efficiently by a Hermitian Fast Fourier Transform (HFFT) of the array y y where 0j<n/2 for n grid points and F i,j is nonzero only when < -+ j i n i min , 1 ( ). The HFFT is taken along the axis of F corresponding to the j index. The resulting Wigner array, W ij =W(x i ,p j ), is defined on the matching momentum grid, defined by while the position grid points, x i , are the same as those used to generate the wavefunction.

Application to a three-dimensional wavefunction
Three-dimensional TDSE simulations are performed using the method described in [1]. This method utilises a spherical harmonic expansion for the wavefunction where the magnetic quantum number is restricted to m=0 because cylindrical symmetry of the wavefunction is guaranteed for linearly-polarised laser fields.
Computation of the full six-dimensional function Wigner function, W(r,p) is both expensive and unnecessary. The semi-classical dynamics may be obtained from a reduced Wigner function in which only the information along the axis of the laser field polarisation is computed. Taking this axis to be the z-axis, this reduced Wigner function is obtained by integrating over all transverse coordinates, yielding The reduced function, W z (z,p z ), captures correlations between the position and momentum along the polarisation axis, as in the one-dimensional simulations, but from quantum mechanical amplitudes derived from the full three-dimensional simulation.
Integration over the momentum variables in (11) can be performed immediately. Taking an intermediate function, W(r,p z ), that is reduced only in the momentum coordinates, we obtain where for convenience, the s denotes the z-component of the vector. The azimuthal symmetry is exploited by casting the remaining integrals over x and y in cylindrical coordinates, so that ò ò Applying the spherical decomposition of the wavefunction in (10), and substituting it into (14) and (16), the reduced Wigner function is written as where 'c.c.' denotes the complex conjugate, arccos , 20 arccos . 21 The integral in (19) is computed directly from the wavefunction. By evaluating the wavefunction separately at (r 1 ,θ 1 ) and (r 2 ,θ 2 ) there is no need to further expand the summations in (19), avoiding any expressions involving products of spherical harmonic functions. However, the integral over ρ cannot be performed while ensuring both r 1 and r 2 correspond to grid points on which the original numerical wavefunction is defined, necessitating the use of interpolation. It is important in this case to ensure that the two wavefunction evaluations are treated symmetrically in order to avoid artefacts in the computed result.
The integral in (18) is evaluated as the hermitian Fourier transform of the function f (z,s), which is the same numerical procedure encountered in the one-dimensional case. The evaluation of the integrals in (19) and (18) are straightforward, but there is a nontrivial computational workload associated with their evaluation due to their nesting. This is a manageable computational task but can be time consuming for large grids. Time can be saved by carefully choosing the domain of s values for which f (z,s) is calculated. This can be achieved, for example, by choosing a vanishing window function (see for example, (23) and (22)), that truncates the integration over s.

Spatial windowing
Since s may be arbitrarily large in (6), the Wigner function couples all parts of the wavefunction. This allows for complete symmetry between the momentum and position aspects of the Wigner function, but in practice it produces undesirable effects because there are parts of the wavefunction that we would like to consider separately. For example, a typical laser-atom interaction allows alarge fraction of the initial electron wavefunction to remain bound to the nucleus and there is no physical interest in its interactions with parts of the time-dependent wavefunction moving far from the nucleus. This is also true of electrons on spatially separated trajectories.
To see why this causes a problem, consider a wavepacket at x=a and another at x=b. This results in a Wigner function that is non-zero and rapidly oscillating in a region around x=(a+b)/2, even though the electron has no 'presence' there. If there is a very large proportion of the wavefunction still in the ground state, as is often the case, this pollutes the Wigner function at all parts of the x,p-plane through which the electron has spread. This is illustrated in panel (a) of figure 1 where the Wigner function is shown for an electron during a few cycle laser pulse. A rapidly oscillating background is clearly visible, extending over 100au from the origin. This is of no physical relevance and obscures the important features, which are the solid red regions appearing as thin lines.
A solution is to use a windowing function, w(s), in (6) The windowing function that is used is a Gaussian function, These show a marked reduction in the background and make the ionised features of the electron significantly easier to identify. The trade-off of this procedure, especially in the σ=20 case, is a slight blurring of the red lines that appear sharp in panel (a). The application of the window function breaks the symmetry in the treatment of the x and p x coordinates. This is evident in the rapid oscillation that persists in the regions between trajectories that are separated in p rather than in x. These oscillations remain present in all panels regardless of the window size that is used.

Comparison to classical trajectories
To compare the trajectories of the semi-classical model to the quantum mechanical behaviour, we apply the Wigner function to one-dimensional TDSE simulations that use an idealised electric field with a short threequarter cycle sine squared ramp-up followed by 1.25 cycles of oscillation at constant amplitude. The form of the field is sin sin for sin for 26 where Ω=ω/2 determines the ramp-up rate and t r =π/2Ω is the ramp duration. The amplitude and frequency are E 0 =0.12 au (5×10 14 W/cm 2 ) and ω=0.03 (1500nm) respectively. The ramp-up is chosen to ensure the wavefunction evolves smoothly, while the 1.25 cycles at constant amplitude are sufficient to ensure that all returning classical trajectories from the first half-cycle can be completed. These laser parameters correspond to the long wavelength and high intensity regime, so the TDSE solution should not differ greatly from the strong field approximation (SFA). The classical and quantum models are compared at three key points relating to the trajectory that gives maximum return momentum:

Ionisation
From the semi-classical model, the electron returns to the nucleus with the maximum possible energy if it is ionised shortly after the electric field reaches a peak value (see inset in figure 2). The main plot in figure 2 shows the Wigner function at this time. The streak in the upper-right quadrant of the plot indicates that parts of the electron wavefunction are ionised as the electric field increases, drawing them away from the nucleus, in the opposite direction to E due to the negative charge. The changing thickness of the streak is determined by the variation in the rate of ionisation. The observed variation in thickness is consistent with the prior evolution of the electric field magnitude, with faster transition to the continuum occurring at times of greater field strength.
Qualitatively, this is very similar to the Wigner function of electron tunnelling in a uniform static field [22], supporting the argument that ionisation can be considered a quasi-static tunnelling process is this regime. Figure 2 demonstrates a well-defined relationship between position and momentum and lends legitimacy to the central feature of the semi-classical model, that ionisation at a particular point in time leads to a deterministic trajectory. This is not surprising from a quantum mechanical point of view, since it can be understood from the principle of stationary phase. The Wigner function, however, allows this feature to be extracted directly from the time-dependent wavefunction.

Continuum evolution
The classical trajectory reaches is maximum distance from the nucleus when the electric field attains the magnitude and opposite direction to what it was at ionisation, for which the net change in electron momentum, ( ) , equals zero. The distance that the classical electron has travelled at this point is shown by the solid grey line in figure 3, which is approximately 150au. The Wigner function intersects with the p=0 axis at almost precisely the same point as the classical model. This is noteworthy given that, unlike the semi-classical model, the TDSE simulation takes into account the atomic potential during the propagation. In the regime considered here, the trajectories are relatively long and take the electron far from the nucleus, a distance of nearly 150au. This corresponds to a regime in which the role of the binding potential in the dynamics is expected to be negligible.
The portion of the Wigner function that remains in the top right quadrant of figure 3 indicates that a significant fraction of the continuum electron retains a positive momentum at this point. These portions of the Wigner function correspond to the long trajectories for which the electron travels further from the nucleus, as well as trajectories that do not lead to re-collision.

Return to the nucleus
The final stage in the semi-classical three-step model is the recombination of the electron with the nucleus. For the optimal classical trajectory, this occurs as the electric field is decreasing in magnitude. The momentum that the classical electron would have at this point is shown by the horizontal grey line in figure 4, with a value of approximately -5 au.
The return momentum of the quantum state can be seen from the intercept of the red line in the Wigner function with the p-axis. Again, there is excellent agreement with the semi-classical model, although the trajectory extracted from the quantum mechanical model gives a slightly higher return energy.

Time-dependent emission
The recombination step of the semi-classical model requires that the photon energy emitted by the returning electron is equal to its kinetic energy plus the ionisation energy of the parent atom. A relationship between return time and return momentum implies, therefore, a time-dependent emission spectrum. We can calculate this emission spectrum from the TDSE simulation and compare it to the spectrum predicted by the intercepts of both the Wigner function and the classical trajectory with the line x=0. To calculate the time-dependent emission we apply a window function to the dipole acceleration prior to obtaining the spectrum. The most satisfactory window function for this purpose was determined through trial and error to be the Gaussian function with τ=3 au. To obtain the frequency at time t, the window function is applied to the dipole acceleration by the substitution This corresponds to a Gabor transform of the dipole acceleration [19]. The resulting spectrum is a function of time, t, and frequency, ω.
The time-dependent spectrum for the TDSE simulation is shown in figure 5. Several distinct regions of emission can be clearly identified, with the boldness of the colour indicating the emission intensity on a logarithmic scale. The arch between 250au and 350au is consistent with the expected emission for a semiclassical electron. This electron would be born into the continuum between the first time the electric field reaches its peak value, which corresponds to a return time at the trailing end of the arch, and the subsequent point at when it next reaches zero value.
The peak of the arch corresponds very closely to the trajectory we have analysed. Figure 5 indicates peak emission occurs just after t=300 au, which is slightly later than the semi-classical prediction at t=297 au. For a return momentum of p and ionisation energy I 0 , the predicted emission energy is For I 0 =0.5 au and a p intercept of 5au, the predicted energy is 13au, which is the same as the peak energy shown in the time-dependent frequency plot. In this regime it is unsurprising that the quantum mechanical simulation validates, to a large extent, the assumptions of the semi-classical model. By analysing trajectories within a quantum mechanical framework, the Wigner function can provide insights into the domain of validity of the semi-classical approximation when the parameters of the laser field are changed.

The Wigner function in the near-threshold regime
In this section, we apply the z-axis formulation of the Wigner function to a 3D TDSE simulation of the generation of harmonics in the near-field regime. In contrast to the simulations described in Section 4, it is not clear that the semi-classical model of trajectories will be as evident in the quantum mechanical description in this regime. Certainly, once the energy of the emitted radiation is less than the ionisation energy, it becomes impossible to view the electron evolution as a classical trajectory through the continuum. This exclusion arises because this model assumes the emitted energy is the sum of the ionisation potential and the kinetic energy of the returning electron. It is also expected that in this regime, the role of the Coulomb potential has a greater influence over the continuum evolution of the active electron then elsewhere.

Simulation parameters
For simplicity, the simulated atom is hydrogen, with the potential V(r)=−1/r and the electron initially occupying the 1s ground state. The pulse duration is 10 cycles, with a peak intensity 6×10 14 W/cm 2 and wavelength 800nm. The electric field is assumed to be of the form This low intensity has the effect that the simulations require only 21 terms in the angular momentum expansion of the time-dependent wavefunction. A fixed grid of radial dimension 100au is used that includes a masking function to minimise reflections from the boundary. For these laser parameters, the maximum distance a classical electron can achieve and still return to the nucleus is less than 25au. The maximum return energy, based on the standard cut-off law [31], is around 0.41au, which is less than the ionisation potential of the atom.

Removal of bound states
In cases involving low intensity interaction such as this, the electron remains almost entirely bound to the atom throughout the interaction. When this happens it is difficult to discern useful information from the Wigner function because the ground state dominates all nearby portions of the wavefunction, resulting in rapid oscillations that mask any underlying structure. In section 3.3, the use of a window function was described that can be used to solve a similar problem, where spatially separated parts of the wavefunction produce unhelpful oscillations in regions where the probability density is negligibly small. In this case, however, a window function is not suitable because the relevant parts of the wavefunction are not necessarily separated in space. A better solution is to remove the ground state and any other problematic bound states, from the wavefunction prior to the computation of the Wigner function. | are the bound states that are to be removed. To understand the effect of this procedure, the Wigner function reduced to its z-components, W(z,p z ), is plotted for the 1s, 2s, and 2p states of a hydrogen atom in figure 6. The ground state Wigner function is non-zero at large p z along the line z=0, which is due to the divergence of the Coulomb potential at the origin. The 2p state has a negative (blue) region near the centre. Often, negative values of the Wigner function are cited as evidence of quantum behaviour [32]. Here, it simply reflects the fact that the Wigner function, by itself, cannot be interpreted as a probability distribution. The procedure can result in a similar negative-value region in the resulting reduced Wigner function because removal of the ground state tends to leave 2s and 2p as the dominant states. It is important to note that the Wigner function is not linear with respect to the wavefunction, so that removal of any state does not have the same effect as the subtraction of that state's Wigner function from the complete Wigner function of the electron.

Analysis of Wigner functions
For the present simulations, more than 99.9% of the probability distribution of the electron remains in the ground state throughout the evolution. Consequently, only the ground state is removed from the wavefunction prior to computing reduced Wigner functions.
The reduced Wigner function at the halfway temporal point of the simulation is shown in figure 7. For comparison, the Wigner function is shown as computed with the full wavefunction (panel a), as well as with the ground state removed (panel b). The electric field at this point in the simulation is zero, with a positive derivative. In panel (a), there is strong variation in the magnitude of the Wigner function, with the most significant variation near the origin. In order to make visible the regions distant from the origin, the colour mapping is on a symmetric logarithmic scale, with the region close to zero treated linearly. In panel (b), there is no such variation in magnitude and a linear colour mapping suffices. The z-axis scale is sufficiently large such that the boundary exceeds the distance an electron can travel and return to the nucleus.
The difference between the plots with and without the ground state is dramatic. If the ground state is retained, the trajectories are almost entirely obscured by the oscillations and the position and momentum distributions of electron wave-packets are difficult to extract. In contrast, panel (b), where the ground state has been removed, shows the trajectories more clearly. The recently-ionised part of the wavefunction appears as the dark red line trailing away from the origin in the upper-right quadrant of panel (b). A streak stretching across each of the upper quadrants and intersecting with the line z=0 at approximately p z =1 au can also be seen clearly. Panel (b) also shows the negative region around the origin, which is indicative of the population of the 2p state.
In comparison with the results from the high intensity simulations, prominent ripples in the Wigner function are evident, which can be attributed to bound state and low-energy continuum parts of the electronic wavefunction. These ripples are strongly time-dependent correspond to the configurations where the belowand near-threshold harmonics are generated. The absence of clear trajectories for these harmonics is entirely reasonable, given there is no semi-classical trajectory that can lead to this emission. A clearer depiction of transitions between bound states is beyond the scope of these simulations as the radial grid spacing is insufficient to provide the requisite energy resolution for states near the ionisation threshold.

Conclusion
We have developed a means by which the trajectories of the semi-classical three-step model of HHG can be analysed by applying the Wigner function to the quantum mechanical wavefunctions obtained from numerical solutions of the TDSE. This method provides a powerful tool for visualising the evolution of the wavefunction, and allows comparison with the trajectory picture while maintaining the generality of the solution of the quantum mechanical equations of motion. By applying the Wigner function to three-dimensional calculations, the pitfalls of one-dimensional simulations are avoided. The analysis can be further generalised to threedimensional electronic motion by reducing, for example, the Wigner function along separate axes.