Attoscience in Phase Space

We provide a brief review of how phase space techniques are explored within strong-field and attosecond science. This includes a broad overview of the existing landscape, with focus on strong-field ionisation and rescattering, high-order harmonic generation, stabilisation and free-electron lasers. Furthermore, using our work on the subject, which deals with ionisation dynamics in atoms and diatomic molecules as well as high-order harmonic generation in inhomogeneous fields, we exemplify how such tools can be employed. One may for instance determine qualitatively different phase space dynamics, explore how bifurcations influence ionisation and high-harmonic generation, establish for which regimes classical and quantum correspondence works or fails, and what role different time scales play. Finally, we conclude the review highlighting the importance of the tools available in quantum optics, quantum information and physical chemistry to strong-field laser-matter interaction.

physics and nonlinear dynamics, and widely employed phase space variables are, for instance, positions and momenta, or angles and angular momenta. Its mathematical origin, dating from 1838, can be attributed to Liouville [2], and its first application to mechanics was made by Jacobi in 1842 [3]. However, the concept of describing the dynamics of a system as a single trajectory moving through multidimensional space was developed many decades later by Poincaré [4] (for a historic review on the subject see [5]).
The quantum phase space was introduced much later, by Wigner, together with the quasiprobability distribution named after him [6]. Since then, quantum phase space distribution functions, constructed using non-commuting operators, have become widespread. A key advantage is that they allow one to employ complex-number functions instead of dealing with operators. Furthermore, they provide valuable insight in quantum-classical correspondence [7], within the constraints posed by the uncertainty principle and its generalisations. However, there are different phase space distribution functions, whose applicability may suit specific problems better than others. This ambiguity stems from the fact that there are different rules for associating non-commuting operators to scalar variables [8,9]; for pioneering work exploring operator ordering in connection with quasiprobability distributions see also [10,11]. For instance, due to their smooth behavior, Husimi distribution functions are popular in the context of nonlinear systems and classical chaos [12], while Wigner quasiprobability distributions, due to the information they provide on non-classical effects and quantum corrections, are widely used in quantum optics [13,14]. Other applications of the Wigner function include optical propagation in waveguides [15], and the computation of angular momentum states [16], which can also be used to model two-level atoms [17][18][19]. The Glauber-Shudarshan P functions [20,21] are also hugely popular in quantum optics, as they are very convenient for normal-ordered products of creation and annihilation operators.
Quantum phase space distribution functions play a major role in quantum optics [13,14] and quantum information [22,23]. This popularity has been triggered by the description of the electromagnetic field modes as quantum harmonic oscillators, for which distribution functions have been especially tailored (see, e.g., the discussion in [8,9,20,21]), and the central interest in the definition of non-classical states of light [24]. Furthermore, due to being formulated in terms of density matrices, quasiprobability densities are well suited for investigating decoherence and the influence of the environment [19,25,26]. Quasiprobability distributions have also been explored in connection with logical gates [27,28] and their classical simulation [29], and coherent-state superposition [30,31].
Other traditional areas in which the quantum phase space is widely used are those dealing with large systems [32], such as chemical physics [33,34] and cold gases [35]. In this case, the huge amount of degrees of freedom makes a full quantum-mechanical treatment prohibitive. Therefore, crucial questions are what degrees of freedom need full treatment and which ones can be approximated, what kind of fluctuations and deviations from the classical picture are expected, and whether there are semiclassical limits one can take into consideration without losing essential information about the system's dynamics.
In the study of complex molecular systems, for instance, it is common to apply mixed classicalquantum methods, which describe less relevant degrees of freedom classically, more relevant degrees of freedom quantum mechanically, and couple them via effective potentials [36,37]. One may also consider systems coupled to baths, whose dynamics are simplified [38]. Alternatively, one may develop semiclassical methods, in which swarms of classical trajectories are employed to construct quantum propagators (see, e.g., [39][40][41][42][43] and the reviews [33,44]). Further approximations may then be applied, such as the smoothening of highly oscillatory terms [7], and linearised semiclassical approximations, in which the main contributions stem from trajectories whose phase space coordinates are close enough [40]. Truncated Wigner approximations are also widely used to describe chemical reactions (for an early example see [40]).
Furthermore, there are many perturbative approximations which incorporate quantum fluctuations around classical limits, such as the semiclassical or truncated Wigner approximation (TWA). For an early discussion of phase space methods in which stochastic evolution equations, including the TWA, are applied to Bose-Einstein Condensates see [45]. The key idea is to embed quantum fluctuations in the initial quasiprobability distributions, which are then evolved classically. This means that quantum corrections appear only through initial conditions. Classical evolution is desirable for large systems as, typically, classical equations of motion scale linearly with its degrees of freedom, while quantum methods scale exponentially. The TWA has been explored in the context of Bose condensed gases perturbed from thermal equilibrium [46] in a wide range of scenarios such as the evolution of macroscopic entangled states [47], quantum fluctuations in condensate oscillations [48], non-classical effects arising from the splitting of condensates [49,50] and vortices in reflecting Bose-Einstein condensates [51]. It works for weakly perturbed systems.
In contrast, in strong-field laser-matter interaction and attosecond science, the phase space picture and its tools have not become mainstream. This is actually surprising for the following reasons: -Classical and semi-classical trajectory-based methods have been used as interpretational tools for quantum effects since over two decades, and helped establish the key paradigms in the field. Therein, laser-induced rescattering and recombination of an electron with its parent ion play a vital role by providing an intuitive interpretation of many strong-field phenomena (see, e.g., [52][53][54] and the special issue [55]). If recombination with a bound state occurs, the kinetic energy acquired by the electron in the continuum is released in form of high-order harmonic radiation, while rescattering will lead to high-energy photoelectrons. If the returning electron rescatters elastically, high-order above-threshold ionisation (ATI) will occur [56,57]. Alternatively, if, upon recollision, it passes on part of its kinetic energy to the core, it may release other electrons, leading to nonsequential double and multiple ionisation [58,59]. This orbit-based picture has been hugely successful in describing the physics of the problem, with far reaching consequences. For instance, the shapes of the high-harmonic and high-order ATI spectra, with a long plateau followed by sharp cutoffs whose energy positions are a multiple of the ponderomotive energy U p [53,57], proportional to the driving-field intensity, can be explained using the laser-induced rescattering picture. Furthermore, because ionisation and recombination occur at very specific times within a field cycle, they can be used for steering subfemtosecond electron dynamics, for generating attosecond pulses or bursts of electrons (see, e.g. [60,61], and the special issue [62]), and for subfemtosecond imaging of matter [63][64][65]. For a given energy, there is usually more than one quantum mechanical pathway along which the electron may interact with the core, so the corresponding transition amplitudes interfere. Quantum interference has many attosecond imaging applications, using high-order harmonic generation or photoelectrons (for reviews see, e.g., [66] or [67]). -The interaction between the continuum and bound states is of vital importance, as strong-field ionisation and laser-induced rescattering or recombination play a key role in explaining strongfield phenomena [53, 54, 56-58, 67, 68]. Yet, tools that could provide rigorous and/or accurate information about non-classicality or boundaries of bound-continuum dynamics are underused.
Overall, the use of phase space has been twofold: either the classical phase space was employed to delimit bound-continuum boundaries, highlight regular or chaotic behavior, and analyse different rescattering regimes, or quantum phase space distributions have been employed to assess classical or non-classical behavior and provide initial conditions for other methods. Often the intuitive picture obtained by classical methods is compared with the outcome of ab-initio computations or other approaches.
In the present article, we provide a review of attoscience in phase space. We will start with a brief overview of its use in strong-field and attosecond physics (Sec. 2), with emphasis on the different phenomena and how specific methods, quantum, classical and semiclassical, have been used over the years. Whilst it is not possible to draw a linear timeline, we have grouped such studies according to common physical or methodological aspects, in order to set an overall landscape. Subsequently, in Sec. 3, we delve deeper into our own work, and use it to exemplify how classical and quantum aspects of the phase space may be used in attoscience. We start by providing a brief statement on the methods (Sec. 3.1) employed and, in the ensuing sections, focus on common strategies rather than phenomenon or publication. This includes using phase space methods to identify different phase space configurations and bifurcations (Sec. 3.2), distinguishing between classical, semiclassical and quantum regimes (Sec. 3.3), and performing an in depth analysis of the different time scales that arise in our investigations (Sec. 3.4). Finally, in Sec. 4, we conclude the article by placing phase space studies into a broader context. In particular, we highlight the importance of bringing the toolkit available in quantum optics to attosecond physics in view of recent trends and developments.

