Floquet Hamiltonian approach for dynamics in short and intense laser pulses

We present a time-dependent Floquet method for short pulses and arbitrary laser frequencies that uses the cycle-averaged Kramers–Henneberger basis. By means of a particular plane-wave expansion we arrive at a time-dependent Schrödinger equation governed by a Floquet Hamiltonian, which consists of convolutions of momentum and Fourier components. A dedicated numerical treatment of these convolutions, based on Toeplitz matrices and fast Fourier transformations, allows for an efficient time-propagation of large Floquet expansions. Three illustrative cases of ionization with different photon energies are analyzed, where the envelope of a short and intense pulse is crucial to the underlying dynamics.


Introduction
Non-perturbative laser-matter interaction provides a rich yet challenging area for theoretical studies. While numerical methods have to deal with large energy bandwidths required to fully account for the dynamics, analytical methods are faced with the challenge of finding an appropriate description of non-perturbative light-matter interaction.
A successful analytical approach to non-perturbative laser-matter interaction is the Kramers-Henneberger (KH) approximation [1,2]. It describes the dynamics in the KH reference frame co-moving with the laser-driven electron(s). Traditionally, the time-dependence of the problem is eliminated by using a Hamiltonian, averaged over one optical cycle. The corresponding potential and eigenenergies are often referred to as the KH approach or the KH 'atom.' For sufficiently large field strengths and high frequencies, the cycle-averaged Hamiltonian () largely determines the properties of the coupled light-matter system, while all higher-order corrections remain small and can be treated perturbatively.
Since its introduction, the KH approach was thoroughly examined and the properties of  are very well known [3,4]. It was applied to a large variety of problems in atomic, molecular [5] and solid state physics [6,7] and in particular, was indispensable in the study of ionization suppression phenomena for atoms in strong and high-frequency fields. Nevertheless, most of the theoretical predictions were not tested experimentally (see [8][9][10] for application for Rydberg state ionization) because high-intensity and high-frequency lasers were not available at that time. Evidence of the existence of 'KH atoms' were found recently in experiments using low-frequency laser fields [11,12].
The situation has, however, changed due to the free-electron lasers (FEL) [13] that are already able to provide pulses of sufficiently high-frequency and intensity to enable the observation of non-perturbative phenomena. The first experimental studies of Raman processes in the VUV and XUV frequency range, which require coherent multiple photon absorption/ emission, have been carried out [14]. It can be expected that FELs will soon reveal high-frequency non-perturbative phenomena which were proposed theoretically, such as adiabatic stabilization [3], dynamic interferences [15][16][17][18][19], Rabi oscillations between core-hole states [20], to name a few.
The KH approach is ideally suited to describe strong-field high-frequency physics to be realized in FEL facilities apart from one crucial aspect: studies so far were mostly limited to continuous-wave laser radiation. Indeed, for a continuous-wave field, a perturbative expansion into Floquet states can be readily developed. On the other hand, for short (FEL) pulses Floquet theory cannot be straightforwardly applied. Clearly, the timedependent aspect is crucial since the short pulses created by FEL sources can lead to additional dynamics driven by the pulse envelope as was recently predicted [21,22], or be necessary to account for phenomena like impulsive Raman scattering [23]. Hence, in order to apply the KH approach to the dynamics involving intense and short pulses, a formulation different from the ones so far known appears to be necessary.
Here we propose a numerical approach for short-pulse non-perturbative laser-matter interaction that is based on a time-dependent Floquet formalism in the KH reference frame. It uses ( ) t which depends on the instantaneous intensity of the laser pulse and relies on time-propagation using the full Floquet Hamiltonian, which is performed with an efficient fast-Fourier-transformation (FFT) based algorithm. Combining these two approaches allows us to obtain both a qualitative and quantitative understanding of the light-matter interaction during the laser pulse, despite treating short laser pulses produced by FEL facilities non-perturbatively.
In section 2 we will present the time-dependent Floquet Hamiltonian approach, followed in section 3 by the introduction of the novel algorithm to solve the Floquet problem in momentum space. The approach is illustrated in section 4, where the role of the envelope of a short and intense laser pulse is investigated for the ionization in 1D potential. By varying the laser frequency, while keeping ( ) t invariant, three parameter ranges are explored: high, intermediate, and low-frequency regimes. We show in section 4.2 that ( ) t provides an excellent approximation of the laser-driven dynamics for frequencies higher than the binding energy of the potential. For intermediate frequencies close to the ionization threshold, discussed in section 4.3, the pulse envelope plays a crucial role in determining the channels involved in the ionization. Finally, the low-frequency regime is discussed in section 4.4; in this case, the photon energy is much smaller than the binding energy of the field-free potential and several hundred Floquet channels are required to fully account for the dynamics. We show that the population is rapidly distributed over many excited states of ( ) t during the rising part of the laser pulse, which has to be considered if one wants to use the KH approach for low-frequency fields.