Overview
Historically, the phase space has been applied to a variety of phenomena in attoscience, along the following research lines: Free-electron lasers and stabilisation, tunneling, rescattering in oneelectron systems and correlated multielectron dynamics. These research lines often overlap, and a key common aspect is to try to understand and control attosecond electron dynamics in greater depth. They are briefly discussed below.

Free-electron lasers and stabilisation
Phase space tools have been first used in attosecond science and related fields in the 1980s, in the context of free-electron lasers (FELs) [113,114]. These seminal papers aimed at bridging a gap that existed between fully quantum electrodynamic descriptions of electrons in a FEL and widespread classical descriptions of these dynamics. This was an important milestone as the classical description of electrons in a FEL are expected to become inaccurate in the XUV/X ray regime, for which quantum fluctuations are important [113,114]. For that purpose, Wigner quasiprobability distributions were used and it was shown that, in the classical limit, the former descriptions were recovered. Following that, there have been phase space studies to determine the boundary between classical and quantum behaviour [115][116][117][118][119]. These studies have been motivated by the prospects of developing a Quantum FEL, which should exhibit a narrower linewidth and better temporal coherence than its classical counterpart [120]. In particular, one- [117] and three-dimensional [118] quantum models based on Wigner functions are presented. In [119], the phase space was used to establish a quantum regime, in which the system may be approximated by a two-level atom by averaging over fast oscillations. Recently, quantum effects in the FEL electrons and its gain were studied using Wigner quasiprobability distributions and the quantum Liouville equation [115].
Further work in the high-frequency regime, in the 1990s, addressed the question of non-classical behavior in atomic stabilisation. Atomic stabilisation stems from the breakdown of Fermi's golden rule for computing ionisation probabilities in strong laser fields. Roughly speaking, stabilisation is the suppression of ionisation with increasing field strength. In the 1990s, it has generated a great deal of controversy, from its definition, to the physical mechanisms behind it and its existence altogether (for reviews see [121][122][123][124]). In this context, the Kramers-Henneberger frame, in which the field time dependence is embedded in the binding potential, is widely used. For high driving-laser frequencies, one may define a double-well effective potential, known as the Kramers-Henneberger potential. Wigner quasiprobability distributions were employed to assess under what conditions stabilisation was classical or quantum [69]. They were compared to classical-trajectory computations and exhibited regions that were attributed to coherent superpositions of a few bound states of the effective Kramers-Henneberger potential, thus highlighting the role of quantum interference [70]. Further work explored how quantum effects in stabilisation depend on the pulse shape and on the effective Kramers-Henneberger potential by using Wigner and Husimi distributions [71]. Much later work relates stabilisation to trapped trajectories and elliptic islands in a chaotic region via a classical phase space analysis, and highlights a hidden short-time nature of stabilisation [72].

Tunneling
In the low-frequency regime, phase space quantum distribution functions have been employed to assess tunneling ionisation dynamics since the early 2000s. Tunneling is crucial for strong-field and attosecond physics, and the question of with what velocity, and at what point in space an electron reaches the continuum, as well as whether one may define finite tunneling times, has attracted a great deal of attention (see, e.g., [125][126][127][128][129][130][131][132][133][134] for a wide range of approaches and points of view). Since the phase space allows for an intuitive view of the classical-quantum correspondence, it is ideally suited to such questions.
Early studies of tunneling ionisation in phase space observed a tail for the Wigner quasiprobability distribution of a system in a static or quasi-static, low-frequency field. This tail has been associated with tunnel ionisation as it crosses from a bound phase space region to the continuum through classically forbidden regions [135]. Its slope has been employed to define a tunnel trajectory, which was first computed by [73] using an analytical model of a zero-range potential in a static field. Further work, a decade later, [79] investigated how the slope of the Wigner function behaved with regard to the potential being short or long range. Both publications focused on the agreement between the tail of phase space quantum distributions and a classical-trajectory picture, which were shown to match far away from the core. Nonetheless, quantum interference fringes associated with tunneling events at different times were observed in both publications. Close to the core, the tail of the Wigner function follows the separatrix and crosses into the continuum either via over the barrier or tunnel ionisation [76]. Recently, Wigner quasiprobability distributions have been employed to reconstruct the tunnel exit, which is an important parameter in determining the tunneling time [136]. Further work by the same group addresses the influence of quantum interference and over-the-barrier ionisation on classical-quantum correspondence when an electron is freed into the continuum [137].
For systems with more than one centre, such as in diatomic molecules, Wigner quasiprobability distributions have been employed in the context of enhanced ionisation [74,75,78]. Roughly speaking, enhanced ionisation means that, for specific internuclear separations, ionisation rates in a stretched molecule are considerably higher than those in an atom with a similar ionisation potential [138]. In [74,75] it has been shown that there are intra-molecular momentum gates in phase space, which facilitate population transfer within the molecule and to the continuum. They have been attributed to the non-adiabatic response of the molecule to a low frequency field. Further work, however, showed that the momentum gates are intrinsic to the molecular system, and exist even for static fields, or no fields at all. Thereby, quantum interference plays an important role by providing a bridge for quasiprobability to flow from one molecular well to the other [78], and the frequencies can be estimated for double-well potential models [139]. Further studies of non-adiabatic effects and bifurcation in strong-field ionisation were conducted in [77].
Tunneling dynamics in strong fields has also been looked at in the context of initial-value representations (IVRs) [76,140]. In initial-value representations, the boundary problems that arise in semiclassical theory are replaced by averages over initial phase space coordinates, which are used to construct wave packets. These wave packets are then evolved in time, guided by ensembles of classical trajectories. IVRs are employed in many areas of science, for instance quantum chemistry, chaos and nonlinear dynamical systems, and are very powerful approaches due to their scalability and absence of cusps and singularities. For key references see, e.g., [39][40][41][42][43] and the reviews [33,44]. However, there has been considerable debate whether these approaches can be used to model tunnel ionisation, as they employ ensembles of real classical trajectories to construct wave packets [141][142][143][144].Tunneling may manifest itself in position space, as transmission, or in momentum space, as above-the-barrier reflection, and it not being accounted correctly will cause semiclassical IVRs to degrade for longer times [144]. Nonetheless, because the top of a potential barrier can be approximated by an inverted harmonic oscillator, tunneling has been found to work well near this threshold [135]. To deal with this, one may either focus on rescattering only and place the initial electronic wave packet away from the core [107,108], or employ short times and IVRs with effective potentials that account for quantum corrections, such as the Coupled Coherent State (CCS) method [76]. The phase space has also been employed to develop path integral approaches [145] that incorporate the residual potentials and the driving field on equal footing, such as the Coulomb Quantum Orbit Strong-Field Approximation (CQSFA) [146][147][148][149] and the semiclassical two-step (SCTS) model [150,151], with emphasis on quantum interference and photoelectron holography. For a review see [67].

Rescattering in one-electron systems
Because most strong-field phenomena can be explained as laser-induced rescattering, one must understand how it manifests itself in phase space. Although structures associated with rescattering have already been identified in [73], closer scrutiny happened only in the 2010s. In [152], distinct interference patterns in Wigner quasiprobability distributions have been associated with different types of intra-cycle electron scattering and above-threshold ionisation. This has been extended in [87] in order to assess lower impact velocities, and to compute the bound-state population using phase space criteria. Therein, phase space signatures of channel closings have also been identified. Further work has investigated the connection between rescattering and entanglement [88].
Rescattering in phase space has also been studied in relation to other phenomena. For instance, in [79], a phase space analysis of rescattering in conjunction with high-order harmonic generation (HHG) was performed employing Wigner and Husimi distribution functions. It was shown that the rescattering wave packet exhibits a chirp, which can be extracted from the Wigner quasiprobability distribution at the position of rescattering. The HHG temporal profile given by the Wigner function strongly resembles that obtained by other means such as windowed Fourier transforms. Recent work has focused on a systematic analysis of the orbit-based rescattering picture for tunneling, rescattering and HHG using Wigner functions with spatial windows in reduced-dimensionality models, and effective Wigner functions for multidimensional systems in order to facilitate the interpretation of more intricate dynamics [89].
Different types of orbits and their role in HHG [107] and ATI [108] have also been investigated using initial-value representations. It was found that irregular orbits play an important role in forming the HHG plateau. Furthermore, phase space tools have been employed to identify regions of dominant, integrable Hamiltonians, which led to HHG spectra with excellent agreement with ab-initio methods [109,110]. A key challenge in modeling HHG is that it is a coherent process that relies on the periodicity of the field. This implies that any dephasing associated with the degradation of the time evolution determined by the IVR will affect the harmonic spectra. This will play a key role if the wave packet is initially bound as tunnel ionisation will be important in this case (for discussions see [76,153]).
Subsequently, a purely classical perspective into how the presence of the Coulomb potential affects recollisions in strong fields, as related to the high-harmonic spectra, is provided in a series of publications [80,81,83]. Therein, classical phase space arguments have been used to show that Coulomb focusing enhances delayed recollisions and increase their energy. These recollisions occur along periodic orbits whose energy are higher than the well-known value of 3.17U p [81]. Nonetheless, in [80], a fully classical method that considers the Coulomb potential in the continuum is employed to explain why the standard cutoff law works. A set of periodic orbits stemming from a resonance with the field are linked to laser-induced recollision, whose maximal energy approaches the standard HHG cutoff in the high-intensity limit. Good agreement with the ab-initio solution of the timedependent Schrödinger equation is observed. Further work explores the extension of the cutoff upon macroscopic propagation [83]. The phase space has also been employed by us in [82] to extract different instantaneous configurations and time scales for HHG in inhomogeneous fields.

Correlated multielectron processes
In addition to one-electron problems, since the early 2000s, the phase space has been used to explore non-trivial features in correlated multielectron processes. This extends from laser-induced nonsequential double ionisation [58], which is the archetypical example of electron-electron correlation in intense laser fields [90], to the temporal profile of autoionisation dynamics in Helium [85]. For instance, in [90] Wigner quasi-probability distributions associated with the centre-of mass coordinates of a two-particle system have been compared to the outcome of a mean-field theory in order to identify signatures of rescattering and electron-electron interaction. NSDI has also been modeled for the Helium atom using IVRs beyond reduced-dimensionality models [154,155]. In particular an alternative version of the Coupled Coherent States (CCS) method that incorporates the exchange symmetry of fermionic particles, the fermionic CCS, has been successfully applied in this context [155].
A whole line of research has been devoted to investigating NSDI in a fully classical framework. In NSDI, an outer electron reaches the continuum, is brought back by the field and transfers part of its kinetic energy to an inner electron. These recollision dynamics are quite rich and can be interpreted using tools from phase space, the theory of non-linear dynamical systems, and effective Hamiltonians for each of the electrons involved. These reduced, integrable Hamiltonians were first defined in [97,98] for NDSI. Subsequently, the role of multiple recollisions on the efficient energy transfer in NSDI has been investigated using symplectic maps and similar approaches to those used in kicked Rydberg atoms, and a road map has been provided for identifying different NSDI mechanisms in [99,100]. Further work by the same group has focused on an in-depth analysis of recollision excitation with subsequent ionisation (RESI) in terms of resonances and their proximity to periodic orbits [103,104]. In particular, a sticky region in phase space arises due to the interplay of the external field and the binding potential. This region traps trajectories before ionisation, leading to time delays for the second electron. A detailed analysis of the types of periodic orbits, resonance conditions and distinct sources of chaos is provided in [104]. Interestingly, if reduced-dimensionality models are used, oscillations in the RESI yield as functions of the laser intensity are reported and attributed to resonances. However, these oscillations are washed out if more degrees of freedom are incorporated, due to chaotic transverse dynamics and additional resonances. These oscillations are distinct from those attributed to quantum interference in RESI [156][157][158][159]. Further work is related to NSDI in bichromatic, linearly polarised fields [102], and the dynamics of recollisions in fields with circular polarisation [101,105,106]. In [101] it is shown that, in contrast to previous assumptions, recollisions may occur in circularly polarised fields, by analysing the system's dynamics in a rotating frame. The physical mechanism is similar to that leading to ionisation in Rydberg atoms in microwave fields. A decade later, this topic is revisited and several associated time scales are analysed in detail [105,106]. In particular, a recollision mechanism taking place over many field cycles is reported.
One may also employ classical models and dynamical systems tools to determine modified threshold laws for correlated multielectron ionisation that account for the presence of an external field [91]. This approach starts by identifying similarities with its field-free counterpart [160], considering the two electrons in an excited compound and identifying the relevant subspaces for which correlated double and multiple ionisation may occur. This information can then be used to construct effective reduced Hamiltonians for the subspace of interest, identify existing symmetries and possible electron escape configurations, and determine which of the latter will prevail. This has been done for nonsequential double [92,95], triple [93] and multiple [94] ionisation. One can also use this method as guidance for defining effective reduced-dimensionality quantum models, so that the actual dynamics are preserved as faithfully as possible, without introducing artificial constraints or correlations [96].
Another manifestation of electron-electron correlation is autoionisation. Thereby, the quantum interference between a direct and a quasi-resonant pathway is of extreme importance. Wigner quasiprobability distributions in the time-energy domain are used to study this interference and disentangle these pathways in a transient process [85]. They provide an advantage over other methods used in time-frequency analysis, for exposing non-classical behavior in a much more explicit way. A further issue is that, in classical-trajectory models, autoionisation manifests itself as an artifact that renders the system unstable. These problems have been overcome in [86,161], which employ quasiprobability densities and phase space arguments to develop classical-trajectory models that do not exhibit these shortcomings and are consistent with their quantum-mechanical counterparts.