Time-dependent Kramers-Henneberger-Floquet approach in momentum representation
The time-dependent Floquet Hamiltonian approach we will formulate in the following is a generalization of the Envelope Hamiltonian introduced in [21] geared towards efficient numerical implementation. The formalism allows one to explore the transformation of the wave function from the field-free to the 'field-dressed' picture while still fully accounting for the effects of a finiteness of a (short) laser pulse.

Kramers-Henneberger transformation
The time-dependent Schrödinger equation (TDSE) within the single active electron approximation in the KH reference frame acquires the form [1,2,24] (in the following, atomic units will be used, unless stated otherwise) where the coupling with the laser field is reduced to the timedependent shift a( ) t of the binding potential a . For simplicity we assume this shift to be of the form corresponding to the classical trajectory of a charged particle in a laser field with vector potential ( ) t A . The KH transformation describes the laser-atom interaction in a frame of reference, where the electron can be considered to be 'stationary,' while the binding potential of the 'atom' is time-dependent. In other words, the electron 'sees the nucleus oscillating back and forth'. The oscillating potential a , assuming a fixed a 0 , can be integrated over a single cycle T ω of the oscillation a( ) t to obtain the averaged potential which is also called the 'KH potential' and the corresponding cycle-averaged Hamiltonian which has been used to describe the properties of atoms in strong and high-frequency laser fields [5,25]. The average potential strongly depends on the electron excursion length a 0 , as illustrated in figure 1, and, for sufficiently large a 0 , transforms from a single-well to a double-well shape. The cycleaveraged potential still depending on the time-variation a ( ) t 0 of the laser pulse envelope, will be central in our exact timedependent Floquet Hamiltonian formulation of the TDSE.

Time-dependent Floquet Hamiltonian approach for short laser pulses
The time-independent KH potential and its properties are analyzed in great detail in the literature using a variety of methods [24][25][26][27][28][29][30][31][32][33]. In practice, however, one needs to deal with finite and often short pulses and considering a static KH potential is not sufficient. Here, we consider a cycle-averaged potential that adjusts to the laser pulse envelope, while still providing an exact description of the dynamics. At first glance it looks cumbersome to perform for each instance of time a full cycle average. However, when switching to momentum space one arrives at a compact and distinct form of the TDSE (see equation (8) below), as we will show briefly here and in detail in appendix A.
The potential in the KH reference frame can be written in a plane-wave expansion  with integer m, anticipating that the potential oscillates with frequency ω. In order to efficiently treat short pulses we will extract the envelope of the laser pulse. This is done by splitting the electron displacement a( ) t into the non-periodic envelope a ( ) t 0 and the periodic oscillation . Thereby, the KH potential becomes a 'two-time potential' which can be used straightaway in expansion (4) to give ò å a + =~w - The components  ( ) V t k, m depend on time through the pulse envelope and, by means of a translation in space and the Jacobi-Anger expansion, can be rewritten as products (see appendix A and [29] for the derivation) where  ( ) V k is the field-free potential and with J m denoting the ordinary Bessel functions of the 1st kind.
Having rewritten the potential as a sum of products we use a similar expansion for the wave function. This allows one to derive an equation for the kth plane-wave and mth Fourier component y  ( ) t k, m of the wave function (see appendix A for details) ò å a y wy The accuracy of its numerical implementation is limited only by the basis and propagation routines, see section 3.1 for more extended discussion. The momentum representation, although not frequently used for solving TDSEs, has some advantages over position space approaches, see e.g. [34].
Here it is used because it reduces (i) time-averaging of the KH potential to multiplication by appropriate factors as in equation (7) and (ii) the TDSE to a convenient form that allows one to use efficient numerical propagation methods, as described in section 3.
There are other approaches that adopt two times [35][36][37][38][39]. Here the two times are used to straightforwardly derive expansion (6), which turn out to be very convenient for short pulses. In contrast to the other approaches, we do not use a Floquet ansatz for the wavefunction or solve the TDSE for two time variables. Instead, we use the Floquet Hamiltonian defined by (8) and solve the TDSE for a general wavefunction. Nevertheless, the Floquet modes associated with the instantaneous field intensity can be obtained by calculating eigenstates of the Hamiltonian. This Floquet Hamiltonian method is similar to the one used for vibrational dynamics of the + H 2 molecule [40]. Our approach, however, utilizes the cycle-averaged potential, which is particularly suited in connection with intense laser fields and is accurate even for very short pulses.