Selected examples
In the following, we provide a few selected examples from our own work of how phase space tools can be used to model strong-field dynamics and extract information that may not be available by other means. These are used to classify regions with qualitatively different behaviors, and understand how processes such tunneling and rescattering unfold. For details, we refer to the original publications [76,78,82,139]. We also intend to go beyond these articles, by bringing together aspects that have not been emphasised so far.

Methods
We employ both classical and quantum-mechanical tools. We focus on simplified, reduced-dimensionality models in which a single spatial dimension is taken into consideration. They provide a transparent, yet accurate picture of the system's dynamics for linearly polarised fields. We also use atomic units. One should note, however, that the choice of lower-dimensional models is not trivial in attosecond physics, as one must ensure that no artifacts such as spurious symmetries and correlations are introduced, with regard to realistic, multidimensional models. For correlated electron dynamics, this is briefly mentioned in Sec. 2.4. For one-electron problems, a discussion of how to keep the onedimensional dynamics and observables as close as possible to a three-dimensional system is given in [162,163].

Classical phase space dynamics
Classically, the phase space dynamics are described by Hamilton's equationṡ where x and p are the position and canonically conjugate momentum, respectively, and the classical Hamiltonian H is defined by In Eq. (3), V eff = V (x) + V l is the effective potential determined by the external electric field acting on the electronic wave packet. The physical picture of a time-dependent effective potential corresponds to the length gauge and the dipole approximation, which are used in this work. Thereby, V l is the potential energy determined by the laser field, which for a field without spatial dependence, reads V l = xE(t).
Here we take the binding potential V (x) to be either soft-core with λ = 1 constant, or as the short-range Gaussian potential where λ = 1/2. For diatomic molecules we consider The field is assumed to be either static, i.e., E(t) = E 0 , or as a linearly polarised monochromatic wave of frequency ω and phase φ Specifically in the work associated with inhomogeneous fields we introduce the spatial dependence of the laser field by considering and the associated inhomogeneous effective potential where β is an inhomogeneity parameter 1 . This approximation has been widely employed in the literature [164][165][166][167][168][169][170] contingent on a small inhomogeneity parameter. This classical approach is used to determine bound and continuum phase space regions, identify fixed points and analyse ensembles of trajectories under different conditions. This is done by employing key concepts of the theory of dynamical systems, some of which are briefly stated here for the sake of self consistency. For a more detailed and rigorous discussion please see e.g., [171]. Solutions of Eqs. (1) and (2) which stay the same for all times, that is, for whichẋ =ṗ = 0, are fixed points. For conservative Hamiltonian systems of the form (3), such as a model atom in a static field (E(t) = E 0 ), one may show that fixed points are centres or saddles. In phase space, these fixed points are located at (x f , p f ) = (x s , 0), where x s is the value of the coordinate x for which the effective potential V eff orṼ eff is stationary. Centres and saddles are given by the minima and maxima of the effective potential, respectively. Centres are attractive and surrounded by closed orbits, while saddles are semi-stable and help delimit qualitatively different dynamical regions in phase space. For that reason, phase space trajectories crossing saddles are called separatrices. For a time-dependent field, this line of argument is approximate, but it can be used if its frequency is low enough for the quasi-static picture to hold. This is the case in the present work.

Time-dependent Schrödinger equation and initial wave packet
As a benchmark, we use the full solution of the time-dependent Schrödinger equation (TDSE), where the length-gauge Hamiltonian readŝ with V eff (x, t) being defined as above and the hats denoting operators. We solve the TSDE in the position space, where Ψ (x, t) is the time-dependent wave function. Depending on the problem at hand, we choose different initial wavefunctions Ψ (x, 0) and potentials V eff (x). We will approximate the initial wave function by Gaussian wave packets of width γ centred at vanishing initial momentum p 0 = 0 and initial coordinate q 0 , or coherent superpositions thereof. Specifically for diatomic molecules, we will consider γ = 0.5 a.u. and q 0 = −R/2 or q 0 = R/2, in which cases the initial wave functions are given by Ψ down (x, 0) or Ψ up (x, 0), respectively. The delocalised wave function is taken to be the symmetric coherent superposition The time-dependent wave function is employed to compute quantum distribution functions and observables discussed in the next section. It will also be used to calculate the time dependent autocorrelation function as well as the ionisation rate from an initial time t = 0 to a final time t = T f , where This definition of ionisation rate was used in the seminal paper [138] in the context of enhanced ionisation of molecules, and in our previous publications [78,139].

Quantum distribution functions
The time-dependent wave function is used as input in the Wigner quasiprobability distribution This function is always real. However, it exhibits both positive and negative values, hence the name "quasiprobability". This, among other features, makes its interpretation as a simple probability distribution difficult. This is why it is often paired with a study of classical trajectories in phase space. On the other hand, its deviations from probability densities can be used to define nonclassicality. For instance, a widespread definition used in quantum optics is to seek regions for which W (x, p, t) is negative and classify them as non-classical [172]. A more restrictive definition is based on the quantum Liouville equation [13]. where are the quantum corrections to the classical Liouville equation. The quantum Liouville equation (also known as Moyal equation) may also be written more compactly as where H(x, p, t) is the system's Hamiltonian and {{·}} give a Moyal bracket [173,174] 2 . In the classical limit (setting = 0), Eqs. (19) and (21) become the classical Liouville equation, so that Q(x, p, t) = 0 and the right-hand side of Eq. (21) will be given by a Poisson bracket. In this case, the Wigner quasiprobability distribution will evolve like a classical entity. This is a useful tool for distinguishing between regimes in which quantum interference is present, but evolves classically by, for instance, following classical separatrices, and those truly quantum regimes with no classical counterpart. A widespread approach in quantum optics and cold gases, known as the truncated Wigner approximation (TWA), is to consider the classical Liouville equation with stochastic quantum corrections. This is many times required in order to deal with large systems. The TWA was first used in the context of Bose-Einstein condensates in [45]; for reviews see [32,35], and is closely related to the linearised semiclassical IVR [33,40]. Nonetheless, in some instances it may be nontrivial to compute a classical limit for the quantum Liouville (Moyal) equation (see [7] for an early discussion). One should note that quantum distribution functions may be defined using any two variables corresponding to incompatible observables, such as time and frequency. For instance, Wigner-type time-frequency distributions were employed to study HHG [176], different regimes in ATI [177] and autoionisation [85]. This approach bears some similarity with the use of windowed Fourier transforms, which is much more widespread and has been used since the 1990s to infer time profiles of harmonics and extract electron return times from ab-initio computations (see, e.g., Refs. [178][179][180][181][182][183] for early studies or Refs. [164,165,[184][185][186] for more recent publications).

Gabor and Fourier transforms
In our work we use the Gabor transform to compute time-resolved HHG. The time-resolved spectra are given by χ G (Ω, t) = |d G (Ω, t)| 2 , where Ω is the harmonic frequency, with σ = 1/3ω, is a windowed Fourier transform with a Gaussian window function, and is the dipole acceleration [187][188][189]. The standard HHG spectrum, for which all temporal information is lost, is computed as χ(Ω) = |d(Ω)| 2 , with and is recovered by setting σ → ∞ in Eq. (22).

Initial value representations (IVRs)
We have also employed initial-value representations (IVRs), where the time-dependent wave function is constructed using a basis of Gaussian wave packets in phase space guided by a trajectorybased grid. IVRs exhibit many advantages, such as providing an intuitive picture in terms of electron orbits, accounting for external fields, binding potentials and quantum interference. Furthermore, they are applicable to large systems due to the numerical effort involved not scaling exponentially with the number of degrees of freedom. Specifically, we employ the Herman Kluk (HK) propagator [42] and the Coupled Coherent State (CCS) representation [43]. Coherent states are widely used in chemical physics and quantum optics for behaving in a quasi-classical way. This allows descriptions in terms of functions of complex numbers, and using the phase space approaches mentioned in this work. The most common coherent states are those of the quantum harmonic oscillator. They have been first introduced by Schrödinger in [190], and systematically studied by Glauber [20] and Sudarshan [21] in the context of quantum optics. Coherent states have also been widely used to construct IVRs guided by Gaussians in phase space, such as the HK propagator [42], the CCS representation [43], the linearised semiclassical IVRs [40] and the Frozen Gaussian Approximation [41]. For a detailed discussion and more coherent-state based IVRS see the recent review [44]. 3 For the Herman Kluk propagator, we consider a state vector |Ψ HK (t) = |q, p R(t, q 0 , p 0 ) q 0 , p 0 |Ψ (0) e iS cl dq 0 dp 0 2π , where |q, p represents a coherent state in phase space and q 0 , p 0 the initial phase space coordinates. This state corresponds to a Gaussian wave packet, see Eq. (13). The prefactor is given in terms of the elements m uv = ∂u/∂v 0 of the monodromy matrix, which is composed of the derivatives of the final phase space variables with regard to their initial values, and is the semiclassical action, where H cl (p, q) is the classical Hamiltonian (see Eq. (3) in this work). However, for clarity, when using IVRs we employ a slightly different notation than in the rest of the paper. This is done in order to distinguish between phase space Gaussian coherent states and their position-space representations (see Eq. (34) below). For a Gaussian initial wave packet, For the CCS method, we write the time-dependent wave function as a superposition of timedependent, nonorthogonal Gaussian coherent states (CS) |z = |z(t) , defined aŝ a|z = z|z and z|â † = z|z * , whereâ † ,â are the creation and annihilation operators, respectively, and whose eigenvalues are parametrised as functions of the phase space coordinate as The time-dependent state reads where denotes the classical action along the trajectory defined with regard to the matrix element H ord (z * , z) = z|Ĥ ord (â † ,â)|z . This represents the diagonal elements of the ordered Hamiltonian matrixĤ ord (â † ,â). In general, z|Ĥ|z = z|z H ord (z * , z ), where |z , |z denotes two arbitrary coherent states. In coordinate space, is a Gaussian wave packet centred at the phase space coordinates q and p. One should note that, for a static field and a Gaussian potential given by Eq. (5), the effective CCS Hamiltonian is given by with η = (λγ/(γ + λ)). Physically, averaging the Hamiltonian over a Gaussian coherent-state basis effectively lowers the potential barrier by introducing an effective, shallower potential. In comparison with the classical Hamiltonian, Eq. (35) shows an effective energy shift γ/4 [198,199].