Physical interpretation of the wave function in the Fourier basis
The physical significance of the index m becomes apparent, if we consider an isolated Fourier subspace m, i.e. ignore the coupling between the wave function coefficients y  ( ) t k, m with different m. In such a case, the only remaining potential coupling terms in (8)  describe the cycle-averaged potential. Therefore, considering a single Fourier subspace in isolation is similar to the original KH approach [1], where only the cycle-averaged potential is considered. The eigenstates of this potential within the Floquet theory correspond to the Floquet states in the infinitefrequency limit [41].
0 couple different Fourier subspaces and lead to transitions between the states of the cycle-averaged Hamiltonian ( ) t . The index m can be interpreted as the number of absorbed/emitted photons [42]. For example, population initially created in the m=0 subspace and ending up in the m=n subspace after the pulse represents n-photon absorption. We will refer to the nth Fourier subspace also as Floquet channel. Our numerical implementation will allow to include enough Floquet channels to achieve numerical convergence so that fields of arbitrary frequency can be considered.

Numerical implementation
To numerically solve the TDSE in (8) we first rewrite it for a discrete momentum k grid, yielding å a y wy where for D dimensions the field-free potential and the wave function are renormalized according to The right-hand side of (10) can be split into two parts. The first part, which using matrix notation is defined by is diagonal and can be easily computed numerically. The computation of the sum requires the main numerical effort as it is associated with the non-diagonal elements of V. In the field-free case (a = 0 0 ), the part (12) describes a convolution between momentum components k of the wave function and the potential. If the laser field is present (a ¹ 0 enter the sum (12). They couple different Floquet channels m and also modify the coupling between momentum components k. Nevertheless, the convolution form of the matrix V in (12) is preserved, since the couplings depend only on the differences -¢ k k and -¢ m m . Note, that V and y depend on time, which will be kept implicit for the brevity of notation.
The convolution form of the matrix V allows one to apply the convolution theorem and to replace the convolution between potential and wave function, described by (12), by their product in the Fourier domain. This greatly increases the speed of computation, in particular if a FFT algorithm is used to convert to and from the Fourier domain.
The convolution theorem strictly holds only for infinite or periodic vectors, which implies an expansion in k and m to infinite order. In practical numerical calculations, the necessity to use a finite size basis will normally violate the conditions for validity of the convolution theorem, consequently causing numerical errors. Therefore, we use an alternative approach that is based on the theory of Toeplitz matrices [43]. It takes advantage of the convolution form of the matrix V and allows one to use the FFT algorithm to accelerate the calculations. However, unlike the direct application of convolution theorem, the method based on the Toeplitz matrix theory is exact for vectors of finite size. This approach is particularly useful to study Floquet systems, as it allows one to truncate the basis to only a few Floquet channels.
1. Description of the algorithm For a single Floquet channel, e.g. = ¢ m m , the elements along the diagonal of the matrix V in (12) are equal, which follows directly from the momentum representation. Such a matrix is called a Toeplitz matrix and its properties are well known in the literature, see, e.g. [43]. It can be fully described by a single row and column only. Furthermore, a product of a finite Toeplitz matrix with any vector can be performed exactly using the FFT algorithm.
The algorithm to multiply a Toeplitz matrix V with a vector y is as follows (see appendix B.1 for a more detailed description):

A circulant vector c is formed from the first column and
the first row of the matrix V. 2. Zeros are appended to the vector y to match the length of c.

A Fourier transformation of both the circulant vector c
and the extended coefficient vector y is performed. 4. The two transformed vectors are multiplied and an inverse Fourier transformation is applied to the product.
The first half of the final vector now stores the matrix-vector multiplication result, while the second half is discarded. If the couplings between different Floquet channels are considered, i.e. -¢ ¹ m m 0, then the matrix V is of block form with all equal blocks on the same diagonal. Additionally, each block is of Toeplitz form. Such a matrix is called a Block Toeplitz matrix with Toeplitz Blocks (BTTB). The product of a BTTB matrix and any vector can be performed using a two-dimensional Fourier transformation algorithm in a similar way as a Toeplitz matrix-vector product, see appendix B.2 for a detailed description. The approach can be further extended to an arbitrary number of dimensions.
The algorithm to calculate the Toeplitz matrix-vector product can be considered as generalization of the wellknown split operator technique (transformation to Fourier domain, multiplication and inverse transformation), see, e.g. [44], which is widely used to solve the TDSE. On the other hand, the algorithm presented here cannot be used to directly evaluate the product of a vector with a function of Toeplitz matrix, e.g.
. Nevertheless, the Toeplitz matrix-vector multiplication algorithm allows us to reduce the number of Fourier and plane-wave components required to achieve high numerical accuracy and allows it to outperform the traditional split-operator technique.

Time propagation
Many different numerical methods to solve (10) could be used, for example explicit Runge-Kutta or Arnoldi-Krylov algorithms. However, to take advantage of the BTTB symmetry of the potential matrix V, the matrix-vector multiplications involving V must be implemented using efficient FFT routines with the method described above. In this work the Taylor expansion propagator is used. This method relies on the expansion of the propagator over a discrete time-step Dt in a Taylor series up to the desired order, so that the wave function expansion coefficients can be computed as where each term in the expansion is evaluated iteratively. Hence, the numerical problem reduces to the evaluation of products , which is evaluated using the Toeplitz matrix-vector multiplication algorithm presented above. The accuracy can be controlled by choosing the order of expansion at each time-step. Although the propagator is not norm-conserving, if enough expansion orders are included norm conservation up to a desired numerical accuracy can be easily achieved. In this work, the expansion was truncated once the norm of the corrections to the wave function coefficients dropped below 10 −16 . The Taylor expansion propagator, combined with the FFT algorithm for matrix-vector multiplication operations, allows to treat large Fourier expansion orders m explicitly. More sophisticated propagation methods that also rely on matrixvector products like Arnoldi-Krylov-propagators may be easily implemented.

Accuracy
The accuracy of the time-dependent Floquet Hamiltonian approach developed in this work is verified by comparing the wave function obtained by directly solving the TDSE in velocity gauge with the solution of the TDSE defined in (10). In both cases, identical plane-wave basis and propagator routines of the TDSE were used.
For all laser pulse parameters that were used in this work, the wave functions obtained from the time-dependent Floquet Hamiltonian approach and by directly solving the TDSE in velocity gauge were found to match up to numerical accuracy, if sufficiently many Fourier orders m were considered. The accuracy of the time-propagation procedure is determined by the time-step and Taylor expansion order. The maximum required number of Fourier components m max can be determined from the plane-wave basis set by requiring that where | | k max is maximum momenta described by the plane-wave basis. In practice, less Fourier components are sufficient.
An illustrative example is provided in figure 2 for ionization from a soft-core potential, which is defined in section 4.1, with ω=1 a.u. photon energy, =´-I 2.4 10 W cm 18 2 intensity and 5 fs full-width at half-maximum (FWHM) duration pulse. The spectra under similar laser pulse parameters were extensively investigated in previous works [15][16][17][18][19] and the calculation is further discussed in section 4.2, therefore here we only note that each Floquet channel provides the m-photon absorption channel, see figure 2(a). The final spectra, obtained by summation over all Floquet channels m, is indistinguishable from the spectra obtained by the direct solution of the TDSE in velocity gauge, see figure 2(b). Note, that the energy resolved spectrum is strongly modulated and does not show a smooth envelope, as might be expected from single-photon ionization. These modulations are the result of a large Stark shift of the initial state during the ionization process [18,19].
The approach was tested to be accurate for photon energies ranging from 0.05 to 1 a.u.. Furthermore, it was accurate for pulses down to single cycle duration for both low and high frequencies. Therefore, the time-dependent Floquet Hamiltonian formalism is capable of fully describing the dynamics driven by intense and short laser pulses using the cycle integrated Hamiltonian for arbitrary laser parameters.
The accuracy of the numerical procedure is further dictated by the quality of the plane-wave basis set. In all the calculations presented in this work, a converged basis set in terms of maximum momenta | | k max and spacing between momenta components Dk is used. Finally, note that atomic potentials with a long-range tail lead to a singularity at the origin in the momentum representation, i.e.
. This singularity could be treated by, for example, a Landè subtraction procedure [34,45]. Alternatively, the potential can be considered in a 'finite box', as in [46]. Finally, since adding a delta function to a momentum representation of a potential leads only to a trivial energy shift of the spectra, the singularity can be removed by introducing a new potential xd k k is the Kronecker delta function. Numerically this corresponds to setting all elements -¢=  V k k 0 to zero. Here we use the latter approach as it is both simple and provides accurate results for all observables studied in this work.

Performance
The Toeplitz matrix approach described above allows to efficiently solve the time-dependent Floquet Hamiltonian formulation of TDSE in (10) using the FFT based matrixvector multiplication. The use of FFT algorithm allows to achieve scaling proportional to ( ) N N log with respect to the size of the basis N=N K ×N F , where N K is the number of plane-waves and N F is the number of Fourier components. This scaling is illustrated as a function of the total number of basis elements N in figure 3 for different calculations with laser pulse parameters used in this work. The size of the plane-waves basis is varied between 512 and 4096 with the maximum momenta kept fixed. The number of Fourier components is varied between 1 and 401. The numerical effort is measured in terms of processor cycles spent solving the TDSE. The number of cycles is then divided by the total number of time-steps used, so that calculations using different pulse lengths could be directly compared. In addition, the expansion order of the Taylor propagator in (13) is kept fixed at 8. In an adaptive expansion scheme, the expansion order mainly depends on´D | | t k max . The main effort required to solve the TDSE stems from updating the wave function at each time-step using the Taylor expansion method, which is illustrated by the black dots in figure 3. It scales proportionally to ( ) N N log as expected. The numerical effort required for different sizes of the plane-wave basis is shown in figure 4. Again, the effort scales proportionally to ( ) N N log F F with the number of Fourier components. Additional numerical effort is required to update the elements of potential energy operator at every time step, since they depend on the laser field. This effort is illustrated by green crosses in figure 3. It can take up to 40% of the total effort. However, it scales linearly with the total number of basis elements since onlyŃ N K F elements are stored in memory.
The time-dependent Floquet Hamiltonian method presented here cannot compete in efficiency with conventional approaches to solve TDSE that do not use Floquet expansion.   The numerical effort required for the latter would be comparable to using just a single Floquet channel, see figure 4. However, the time-dependent Floquet Hamiltonian method does not aim to compete with the established approaches in terms of speed or accuracy, but rather to provide an efficient way to tackle large-scale Floquet problems in for short laser pulses. Therefore, the strength of the current approach is its ability to provide insight into the dynamics during the laser pulse, which is possible due to Floquet-like approach only.
1. Extension to more spatial dimensions Although this work is limited to one-dimensional potentials, the generalization to more dimensions D for a plane-wave basis in Cartesian coordinates is straightforward. However, such an approach does not take advantage of the symmetry of the potential and therefore in general requires a large number of basis elements to be included into the Hamiltonian, which scales asŃ N K D F . The corresponding increase of numerical effort can be extrapolated from figures 3 and 4.
A direct extension of the method to, e.g. a spherical coordinate system is not straightforward. The advantage of the plane-wave basis in Cartesian coordinates is the separation of any arbitrary potential in the KH reference frame into time-independent and time-dependent parts, as can be seen from (7), which allows us to calculate the coupling between the plane-wave components at each timestep efficiently. We did not find such a simple form for the expansion of the KH potential into spherical harmonics for linearly polarized fields.
A possible alternative approach to describe atoms in linearly polarized laser field beyond a single dimension is to use cylindrical coordinate system (see, e.g. [47]). Since the KH potential is symmetric around the laser polarization axis, a plane-wave expansion can be applied along this direction. The KH approach can also be formulated for a circularly polarized field, see, e.g. [47]. Finally, multi-pole expansion of the KH potential can be used, which allows for an efficient description using conventional quantum chemistry methods [31].

Model system
The time-dependent Floquet Hamiltonian approach is illustrated using the example of a 1D model atom, described by a soft-core potential which has been widely used to study the dynamics of atoms in high intensity laser fields both analytically and numerically [27,48]. In this work the softening parameter is chosen to be = x 2 0 2 , which leads to a binding energy equal to that of a hydrogen atom I p =0.5a.u.. The laser pulse with a peak electric field F 0 is defined in terms of the classical electron trajectory introduced in (2) with a Gaussian envelope function [21] a w w An envelope of T=5 fs FWHM is used, unless specified otherwise. For all except few cycle pulses w( ) P T 1 holds.
Furthermore, laser intensity and frequency are chosen such that the maximum classical excursion length a 0 is equal to 10 a.u. 16 0 0 2 for all laser frequencies investigated. Since the KH potential depends only on the classical excursion length a 0 the eigenenergies of the cycle-averaged potential will have identical time-dependence. Nevertheless, the dynamics will still depend on the frequency via the spacing between Floquet channels. Therefore, the choice of a constant maximal a 0 will allow one to clearly separate the role of the cycle-averaged potential from the role of the laser frequency.
The typical eigenenergy spectra of the 1D cycle-averaged potential a ( ) V x, 0 0 are depicted in figure 5(a) as a function of time during the pulse. They are obtained by diagonalizing the Hamiltonian in a single Floquet channel. The eigenenergies of ( ) t strongly depend on the instantaneous intensity of the laser pulse due to the widening and the formation of the dichotomy of the cycle-averaged potential in figure 1 for increasing electron excursion a 0 .
In figure 5(b) the effective quantum numbers * = n -E 0.5 n are plotted for bound states of ( ) t as a function of time along the laser pulse. Bound states of a hydrogenic potential would lead to an infinite series of equally spaced n * , which is the case at early times in figure 5(b). Deviation from a pure hydrogenic potential lead to an uneven spacing of n * , referred to as quantum defect. In the case of the cycle-averaged potential, the quantum defect is a result of the dichotomy of the potential and is clearly visible for the lowest eigenstates. On the other hand, although the energy of the higher (n>3) eigenstates is lowered due to the widening of the cycle-averaged potential, n * stay approximately equidistant, indicating that these states are determined by the long-range Coulomb tail of the cycle-averaged potential and are not strongly influenced by its dichotomy.
All the calculations presented further in this work were done using 2048 plane-wave basis states with momenta equidistantly spaced by p D = k 2 2000 a.u.. A time-step of Δ t=0.1 a.u. was used for the propagation, which we found sufficient to obtain converged ionization probabilities. The ground state was obtained by diagonalizing the Hamiltonian defined by (10) with a = 0 0 for a single m=0 Floquet channel. Finally, the carrier-envelope phase was set to f = 0 in all calculations.

High-frequency
The KH approach was originally proposed in the context of high-laser frequencies, for which the underlying dynamics is by now mostly well understood, see [3,4] for comprehensive reviews. Therefore, the high-frequency case provides a good reference point to illustrate the influence of the pulse envelope on the dynamics induced by high-intensity laser fields using the time-dependent Floquet Hamiltonian approach developed here. We choose a laser frequency of w =1 a.u. 27 eV, substantially larger than the field-free ionization potential. Accordingly, the peak laser intensity is set to =´-I 2.4 10 W cm 18 2 , so that the maximum electron excursion length is a = = ( ) t 0 10 0 a.u.. Note that for these laser parameters, non-dipole effects contribute negligibly to the dynamics [49]. Hence, we safely work in the dipole approximation.
The population in the Floquet channels m=0 and 1 as a function of time is depicted in figure 6(a). During the initial part of the pulse around 30% of the population is transferred from m=0 to the m=1 Floquet channel, indicating a onephoton absorption process. The population transfer stops, when the adiabatic stabilization regime is reached. Around the peak of the pulse, despite the rapid increase of field strength, population in each Floquet channel m stays approximately constant, implying that the Floquet channels are decoupled as predicted by the high-frequency Floquet theory [50]. As the field intensity decreases at the end of the pulse, another 20% of the population is transferred to the m=1 Floquet channel by single photon absorption. The remaining Floquet channels (m>1) contain =1% of the population after the pulse.

Non-adiabatic excitations
Projecting the population in each Floquet channel m onto the eigenstates of ( ) t reveals that the ionization process is adiabatic, see figure 6(b). The population in the m=0 subspace stays in the ground state throughout the dynamics and no substantial transitions due to the time-dependence of ( ) t takes place. In this case, the ionization process can be described in terms of a discrete state that belongs to the m=0 Floquet channel embedded into the continuum of  (14)) as a function of time along the pulse envelope (bottom axis) and classical excursion length a 0 (top axis) for a maximal excursion length of a = 10 0 a.u. and a pulse defined in (15).
states that belong to the m=1 channel, as assumed within the high-frequency approximation [50]. Furthermore, once adiabatic stabilization sets in, the envelope of the laser pulse plays a minor role.
The adiabatic picture is not applicable for a shorter laser pulse with the same peak intensity. In this case, due to the rapid change of the eigenstates, non-adiabatic excitations from the ground to the excited states occur as is shown in figure 7 for a T=1fs FWHM pulse. However, the population stays in the m=0 Floquet channels, i.e. no photons are absorbed from the field indicating that the excitations are induced by the envelope of the pulse. This is confirmed by excitations of even-parity states only, as absorption of a photon would lead to the excitation of odd-parity states. At the end of the pulse, the excited states of ( ) t transform to the corresponding field-free states. Such non-adiabatic transitions [51] were investigated in [21] using the envelope Hamiltonian formalism, where it was shown that they can be quantified using time-dependent perturbation theory.