Identifying bound and continuum regions
A powerful reason to study phase space dynamics in attoscience is that trajectories can easily be understood as bound or unbound depending on the phase space configuration. This is crucial when looking at signatures of strong field tunneling and over-the-barrier ionisation. Bound and continuum regions can be defined by inspecting the phase portrait of the system. This is exemplified with Fig. 1, in which we consider a model atom in a static field. Similar pictures have been provided in [73,200]. The figure shows two fixed points according to the definition provided in Sec. 3.1.1: a centre at the origin, and the Stark saddle to the left, whose positions are determined by the minimum and the maximum of the effective potential V eff (x), respectively. The figure also shows a separatrix, which corresponds to the stable and unstable manifolds of the Stark saddle. It is associated with the minimum energy to undergo over-the-barrier ionisation. The closed region to the right of the saddle is bound: trajectories within it will propagate along closed orbits. If the particle starts on the left of the Stark saddle, or if it has an energy higher than that of the separatrix, it will be free. This clearly shows that the particle's energy is paramount to defining the continuum regions; the particle being close to the core is not sufficient for describing its dynamics. Quantum mechanically, there will always be an uncertainty for the initial wave packet, which will manifest itself as a phase space spread (see Fig. 1, upper right panel). In trajectory-based grid methods, such as IVRs, one may employ an ensemble of classical trajectories to mimic the initial spread (see Fig. 1, lower left panel). One should note that, due to the initial uncertainty, some trajectories belonging to the ensemble will be in the continuum from the start. Although they cannot cross separatrices, they do from a tail whose energy is above that of the saddle. Energies above the saddle indicate above-the-barrier ionisation (see Fig. 1, lower right panel) [76].
For a time-dependent field of the form (7), in the parameter range considered in this work, namely low frequencies and high driving-field intensities, one may assume that the system behaves quasi-statically. This means that one may use the approximation that the phase space configurations discussed for static fields and shown in Fig. 1 hold for each instant of time. They will change instantaneously, such that an electron reaching the continuum at different times and propagating in the field will be exposed to a wide range of transient bound and continuum regions. This implies that the time should be included as an additional variable in an extended phase space.
The Stark saddles will be symmetric with regard to the origin for consecutive half cycles. An illustration of a time-dependent separatrix for a field given by Eq. (7) is displayed in Fig. 2. The bound region enclosed by the separatrix increases with time in the first quarter cycle of the field, as it evolves from a maximum to a crossing. This evolution is indicated up to t = 30 a.u. (see upper color bar), which is less than a quarter of a cycle. Thereafter, sign of the field will reverse and the bound region will decrease until the minimum of the field, in the subsequent half cycle, is reached.
The picture also shows a sample phase space trajectory, in which an electron performs two revolutions around the core until it eventually escapes. This orbit is indicated by the dashed line and should be followed in the clockwise direction, starting from the black dot. The electron is initially in the continuum, that is, outside the region enclosed by the separatrix at t = 0 (blue curve). It follows this separatrix closely, but is captured as the bound region increases. At t = 20 a.u. (solid square), it is in a bound region far away from the saddle indicated by the dark yellow curve. This region increases up to a quarter of a cycle. Thereafter, the oscillating field reverses its sign, and the electron performs a second revolution around the core, which, in phase space, is roughly the mirror image of the first with regard to the momentum axis. Eventually, the electron escapes as it is once more in the continuum (blue triangle).

Inhomogeneous fields and phase space configurations
Time-dependent separatrices and transient phase space configurations play a key role for inhomogeneous fields, which are employed to model HHG in alternative media such as nanostructures [82]. The laser field becomes enhanced by plasmonic resonances and the external field has to be regarded as inhomogeneous. Whilst it is well known in the literature that even small inhomogeneities may  From [82] lead to huge changes in the high-order harmonic spectra [164][165][166][167][168][169][170], such as an increased cutoff energy (see Fig. 3), the analysis of phase space regions provides further insight. An inhomogeneous field approximated by Eq. (8) introduces an additional fixed point with regard to its homogeneous counterpart, whose nature changes with the sign of the laser field. For E(t) > 0 or E(t) < 0, where E(t) is given by Eq. (7), it will be a centre or a saddle, respectively. The two configurations introduce either a concavity, see Fig. 4(a),(c) or a convexity, see Fig. 4(b), (d), to the dynamical system and are a crucial tool in understanding the increase in the HHG cutoff energy, previously attributed to the higher electron momentum for the convex system at the instant of ionisation [166].
However, by separating our configurations into two toy potentials (one always convex, one always concave), we show that, instead, the increased cutoff energies are due to the concavity providing [Lower row] High-order harmonic spectra computed using the dipole acceleration (23) for the same parameters as Fig. 3 using two symmetric toy models which we artificially constructed from [82]. Panel (c) corresponds to the two-saddle potential and panel (d) to the two-centre potential. The gray dots give the spectra computed for the homogeneous case. From [82].
additional confinement and forcing high-energy orbits back to the core. This is a similar mechanism to that in [201], in which an additional confining potential led to a substantial increase in the cutoff energy without loss of intensity. Indeed, Fig. 4(c),(d), shows that only the spectrum obtained from a concave potential leads to an increase in cutoff energy.

Enhanced ionisation and nested separatrices
The change in the nature of the fixed points due to the inhomogeneous field leads to different configurations, depending on the external laser field. This is a clear case of a bifurcation. Studying the H + 2 molecular system in [78] using phase space diagrams also leads to several fixed points, shown in Fig. 5. Thereby, one may identify a Stark saddle to the left, a central saddle between the two molecular wells and two molecular centres separated by R. We again have two different configurations, not due to a time dependent field, but to a bifurcation for increasing internuclear separation R.  (16) as a function of the inter-nuclear distance R, using different starting wave packets: delocalised (red), localised upfield (orange) and localised downfield (purple). The vertical line indicates the inter-nuclear separation for which the phase space configuration changes. From [78].
For R below the critical (bifurcation) value, see Fig. 5(a), the separatrices are nested and classical trajectories close to the core stay bounded. If the energy of the central saddle becomes greater than that of the Stark saddle, the separatrices open up, see Fig. 5(b). Then, a trajectory in the upfield molecular well only needs to pass the upfield potential barrier in order to reach the continuum. This explains the behaviour of the ionisation rate with respect to the internuclear distance, shown in Fig. 5(c). Indeed, after the bifurcation and as the separatrices open up there is a sharp increase in the ionisation rate. Additionaly, the ionisation rate for a wave packet localised around the upfield molecular well, Eq. (13), is about double that of a delocalised wave packet, Eq. (14), and the ionisation rate for a wave packet localised downfield, Eq. (13), is suppressed. From this it is deduced that maximum enhancement is achieved if the energy of the saddle is high enough for the tunnelling electron coming from the upfield population not to be trapped by the downfield centre, but low enough for the effective potential barrier to be narrower than that of a single atom.

Semi-classical versus quantum regimes
In the following, we will illustrate different physical regimes, which may or may not have a classical counterpart. Thereby, we will also exemplify roles that trajectory-ensembles may play. Classical-trajectory ensembles allow for an intuitive understanding of an electronic wave packet's evolution when paired with quantum mechanical methods. They may, for instance, be employed to link the outcome of quantum-mechanical computations, such as the TDSE, to the picture of a returning electron, or be used to construct grids for initial-value representations. A widespread example is to use trajectories to infer electron return times from HHG spectra. This is known since the 1990s [202], and the fact that these times are well specified within a field cycle has paved the ground for attosecond-pulse generation. It is also well known that windowed Fourier transforms of the time-dependent dipole acceleration [Eq.(23)] lead to periodic arch-like structures which can be explained using classical-trajectory ensembles and give pairs of return times merging at the cutoff [179,181].

Pairing trajectories with quantum mechanical methods
Some of these structures are illustrated in Fig. 6, in which the return times of an ensemble of classical trajectories are used to understand the cause of plasmonically enhanced HHG. While this is a quantum-mechanical process, when the trajectories are paired with time-frequency (Gabor) maps computed from the TDSE, they can offer very useful insight. There is very good agreement between the two, and we discern peculiar features, which increase with the inhomogeneity of the field. First, there is a suppression in the arch-like structures that occurs with the same periodicity as the field and happens for times after each field crossing. This can be explained by the phase space dynamics in Fig. 4: For times prior to the crossing, the prevalent configuration, with two centres, forces the electron to return, while the phase space configurations subsequent to it, with two saddles, may contribute to irreversible ionisation. This confirms that confinement plays a more important role than the electron reaching the continuum with a higher velocity. A second feature are structures occurring over several field cycles that extend up to very high frequencies. These structures are associated with a second time scale introduced by the additional centre, which will be briefly discussed in Sec. 3.4. The frequency with which they occur increases with the inhomogeneity parameter β. Further details are provided in our original publication [82].
In addition to that, classical-trajectory ensembles are employed to construct time-dependent grids which guide initial-value representations (IVRs). In [76], we assess the suitability of IVRs for modeling tunnel ionisation in an unusual setting, namely starting from a bound state. We perform these studies in phase space using the Wigner quasiprobability distribution constructed from IVRs and from the TDSE. Therein, we consider the Herman Kluk propagator, which is a semiclassical IVR, and the Coupled Coherent States (CCS) method, which is a quantum IVR solving the TDSE in a coupled coherent state basis. The two Hamiltonians are stated in Sec. 3.1.5 and differ, because the quantum corrections due to the coherent-state averaging introduce energy shifts and effectively lower the potential barrier. This implies that CCS orbits may cross classical separatrices, while the orbits employed in the HK propagator may not. A detailed discussion of the differences between the two propagators is provided in [198,199] and in our previous work [76]. Fig. 7 illustrates how specific trajectories in position and momentum space differ when using the two IVRs. An important issue is that the quantum corrections will allow the CCS orbits to cross classical separatrices. This leads to the CCS being more accurate at reproducing tunneling, and there will be a degradation of the outcome of the HK propagator for longer times. This degradation has been observed in our previous work [76] for high-harmonic spectra.
The agreement between the different methods is illustrated in Fig. 8, in which Wigner quasiprobability distributions were calculated for a Gaussian potential using the CCS and the HK propagators. In all cases, the Wigner function presents the signature semiclassical tail discussed in [73] (see also [200] for a seminal discussion of such features in the context of quantum localisation). The tail follows the separatrix and crosses from the bound to continuum region around the Stark saddle. There are also interference fringes on the left-hand side of the saddle, which may be associated with For the HK propagator and the CCS method we use 10 7 and 1600 trajectories, respectively. The color bars give the square of the Wigner quasiprobability density. From [76].
ionisation events in previous times. The agreement between the CCS and the TDSE is much better, with practically identical quasiprobability distributions. Furthermore, significant fewer trajectories are required than in the HK case. This is due to the quantum corrections in the CCS Hamiltonian H ord (p, q). Nonetheless, it is remarkable that these quantum features are present despite using HK wave packets and a real-trajectory grid. This is in strong contrast to how the trajectories guiding the grid behave, exemplified in Fig. 1: the HK trajectories with over-the-barrier energy never cross phase space barriers, while the Wigner function does. This discrepancy is justified by the non-locality of the Wigner function near separatrices, and by the fact that, near the Stark saddle, the barrier is approximately parabolic [135]. This means that, at least for short times, the HK propagator is applicable. One should note, however, that the HK trajectories contribute to forming the tail. Indeed, if the classical trajectories with initial positions outside the bound region are removed, the tail as well the interference fringes are absent [76]. This is consistent with the findings of [204], who show that the total weight of the classical phase space trajectories corresponding to energies below, or above the top of a parabolic barrier give the reflection or transmission coefficients.

Non-classical behaviour in phase space
In the examples provided above, one clearly sees a classical-quantum correlation, either in structures from time-frequency maps being related to classical return times, or in Wigner function tails following classical separatrices. However, what about situations in which the quantum phase space evolution has no classical counterpart? Below we illustrate the difference between semiclassical and quantum pathways, in the context of molecular enhanced ionisation [78].
The Wigner function nicely illustrates the physical mechanisms behind molecular enhanced ionisation, and these different types of behavior. From this quasiprobability distribution analysis in Fig. 9(a ),(b ), the signature semiclassical escape tail associated with over-the barrier ionisation is identified, as well as the oscillatory behaviour around the classical separatrix. However, intriguing structures that cycle through the momentum space, whose behavior does not follow separatrices, are also present. In Fig. 9(a)(b), the quasiprobability distribution starts a cycle from the upfield centre to the downfield well. Using the interference fringes around the central saddle as a quantum bridge, this flow does not follow any of the classical separatrices. As shown in Fig. 9(a ),(b ), the  cyclical movement then transfers part of the downfield population back to the upfield molecular well, while the semiclassical tail starts to form. In [74,75], those structures were called momentum gates and were related to the non-adiabatic following of the time dependent field. However, as is seen in Fig. 9, these are present even when using a static field, and for various internuclear distances (before and after the bifurcation in Fig. 5). This provides evidence that the actual mechanism is intrinsic to the molecule, instead of a response to an external field. In [78], we show that quantum interference builds a bridge between the two centres in the molecule, which supports this direct intra-molecular population flow. For that reason, we refer to the mechanism behind the momentum gates as "quantum bridges".
From studying the effect of different internuclear distances as well as different initial wave packet configurations, the optimal configuration for a static field stems from using a localised upfield initial wave packet, see Fig. 10(a). Indeed, the cyclical motion is absent, meaning the subsequent quantum bridge bringing the population back to the upfield well does not form. Therefore, the initial quantum bridge forms, leading the upfield population through the downfield centre, to the semiclassical escape pathway and to the continuum.
Understanding the formation of the quantum bridges seems paramount to optimising enhanced ionisation. Unfortunately, classical arguments are not sufficient in explaining the time evolution of the Wigner function during molecular enhanced ionisation. They fail to predict the cyclical motion of the quantum bridge after the bifurcation as well as its frequency. This all seems to point to that the evolution of the quantum bridge is inherently nonclassical. In order to quantify this we use the quantum Liouville equation, Eq. (19) to define the amount of quantum corrections to the evolution of the Wigner function. The result is quite staggering and quantum corrections along with the corresponding Wigner quasiprobability distribution are shown in Fig. 10. If the Wigner function has a fully classical time evolution, Q(x, p, t) vanishes everywhere. As shown in Fig. 10 (b), the quantum corrections become very strong around the quantum bridge as it starts to build up. In contrast, they are completely absent along the semiclassical tail. For this we conclude that enhanced ionisation is due to the interplay of two pathways: The semiclassical escape pathways associated with tunneling mechanisms that follows the separatrix with a classical Wigner function evolution (despite its description of an inherently quantum mechanical process) and the quantum bridge. The latter stems from the interference between the two molecular wells, breaks all phase space constraints, presents very strong quantum corrections and cannot have a classical analogue.