Intermediate-frequency
A range of 'intermediate' frequencies can be defined, where the laser frequency is smaller than the binding energy of the field-free potential, but larger than the binding energy of the cycle-averaged potential at peak intensity. We will show, that for such 'intermediate' frequencies, the field-free ground state does not simply adiabatically connect to the ground state of ( ) t , as in the high-frequency case. Instead, it undergoes a series of crossings with excited states that belong to higher Floquet channels.
We choose a laser frequency of ω=0.4 a.u.∼10.9 eV and an intensity of =´-I 9 10 W cm 16 2 . Therefore, two photons are required for ionization, however the photon energy is still twice the binding energy of the cycle-averaged potential at the peak of the pulse, see figure 5(a).
The time-dependent populations in the m=0, 1 and 2 Floquet channels, which are shown in figure 8(a), immediately suggest that ionization proceeds in a sequential manner via the intermediate m=1 channel. This is confirmed by the time-dependent population in the eigenstates of ( ) t depicted in figure 8(b). While the ground state of the m=0 Floquet channel is rapidly depopulated, the population is transferred to the odd-parity excited states of the m=1 channel, i.e. via one-photon transition. From these states, the population is slowly transferred to the m=2 channel, i.e. ionized via absorption of a further photon.   figure 8(c). At the beginning of the pulse the energy of the ground states rapidly increases and undergoes a series of crossings with the excited states that belong to the m=1 Floquet channel. At each of these crossings, a fraction of ground state population is transferred to the excited state. As the energy of the ground state decreases at the end of the pulse, the population is exchanged again at the second crossing. Between these crossings, a small but significant coupling of the states leads to small Rabi oscillations that are seen around t=0 in figure 8(b).
The two transitions that occur at the crossings of the ground and excited states of ( ) t lead to interference that depends on the phase accumulated in each state in between. This phase in turn depends on both the energy differences and the couplings between the states. Since the time between two crossings depends on the pulse duration it strongly influences the final population after the pulse, as seen in figure 9(a).
The origin of oscillations in figure 9(a) is further elucidated by the evolution of excited-state populations during the pulse. In figures 9(c) and (d) these populations are shown for pulse durations, when either n=0 or n=3 state is predominantly populated after the pulse. Initially the dynamics in both cases is very similar. Clear differences emerge only after the second crossing between the states, indicating that interference effects determine the final populations.
The final populations oscillate with well-defined frequencies as the pulse duration changes, see figure 9(b). The biggest amplitude oscillation is between the ground and n=3 state, which is to be expected since the coupling between these states is at least a factor of two larger than between any other states and n=3 state is the first one to undergo a crossing with the ground state. The time-window of strong interaction with the ground state is also longest for the n=3 state, since the crossing occurs at the beginning of the pulse, where the energy-time gradient is not as steep as for higher-energy states.
Unlike the final population of each state in figure 9(a), which requires one to consider all interactions, we found that the frequencies of oscillation of final populations in figure 9(b) are determined mainly by the dynamics of the ground and a single excited state. The presence of higherenergy states does not significantly perturb these frequencies, since they mainly depend on the phase difference accumulated between the times of crossing of the two states. These times in turn depend on their energy difference-the higher the energy of the excited state, the later the first crossing will occur. Higher-energy states contribute to smaller frequencies reducing their influence.
Interaction of the ground state with any individual excited state can be readily described by a Landau-Zener-Stückelberg (LZS) interference process [52,53]. In case of multiple states with non-trivial time-dependence of energies and couplings, the interconnected LZS transitions lead to complicated and rich dynamics. Nevertheless, characteristic features prevail. The ground state will be depleted sequentially transferring population to higher excited states at later times. Therefore later crossings will become less important due to weaker couplings and the smaller population available for transfer. Hence, the traces of single state dynamics show up in figure 9(b) even for very high laser intensities.