Understanding time scales
Phase space methods are also helpful for identifying time scales and the physics behind them. This is extremely important when working with time-dependent fields. For instance, in molecular enhanced ionisation, it is paramount to understand to timescale of formation of the quantum bridge as well as its cyclical motion, and how it relates to the frequency of the time-dependent field. Indeed, the quantum bridge can both aid or hinder ionisation and by changing parameters such as the type of field, the initial wave packet or the internuclear distance, enhanced ionisation can be optimised. As seen in [78], the quantum bridge and its frequency is inherent to the molecular system, and is present even in the absence of an external field. In [139] an analytical method based upon those employed in quasi-solvable models [205][206][207] is used to determine the origin and the value of the frequency of the quantum bridge in a field-free system. Using the autocorrelation function of different initial wave packets, it is clear that this frequency is due to the coupling of different eigenstates. This leads to a very interesting situation in [78], since the frequency of the quantum bridge is higher than the frequency of the time-dependent laser field. The quantum bridges could be controlled with the appropriate coherent superposition of states or different driving fields, which opens a wide range of possibilities for studying quantum effects in enhanced ionisation. We have verified that, for the model potential in [139], the frequencies computed analytically are quite robust upon inclusion of an external static field. This is illustrated in Fig. 11 4 . Fig. 11 (a) Comparison between the absolute value of the autocorrelation function |a(t)| 2 , see Eq. (15), calculated using the analytical method in [139] in a field free system (red, dotted line) and numerical computations from [78] using the same parameters but with a static field of strength E 0 = 0.01 a.u. (I = 3.51 × 10 12 W/cm 2 ) (dark blue solid line) and E 0 = 0.05 a.u. (I = 1.72 × 10 14 W/cm 2 )(light blue solid line). [Bottom row] Wigner quasiprobability distributions using the same parameters as (a) at time t = 0.7 a.u. using in panel (b) the analytical method in [139] for a field free system and in panel (c) the numerical method in [78] for a static field of strength E 0 = 0.05 a.u. (I = 1.7 × 10 14 W/cm 2 ).
Another example are the several timescales that occur for HHG in inhomogeneous media [82]. The very good agreement between the quantum time-frequency maps and the classical-trajectory computations means that we can explore the dynamical aspects of the latter. If the atomic potential V sc (x) is neglected, the equation of motion describing the trajectory ensemble in the inhomogeneous field can be re-written as Mathieu's equation. Explicitly, where Υ = βx + 1, τ = ωt, = βE 0 /ω 2 . This equation is widely used to study particles in ion traps [208][209][210][211], and provides a wide range of dynamic information, from the values of the inhomogeneity parameter β for which the system is stable, to the time scales involved. From the phase space study of individual trajectories, shown in Fig. 12, one can see that they experience two differing motions: a slow and large oscillation and a small and rapid one. This allows us to apply Dehmelt's approximation within the stability region, and its accuracy is shown in Fig. 12. From this we conclude that the secular oscillations are indeed responsible for the high-frequency structures in the time-resolved spectra seen in Fig. 6.

Conclusions
The take home message of the present review is that quantum optics, quantum information, chemical physics and the theory of dynamical systems have developed powerful toolkits that are underused in strong-field and attosecond physics, among them classical and quantum phase space. We have provided a few examples of how phase space arguments and/or quantum quasiprobability distributions can be employed in the context of high-order harmonic generation and strong-field ionisation, be it for establishing constraints and determining different dynamical regimes, or for studying non-classical effects. Moreover, the phase space can also provide guidance for constructing effective Hamiltonians, determining relevant subspaces and understanding the interplay between the residual binding potentials and the laser field in greater depth. This is particularly important in the context of correlated multielectron dynamics, for which the large number of degrees of freedom may pose additional difficulties. Some of these techniques have been referred to in the traditional setting of laser-induced nonsequential double ionisation (see Sec. 2.4), and will become increasingly necessary for extended systems such as large molecules, solids and nanostructures. For instance, being able to select the relevant degrees of freedom and treating them quantum mechanically, while describing the less relevant ones classically, or incorporating quantum corrections around classical evolution are widespread strategies in quantum chemistry [33,34], cold gases [32,35] and in recent years photosynthetic compounds [212,213]. Thereby, a key question is how to adapt these techniques to a highly transient, subfemtosecond regime, in the context of attochemistry. Recently, quantum and classical approaches have been combined to investigate electron and nuclear dynamics in pumpprobe experiments in glycine, and the initial state was computed using phase space techniques [214]. Even in the single-electron regime, it is desirable to move from one-dimensional models. This interest ranges from a more accurate description of the dynamics [162,163] to studies in orthogonally polarised fields [215].
In addition to that, it is clear that one must go beyond traditional modelling in strong-field laser-matter interaction, which employs either pure quantum states or classical methods, and often takes the system to be initially in the ground state. One then assumes that the ionising electron creates a hole in the orbital with the lowest binding energy. Nonetheless, one must bear in mind that, in real-life situations, the atomic and molecular ions generated by a strong laser field will be in a coherent superposition of states. This means that it is important to establish whether there will be coherences between the ionisation channels [216]. Furthermore, the outgoing electron is expected to exhibit a degree of entanglement with its parent ion, which may harm coherence. This has been discussed in recent theoretical work, in which two time delayed XUV pulses are employed to control the degree of entanglement between vibrational and electronic degrees of freedom in H 2 . Entanglement prevents coherent superpositions of states, which would lead to vibrational wave packets, from forming [217]. To be able to determine and prepare the ion in appropriate coherent superposition is really important in the context of attosecond hole migration [218,219] and pumpprobe schemes [214,217,220], which rely on well-defined phase relationships. Electronic coherences also play a key role in XUV induced bond formation [221]. Thereby, one must assess how nuclear and electronic degrees of freedom couple, with the aim of maximising coherence [222]. This also implies that a density-matrix formalism and effective Hamiltonians [216,223], which are more suitable for open quantum systems, are required. This is particularly true if one takes into consideration the current trends, towards extended systems such as large molecules [214,[224][225][226][227][228] and nanostructures [229,230], for which overcoming decoherence, quantifying entanglement and non-classical behavior, and controlling the coupling with the environment will pose major challenges. There is theoretical evidence that nuclear degrees of freedom in a large molecule weaken coherence. Still, it is remarkable that even in the worst-case scenario phase relations may survive [214]. The present review is a brief illustration of how phase space tools widely used in other research areas, such as quasiprobability densities, may be employed in attosecond physics.