Low frequencies
Although the KH approach was originally proposed to study the interaction of atoms with high-frequency laser fields, it was speculated that it could also be applicable for low-frequency and high intensity radiation [32,33,54,55]. More recently, the KH approach and in particular the properties of the cycle-averaged potential was used to explain the nonlinear Kerr effect in laser filamentation [56] and acceleration of neutral atoms in laser fields [57]. The time-dependent Floquet Hamiltonian formalism developed here allows one to directly investigate the KH approach for frequencies that are much smaller than the ionization potential.
We choose the laser frequency of ω=0.057 a.u., which corresponds to λ=800 nm wavelength radiation, and I= 3.7×10 13 W cm −2 intensity, so that the maximal electron excursion length is again a = 10 0 a.u.. Therefore, the eigenenergies of ( ) t and its dependence on the pulse shape shown in figure 5 is identical to the high and intermediatefrequency cases analyzed above. The FWHM duration of the pulses is set to T=30fs. However, the essential results presented here do not depend sensitively on the duration of the pulse. Note that in order to obtain converged results, 201 Floquet channels (±100 ω) are treated explicitly in the numerical calculation. The population in Floquet channels as a function of time is plotted in figure 10. Almost all of the population is transferred to higher Floquet channels at the peak of the pulse and then returns to the m=0 channel at the end of the pulse. Therefore, these transitions are virtual, which is in contrast to high and intermediate-frequency cases. This is not surprising, however, since for the given laser field parameters the total ionization is less than 1% and most of the population is expected to stay in the ground state after the pulse.
Projecting the population at the peak of the pulse onto the instantaneous eigenstates of ( ) t in figure 11(a) reveals that the population is distributed over many states. Crucially, no single state is dominating. Also, many Floquet channels are populated during the peak of the pulse, as is shown in figure 11(b). An increase of the field intensity leads to a broadening of the distribution over both the excited states n and also over Floquet channels m.
Virtual excitations created in multiple Floquet channels can be understood qualitatively within the Floquet picture. Since the interaction strength between states that belong to different channels is much larger than the energy spacing between them, a quasi-continuum of states is created. In this situation, the eigenstates of ( ) t do not correspond to any adiabatic or nearly adiabatic states of the field driven system. Therefore, the wave function, which stays nearly identical to that of the field-free ground state, is distributed over many excited states in the KH reference frame. The redistribution occurs at avoided crossings between states that belong to different Floquet channels, similarly as in the intermediatefrequency case. However, for low frequencies many more avoided crossings become important.
For sufficiently large peak field strengths a regime may exist, where ( ) t becomes applicable [32,54,55]. However, by the time this intensity of the laser pulse is reached, the wave function is already distributed over the excited states of ( ) t . Therefore, in order to apply the KH approach at low laser frequencies, it is essential to consider the transformation of the field-free ground state wave function to the 'fielddressed' KH picture during the switching-on of the pulse.

Summary and conclusions
We have developed a time-dependent Floquet Hamiltonian approach formulated in the KH reference frame to study dynamics driven by short and intense laser pulses. It constitutes a systematic and flexible extension of the Envelope Hamiltonian [21] applicable for arbitrary frequencies and provide a convenient and efficient way to propagate Floquet Hamiltonians. Numerical application of Floquet approaches is often hampered by the rapid increase of the number of Fourier components required to describe the Hamiltonian. Therefore, we have devised an efficient numerical procedure to propagate the Floquet Hamiltonian which is able to overcome the hurdle of large expansions. Indeed, we have performed calculations for laser parameters, for which several hundred Floquet channels had to be considered explicitly. Key element is the formulation of the problem in the momentum representation, which is particularly suited for the KH reference frame as it allows us to separate the components of the fieldfree potential from the field-dependent ones. We further use the formalism of Toeplitz matrices and the FFT algorithm to achieve a favorable scaling with Floquet channels. However, unlike the split-operator methods that also rely on FFT, the Toeplitz approach is numerically exact for finite-size matrices. For Floquet problems it allows us to truncate the basis to only several Floquet channels. Yet, the method can be applied to any other Fourier basis to achieve an efficient and accurate propagation.
The main advantage of the method is its ability to extend the KH approach, which is particularly suited for high-frequency and high-intensity fields, to the limit of very short pulses. Thereby, we can investigate physical effects that emerge at high intensities and can only be understood by explicitly considering the time-evolution of the pulse envelope.
We have shown that the pulse envelope exerts control over two types of dynamics. For very short pulses, the rapid change of the eigenstates of ( ) t over time leads to nonadiabatic excitations. They are induced by the pulse envelope and therefore do not involve the absorption of any photons from the field. Thus, even for very high frequencies they can lead to significant population in low lying bound states.
The second type of transitions, sensitive to the pulse envelope, occurs at the crossings between the discrete eigenstates of ( ) t that belong to different Floquet channels, as their energy changes along the pulse. Although dynamics at each crossing can be easily understood in terms of LZS theory, in our case multiple states are strongly coupled evading simple interpretations. Nevertheless, strong features due to the coupling between individual states can be discerned, which is quite remarkable, considering the high intensities used. They lead to a large sensitivity of final state populations to the pulse duration providing a possible route for their coherent control.
An extreme case for our KH approach is the low-frequency limit, when the photon energy is much smaller than the binding energy of the electron. In this case, the population is transferred between the bound states of ( ) t at their crossings, which are very dense due to the small energy spacing between the Floquet channels, resulting in a rapid distribution of the population over many eigenstates of ( ) t before any populations has a chance to reach continuum states.
To summarize, the time-dependent Floquet Hamiltonian approach presented here provides a convenient basis for short laser pulses and for all but the smallest photon energies. In all investigated cases, including few-cycle pulses, the approach was able to provide accurate numerical results indistinguishable from the ones obtained using conventional techniques of propagating the TDSE. However, unlike the conventional TDSE propagators, the time-dependent Floquet Hamiltonian method allows one to obtain insight into dynamics that crucially depend on the pulse envelope. Such dynamics will become particularly important for short and intense pulses generated by FEL facilities, which often devise unusual pulse shapes, or for investigating the dynamics of states dressed by unusually shaped pulses such as used in [12].
The richest envelope-dependent dynamics is observed in the intermediate-frequency range. For multi-electron systems, this energy range will be much more extended due to the ubiquitous presence of core and double excitations, which lead to a rich energy structure even for very high photon energies. Therefore, we expect that for such systems the timedependent Floquet Hamiltonian formalism presented here will be even more valuable. was used. Note that d ¢w ( ) t t is periodic, with the period T ω =2π/ω. However, its use is justified since , , and the integration in (A5) can be limited to the range ¢ Î where zero was appended in each row to concatenate a Toeplitz matrix row with its column, and an arbitrary number of zeros can be used when forming a circulant matrix. Multiplication of a circulant matrix C and a vectorx is equal to a convolution between the first column of the circulant matrix º C c n0 and the vectorx. Hence, from the convolution theorem it follows that       where å denotes the convolution operation and  is the discrete Fourier transformation.
To efficiently multiply a vector x by a general Toeplitz matrix T we need to: (i) form a column of the circulant matrix c from the 1st column and the 1st row of the Toeplitz matrix T; (ii) append the vector x with zeros forming an extended vectorx that has the same length as the vector c; (iii) perform the Discrete Fourier Transformation of the vectors c andx, multiply them together element-by-element and perform the inverse discrete Fourier transformation of the product. The result will be stored in the first N elements, where N is the size of the original vector.