Exploring laser-driven quantum phenomena from a time-frequency analysis perspective: A comprehensive study

Time-frequency (TF) analysis is a powerful tool for exploring ultrafast dynamics in atoms and molecules. While some TF methods have demonstrated their usefulness and potential in several of quantum systems, a systematic comparison among these methods is still lacking. To this end, we compare a series of classical and contemporary TF methods by taking hydrogen atom in a strong laser field as a benchmark. In addition, several TF methods such as Cohen class distribution other than the Wigner-Ville distribution, reassignment methods, and the empirical mode decomposition method are first introduced to exploration of ultrafast dynamics. Among these TF methods, the synchrosqueezing transform successfully illustrates the physical mechanisms in the multiphoton ionization regime and in the tunneling ionization regime. Furthermore, an empirical procedure to analyze an unknown complicated quantum system is provided, indicating the versatility of TF analysis as a new viable venue for exploring quantum dynamics.


Introduction
Time-dependent quantum mechanics is a fundamental topic in physics, chemistry, and engineering. Traditionally, dynamics of a quantum system such as lifetime and energy difference between states can be revealed by spectral lines in the frequency domain by using the Fourier analysis. However, the Fourier transform presents limited chronological information of a dynamical system. To explore chronological information in a quantum dynamical system, time-frequency (TF) methods are necessary.
In recent experiment advances, TF methods such as wavelet transform combined with ultrafast spectroscopy techniques were used to probe dynamics in molecules, solids and liquids [1,2,3,4,5,6,7], particular in light harvesting complexes [8,9,10]. In attosecond physics, the short-time Fourier transform (STFT), the continuous wavelet transform (CWT), Wigner-Ville distribution (WVD) (which belongs to the Cohen class distribution), were employed to uncover dynamical mechanisms including high-order harmonic generation (HHG) [11,12,13,14,15,16,17]. Lately, Hilbert transform and the synchrosqueezing transform (SST) were applied to investigate multiple scattering [18] and the dynamic origin of near-and below-threshold harmonic generation, respectively [19,20]. While these studies show the usefulness and potential of each individual TF analysis for probing dynamics of a quantum system, several crucial issues remain when it comes to analysis of a quantum system with unknown complicated dynamics. First, different types of TF methods and the choice of window functions may provide extremely different TF representation, leading to conflicting physical interpretation. Second, it is difficult to select a particular TF method based on previous studies, in which the physical models are different, e.g., 1D [12,18], 2D [12] and 3D atoms [12,15,16,19,20], and solved by different numerical schemes [12,15,20,18]. Clearly, TF representations derived from different TF methods based on different physical models cannot provide an impartial comparison. Third, in the past decades, several modern TF methods, e.g., the Cohen class distribution, empirical mode decomposition with Hilbert spectrum (EMD-HS), reassignment methods (RM), and SST have been proposed and successfully applied to classical macroscopic dynamical systems including molecular dynamics [21], cardiopulmonary coupling phenomena [22], chronotaxic systems [23], lamb wave propagation [24] and seismic data [25]. However, Cohen class distribution other than the WVD, EMD-HS, RM, and different forms of SST have not been discussed in quantum dynamical systems. As a consequence, a comprehensive TF study in the same benchmark is essential.
To this end, we take 3D hydrogen atom as a benchmark because hydrogen is one of the most representative systems in quantum mechanics. All simulations are based on the same numerical method (time-dependent generalized pseudospectral method), and different contemporary TF analysis, including the STFT, CWT, Cohen class distributions, RM, SSTs and EMD-HS, are applied to study the simulated signals. The systematic comparison of the TF methods enables to interpret physical processes in a quantum system. This article is organized as follows. In Section 2, we summarize the physical model of the hydrogen system in a strong laser field and the numerical simulation. In Section 3, we provide an overview of the TF methods, as well as a discussion of their pros and cons. In order to show that quantum dynamical phenomena can be depicted by the TF methods, we perform simulations with two different sets of laser parameters, including hydrogen in the regime of multiphoton ionization and tunneling ionization. The results of TF representations by different TF methods and a discussion of their limitations and potential applications are given in Section 4. In Section 5 we conclude the paper and provide an empirical procedure to analyze an unknown complicated quantum system.

Theoretical modeling
We simulate hydrogen dynamics in a linear polarized laser field in the framework of the electric dipole approximation and non-relativistic quantum mechanics. Note that for the wavelength range of approximately 500 to 1200 nm, the electric dipole approximation is valid [26] only when the laser intensity I 0 is smaller than 10 15 ∼ 10 16 W/cm 2 . The time-dependent Schrödinger equation for atomic hydrogen interacting with a linearly polarized field along the z-axis in atomic units can be expressed as where ψ(r, t) is the wave function at position r = (x, y, z) and at time t, z = r cos(θ), and E(t) is the external laser field.
To numerically solve this equation, we adopt a time-dependent generalized pseudospectral method [27,28] which consists of two essential steps: (1) The spatial coordinates are optimally discretized in a nonuniform fashion by means of the generalized pseudospectral technique: the grid is denser near the origin and sparser away from the origin; (2) A second order split operator technique in the energy representation, which allows the explicit limitation of undesirable fast-oscillating high energy components, is used to obtain an efficient and accurate time propagation of the wave function.
According to previous studies [29,30], dynamical phenomena such as the HHG is associated with the electric dipoles, i.e., the laser-driven electron oscillating around the stationary nucleus, which could be expressed as, respectively, the time-dependent induced dipole in the length form, denoted as d L (t), and in the acceleration form, denoted as d A (t) [28]: (2) d L (t) = ψ * (r, t)zψ(r, t)dr.

Time-frequency Analysis
Suppose the signal x(t) is composed of finite oscillatory components, that is, x(t) = L l=1 a l (t) cos(2πφ l (t)), and each component has a time-varying amplitude modulation (AM), a l (t) > 0, and time-varying instantaneous frequency (IF) φ l (t), where stands for the first derivative, then the TF representation reflects the IF, which is a generalization of the notation frequency, and the localized phase information could be extracted. We call a l (t) cos(2πφ l (t)) an intrinsic mode type (IMT) function. More details about this kind of function could be found in Appendix A.
To disclose the time-varying nature of this kind of signal, in particular the IF and AM, the TF analysis is a powerful approach. In the following subsections we provide an overview on classical and contemporary TF methods.
3.1. Linear TF methods. The global analysis nature of the Fourier transform is responsible for its limitation in extracting time-varying dynamics inside an oscillatory signal. An intuitive way to resolve this limitation is analyzing the signal locally; that is, we could crop a segment of finite length and apply the Fourier transform, and expect to observe how a signal oscillates locally. This idea leads to the STFT. We introduce a window function g to crop the signal at different time and perform the Fourier analysis. This analysis results in a time-frequency (TF) representation of the signal [31]: We call STFT g x (t, ω) a TF representation of the given function x and |STFT g x (t, ω)| 2 the spectrogram of x. Here the STFT by Eq. (4) is different from the conventional definition by an additional modulation factor e −iωt .
Nevertheless, temporal and frequency resolutions cannot be achieved simultaneously by the STFT, according to the Heisenberg uncertainty principal [31]. In other words, a wide window provides a good frequency estimation in the STFT at the cost of poor temporal resolution, while a narrow window has the opposite trade-off. In this research we follow the tradition and choose the Gaussian function with a standard deviation of σ (the full width at half maximum of the Gaussian is 2 √ 2ln2σ) as the window function, i.e., g(ζ) = e −ζ 2 2σ 2 , which leads to the Gabor transform (GT): Here the GT we apply here differs from the usual one by an additional modulation factor e iωt .
Based on the same local analysis idea as that of the STFT, we could analyze the momentary behavior of the signal x(t) by the CWT: where g is the chosen mother wavelet and the scale parameter a controls the dilation of the window, and * denote the complex conjugate [31,32]. Due to the dilation nature of the transform, the TF representation analyzed by the CWT has a good time resolution and poor frequency resolution at the high frequency region and a good frequency resolution and poor time resolution at the low frequency region. One of the commonly applied mother wavelet is the Morlet wavelet, and the CWT based on the Morlet wavelet is called the Morlet wavelet transform (MWT): where 2τ 2 e iζ is the Morlet wavelet and τ > 0. The standard deviation for the dilated Morlet wavelet is σ = τ /ω. In Eq. (7), we follow the convention and use ω, which is the inverse of the scale parameter a in Eq. (6). Note the fundamental difference between Eq. (7) and Eq. (5) -in the GT, the Gaussian function is of a fixed width but in the MWT the width varies. In the MWT the window width σ varies with ω such that the quality factor τ (the inverse of the relative bandwidth) is invariant on the TF plane [31]. In other words, the window width σ becomes smaller as ω increases, and vice versa. Since the frequency resolution depends on the scale, we say that the TF representation by CWT has an adaptive frequency resolution.
Both the STFT and the CWT belong to the linear type TF analysis, in which the signal is characterized by their inner products with a preassigned family of templates with free parameters. Clearly, these TF representations depend on the chosen window, which might cause artificial patterns on the analysis result. Further, while in general the underlying structure of the signal under analysis is not known, there is no systematic way to design the window which faithfully reflects the structure. The above facts render the linear type TF analysis non-adaptive to the signal under analysis.

Quadratic TF methods.
To resolve the non-adaptivity issue of the linear type TF analysis as well as finding the higher order structure, the Wigner-Ville distribution (WVD) [31] was proposed. The WVD is based on the concept of the autocorrelation function and is defined as Note that by a slight change of variable, we could view the WVD as a variant STFT which is free of the window choice issue. In this sense, WVD is adaptive to the signal under analysis. Also note that by a direct derivation, the spectrogram could be understood as follows [31]: The WVD has several good mathematical properties, e.g., the signal energy is preserved in the WVD, the WVD is a real-valued function on the TF-plane and it provides a precise information about the chirp signal, like x(t) = e i2π(β0+β1t+β2t 2 ) , where β 0 ∈ R, β 1 > 0 and β 2 ∈ R. However, choosing the signal itself as the window function causes specific artifacts depending on the signal type. For example, when there are more than one oscillatory component in the signal, the interference patterns in the TF representation is inevitable, which might lead to mis-interpretation of the signal.
By construction, the WVD is quadratic in the signal which could be viewed as an energy distribution. A direct generalization of the WVD based on imposing some constraints on the covariance structure [31] of the signal leads to the Cohen's class. In other words, the Cohen's class , which comprises all bilinear TF representations that are covariant to shifts in both time and frequency [31] and has the following form: where f (ξ, ζ) is an arbitrary parameter function. Different parameter functions lead to different TF representations. Note that the WVD is a member of the Cohen's class when f (ξ, ζ) = 1. Different parameter functions lead to different TF techniques.
One particular technique is eliminating the interferences in the WVD when there are more than one oscillatory component. To be more specific, we can employ a separable parameter function f (ξ, ζ) = G(ξ)h(ζ) ≡ W h (ξ, ζ) in Eq. (10), where suitably chosen G(ξ) and h(ζ) permit a continuous and independent control of the interferences in time and frequency, respectively. The corresponding representation is called the smoothed pseudo-Wigner-Ville distribution (SPWVD) [31] and could be rewritten as where g(ζ) is the inverse Fourier transform of G(ξ). Particularly, G(ξ) and h(ζ) are even functions with g(0) = 1 and H(0) = 1, where H(ξ) is the Fourier transform of h(ζ).
We mention that we could consider a variation of the Cohen's class to get the family of the affine class [31], e.g., the affine WVD or the affine SPWVD. The transform considered in the affine class is the time-scale covariant and the TF representation also has the adaptive resolution feature. To simplify the discussion, we do not study the affine class transforms, but we mention that their performance is similar to those in the Cohen's Class. The above-mentioned transforms are overall called the quadratic TF analysis. A higher-order generalization is also possible, and we refer the reader of interest to, for example, [33,34] 3.3. Reassignment Methods. As features of a TF representation computed by a linear type of TF method are smeared by the introduced window functions and those by the quadratic type ones are obscured by interferences, the RM was proposed to sharpen the resolution of the TF representation. Generally speaking, the coefficients of the TF representation at (t, ω) is reallocated to a different point (t,ω) according to a predefined reallocation rule [35,36,37]. A common choice of the reallocation rule is to assign values of a TF representation to the local centroids. Note that this is different from the averaging idea behind the STFT (Eq. (9)). In this study we apply the RM to the STFT and SPWVD, and called the methods RM-STFT and RM-SPWVD, respectively.
The reallocation rule for the RM-STFT is derived by estimating the local centers of gravity, denoted as: The notations tg and dg stand for the window function tg(t) and the first derivative of g(t) in the STFT, respectively. The RM-STFT is thus defined as Essentially, this formula reassign the energetic contents of the spectrogram to the new location (t x ,ω x ). As a consequence, the RM leads to a substantially improved resolution in the TF representation. Note that the reassignment rules Eq. (13) and Eq. (14) can lead to the group delay and the IF of the bandpass filtered signal y(t) = STFT g x (t, ω) [37], using only the unwrapped phase of the STFT g x (t, ω). However these physical meaningful expressions are numerically inefficient. Other forms of reassignment operators can be found in [31].
Next we consider SPWVD. Although the SPWVD can smooth out the interferences in the WVD, the smoothing functions introduce artificial broadening. Applying the reassignment technique on the TF representation of the SPWVD can reduce such artifacts. Similar to Eq. (15), the reassigned representation of the SPWVD is where the reassignment rule (t x ,ω x ) is determined by the concept of expectation: Despite RMs are intuitive techniques to sharpen the linear type and quadratic type TF representations, inverse routines as well as mode reconstruction are not available.

Synchrosqueezing transform.
In this section we describe the SST, which is a special RM aiming to address the intrinsic blurring issue in the linear TF methods. To be more precise, the SST manifests IF characteristic according to the reallocation rule that consists of solely the frequency information, rather than the centroid of the TF representation [38,39]. In addition to sharpening the TF representation, the causality property of the signal is preserved in the SST, which allows the decomposition of oscillatory components when the signal is composed of several oscillatory components. We refer the reader who has interest in the theoretical analysis and other details, like reconstructing each oscillatory component and robust to noise, of SST to [38,39,40]. In this section we describe the synchrosqueezed STFT (SST-STFT) and the synchrosqueezed CWT (SST-CWT).
is the Fourier transform of g and d is the smallest gap of IFs of any two consecutive oscillatory components. The SST-STFT is given as whereω x (t, η) is the reallocation rule and α > 0 is a controllable smoothing parameter for the resolution, which in practice is chosen to be small. The reallocation rule given by the following equation utilizes the phase information hidden inside the smeared TF representation: Note that only the frequency reassignment operator is considered in the SST-STFT, so that the causality property of the signal can be preserved. A slight modification ofω(t, η) based on the second-order information of the IF [41] might further improve the TF representation. Indeed, when the interested oscillatory component could be approximated by a chirp function [41,36], although the TF representation is sharpened by the SST-STFT, a mild spreading is inevitable. We call this mild spreading the diffusive pattern. To cope with the diffusive pattern, we could consider the following reassignment rule: The reconstruction formulas [39] for each IMT function x k (t) from SST g,STFT where 1, denotes taking the real part, and C g ≡ g(0).

3.4.2.
SST-CWT. The definition of the SST-CWT is similar to that of the SST-STFT. For a CWT with a mother wavelet g such that supp(G(ω)) ⊂ [1 − ∆, 1 + ∆], with ∆ < d 1+d , the SST-CWT is given by where α > 0 is a smoothing parameter andω x (t, η) is the reallocation rule defined bŷ When describing an oscillatory component with a fast varying IF, the TF representation in SST-CWT may show a slight spreading. While the second-order information of the IF has been discussed in the SST-STFT [41], it is not yet been studied.
The reconstruction formula [39,40] for each IMT function x k (t) under SST-CWT is where 1, and R g := η dη. The results of the SST are "adaptive" to the signal in the sense that the error in the estimation depends only on the first three moments of the mother wavelet instead of the profile of the mother wavelet [38,39,40]. In other words, the influence of the chosen window on the associated TF representation is minimized compared with the linear TF methods. In addition to the above properties, the SST is robust to several different kinds of noise, which might be slightly non-stationary [40]. An important property shared by the SST-STFT, SST-CWT and RM-STFT is that by taking a short window, the fast varying IF could be well captured.
At the first glance, the results of the reassignment technique and SST seem to break the well-known Heisenberg uncertain principle, which says that the temporal and frequency resolution cannot be achieved simultaneously. However, in the reassignment technique and SST, we have shown that we could obtain a TF representation with an almost perfect time and frequency resolution. The main reason is that we focus on the oscillatory signals, in particular the adaptive harmonic model but not the whole L 2 space. We mention that the behavior of the RM and SST on the general L 2 function is an open question.
3.5. EMD. One popular way to define the IF of a given signal x(t) is via finding the analytic representation of x(t) with the aid of the Hilbert transform: where H(x (t)) denotes the Hilbert transform of x(t). In this equation, a(t) and Φ(t) are the modulus and phase ofx(t), respectively [42,43]. The IF of x(t) is thus defined as the rate of the varying phase: However, the IF provided by Eq. (28) is not always meaningful. To deal with this problem, the empirical mode decomposition algorithm (EMD) [44] was proposed to decompose x(t) into a series of functions called the intrinsic mode function (IMF), x = K k=1 x k (t), on which the Hilbert transform can be applied subsequently, via an algorithm called the sifting process. IMFs must satisfy the following two conditions: (i) In the whole data set of a signal, the number of extrema and the number of zero crossings must either be equal or differ at most by one; (narrow band) (ii) At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero (adoption of local properties). Note that in general an IMF is different from an IMT function. Details of the EMD algorithm could be found in, for example, [44].
For each IMF x k (t), the corresponding IF and AM amplitudes can be estimated by Eq. (27) via the Hilbert transform and Eq. (28). By assigning the IF and amplitudes of all IMFs on a TF plane, we obtain the Hilbert spectrum (HS) HS x (t, ω) (HS), which is a TF representation for x(t) determined by where δ(ω) denotes the Dirac delta measure. Although the EMD along with HS has been applied to several fields, due to a number of heuristic and ad hoc elements in the EMD algorithm, it is difficult to analyze its accuracy and limitation. In addition, despite the solid mathematical support, the Hilbert transform might be limited when applied to analyze the momentary dynamics of an oscillatory signal. First, due to the slow decay nature of the kernel in the Hilbert transform, keeping the causality of the signal structure might be difficult. Second, when the signal has time-varying amplitude and frequency, in general there is no guarantee to get the correct analytic signal from a given real oscillatory signal by the Hilbert transform in Eq. (27).
We mention that the part of using the Hilbert transform to obtain the AM and IF of each IMF could be replaced by using the SST-STFT or SST-CWT. In this case, the method is regarded as the EMD-SST-STFT and EMD-SST-CWT.
A comparison of different TF methods discussed in this study is summarized in Table 1.

Results and Discussions
In this section we employ the aforementioned TF methods on the time-dependent dipole computed by the time-dependent generalized pseudospectral method at the ab initio level. In a strong laser field, a variety of dynamical processes can occur in atoms and molecules, such as ponderomotive effect, the AC Stark effect, HHG [46], multiphoton ionization [47,48,49], above threshold ionization, tunneling ionization, rescattering of electron wavepacket [50,51], etc. Up to date, the STFT [11,12,13], CWT [14,16,15], WVD [17], and Hilbert transform [18] have been engaged in investigation of estimation of emission time and multiple The multiresolution-like behavior of the sifting process was studied in [45]. ‡: The stopping criteria of the sifting process depends on tuning several parameters. : The inversion from the HS is not possible for the definition (29). It is possible if we define the HS as k a k (t)e iΦ k (t) δ(ω − ω k (t)). * : It is sensitive to noise; a shortly existing oscillatory component destroys the whole analysis.
scattering. However, to the best of our understanding there is no comprehensive study comparing these TF methods in a quantum dynamical system.
By taking the well-studied atomic hydrogen as a benchmark, here we focus on the laser-driven hydrogen in the multiphoton ionization regime and in the tunneling ionization regime. When the laser photon energy ω is much smaller than the ionization potential I p , the multiphoton ionization and the tunneling ionization can be described by a dimensionless Keldysh parameter γ K = I p /2U P , where U P = E 2 0 /4ω 2 0 is the ponderomotive potential [29,48]. Similar physical dynamics can be achieved for any atom-field interaction given a fixed γ K [52]. Generally speaking, multiphoton ionization of atoms become dominant when γ K 1, while tunneling ionization is predominant when γ K 1 [29]. In the following we perform simulations with laser parameters in these two regimes.

Multiphoton Ionization Regime.
In the first simulation, the laser field parameters are arranged such that the Keldysh parameter is γ K = 3.07, indicating that multiphoton ionization is the major mechanism. The laser wavelength is 880 nm, which corresponds to ω 0 ≈ 0.05178273 in the atomic unit (a.u.), and the laser intensity of I 0 = 10 13 W/cm 2 , which corresponds to a field amplitude E 0 ≈ 0.0169 (a.u.). Note that for this γ K value, both tunneling and multiphoton ionization can occur, but the latter is dominant. The laser field ( Fig. 1(a)) E(t) = E 0 F (t) sin(ω 0 t) has a profile of F (t) = sin 2 (πt/(nT )), where n = 60 is the pulse length measured in the optical cycle (T = 2π/ω 0 ).
The computed induced dipole in the length form d L (t) is presented in Fig. 1(b). The power spectrum [28], computed by the squared Fourier spectrum, of the laser field and d L (t) are shown in Fig. 1(c)and 1(d), respectively. Both the length and the acceleration forms of dipole moment present the same detail structures in their power spectra [28]. Here we present the analysis using d L (t) for the multiphoton ionization case.
Although the profiles of d L (t) and the applied laser field look similar in the time domain, they are different in the frequency domain. For example, while the power spectrum of the laser field has only one peak located at ω 0 , that of d L (t) reveals odd harmonics due to the parity symmetry [29]. However, the meaning of the substructures within the odd harmonics and their corresponding dynamics are unclear.
To unveil the dynamics of d L (t), first we apply the conventional linear and quadratic TF methods, as shown in Fig. 2. The TF representation of the GT, computed by Eq. (5) with σ = 57.94 a.u., and the MWT, computed by Eq. (5) with τ = 30, are displayed in Fig. 2(a) and Fig. 2(b), respectively. Both the GT and MWT depict separate broad lines regarding the odd harmonics in the HHG process. While frequencies of the extracted components are consistent with the information shown in the Fourier spectrum in Fig. 1(d), the GT and MWT further capture the momentary behavior of each component. For a fixed resolution in the GT, low frequency components could be lost given an insufficient window width. In the TF representation of MWT, the same components are expected to be captured but the representation should be different from the GT due to the dilation nature of the transform; that is, the adaptive frequency resolution of the transform. Indeed, the frequency resolution below the 7th harmonic is improved compared with Fig. 2(b), while harmonics on the upper TF plane remain broaden. In addition, the frequency-dependent weighting in √ ω enhances the representation for high order harmonics. Note the width of the window function g(t) is uniform in the GT, and increases as ω increases in the MWT, as shown by the vertical lines between neighboring harmonics in the upper TF plane. The crescent shape distribution near the boundaries of the TF plane is caused by the boundary effect in CWT -for a low frequency, the corresponding scale is large and the influence of the window cut-off near the boundary becomes apparent. The boundary effect is also inevitable in the GT, like the vertical artificial pattern in the very beginning and end. However it is less dominant.
As discussed above, while the window is avoided in the WVD, a strong interference pattern is inevitable [31]. The TF representation of the WVD, shown in Fig. 2(c), reveals components between the odd harmonics, which is inconsistent with the features in the Fourier spectrum and violates the parity symmetry in physics. Moreover, harmonics above 15, of which intensities are weak, are not observed within the range of colorbar. Interference along the time-direction is also observed. As there are more oscillatory components and more overlaps caused by x(s + ζ/2)x * (s − ζ/2) in Eq. (8) as time increased, interference becomes stronger, as shown on the right-half TF plane in Fig. 2(c). Despite generating inaccurate information, note that in this simulation the WVD suggest the features of the TF representation, i.e., slow-varying components, regardless of a window parameter. This suggest that the window width in the linear type TF methods should be large enough such that patterns depict the harmonics as that performed by the WVD.
In order to remove the interference, we employ the SPWVD with the filtering functions g(s) = e −s 2 2σg 2 with σ g = 100 and H(ξ) = e −ξ 2 2σ H 2 with σ H = 0.025. Note that g(0) = 1 and H(0) = 1. The filtered result is given in Fig. 1(d). Although the even harmonics are eliminated, temporal interference cannot be fully removed. The introduced filters also worsen both temporal and frequency resolutions, as well as generating unexpected features. Note that the dynamic range (colorbar) is increased to reveal weak details.
By reassigning the energy distribution of the linear type TF methods and the SPWVD to their local centroids, the TF representations are sharpened. The results of reassigned GT (RM-GT) and reassigned SPWVD (RM-SPWVD) are shown in Fig. 3(a) and (b). In Fig. 3(a), the TF representation of the RM-GT significantly improved the resolution, revealing clear and distinct odd harmonics, as well as their subtle variation. Note that it is the square of the GT representation that is reassigned, therefore the range of the color bar differs from Fig. 2(a).
In Fig. 3(b), although the broadening that comes with filtering functions are removed and the TF representation depicts odd harmonics similar to those in Fig. 3(a), the interference artifact cannot be eradicated. We further employ the SST on the GT (SST-GT) and MWT (SST-MWT), in which the allocation rule is a good approximation of the IF of each IMT. The results shown in Fig. 3(c) and (d) also depict similar features as that in the RM-GT. Note that the SST-MWT preserves the adaptive resolution, weighted frequencies, and the boundary effect of the CWT.
In Fig. 3(a), (c) and (d), we show that the IMT functions can be decomposed without ambiguity by the RM-GT, SST-GT, and SST-MWT. Although the RM can illustrate the distinct odd harmonics and frequency shift of the AC Stark effect by finding the local centroids, there is no routine for reconstructing the harmonics. SST, on the other hand, comes equipped with an inverse method (See Eq. (23) and Eq. (26)).
To explore the physical meaning of the shifting in Fig. 4(a)-(c), we analyze it using the Floquet method [53], which has been extensively used in chemical physics [54,55]. Details of the Floquet method can be found in Appendix B. In Fig. 4(d), the blue, green, red, and cyan lines denote the energy difference of 1s-2s, 1s-2p x , 1s-2p z , and 1s-2p y , respectively, computed by the Floquet method, in the unit of ω 0 . Due to the breaking of the spherical symmetry by the laser field, the energy levels of 2s, 2p x , 2p y , and 2p z shift and split. The quasi-energies derived by the Floquet computation and the curve by the SST-GT are superimposed for comparison. We found that the spectral line on the TF representations obtained by the SSTs not only demonstrate the AC Stark effect but also depict the selection rule, i.e., only the 1s-2p z energy difference is presence.
It is worth mentioning that figures in Fig. 3 illustrate a frequency shift at the beginning cycles of the 7th harmonic, which could be the AC Stark effect. In the following context, we discuss about whether the IF components obtained by the reassigned values have a physical meaning, or they are simply clamped by artificial processing.
Details around the 7th harmonic are enlarged in Fig. 4 for the (a) RM-GT, (b) SST-GT and (c) SST-MWT. Note that the colorbar in Fig. 4(a) is different from from that in Fig. 4(b) and (c), since in the RM-GT the squared TF representation is reassigned. Despite different intensities, the RM-GT, SST-GT and SST-MWT depict a similar shifting trend descending from the 7.241-th harmonic to the 7.200-th harmonic, which corresponds to the AC Stark effect. The 7.241-th harmonic at the beginning cycles corresponds to the energies for 1s-2p transition ( 1 2 (1 − 1 2 2 )/ω 0 = 7.241 in a.u.), manifested as a small peak overlapping the H7 of the power spectrum in Fig. 1(d). The intensity of such spectral line is small because it arises from the near resonance absorption (not resonance absorption).   Fig. 4(d), we conclude that the decomposed IMTs are physically meaningful.
The ionization process for the laser field in Fig. 1(a) is summarized in Fig. 5. In Fig. 5(a), when the laser field is small, the spectral line 7.241 aroused from the 1s-2p near resonance absorption caused by the atomic structure. As the laser intensity gradually increases, the energy levels of 2s, 2p x , 2p y , and 2p z shift and split due to the breaking of the symmetry by the electric field, which is regarded as the AC Stark effect. When the laser intensity increases furthermore, the high-order harmonic process is induced and the dynamics is dominated by the transition between dressed states |M, N formed by the electron state M and the photon state N , as illustrated in Fig. 5(b). A summary for comparison of TF methods in the multiphoton ionization regime is provided in Table 2. Note that in this table we mention that the odd harmonics and the AC Stark effect cannot be obtained by EMD with the HS. In fact, we could not yet find physical meaning for the series of IMFs decomposed by the EMD in the multiphoton ionization regime.

4.2.
Tunneling Ionization Regime. In the second simulation, we discuss the TF representations for signals in the tunneling ionization regime. The laser wavelength in this subsection is ω 0 ≈ 0.05696100 in atomic units (a.u.), corresponding to 800 nm, and the laser intensity is I 0 = 3.5 × 10 14 W/cm 2 . The Keldysh parameter is γ K = 0.57, suggesting that the tunneling ionization is dominant. The laser field has 10 cycles    and its profile is ramped on according to which leads to the laser field E 0 F (t) sin(ω 0 t), as shown in Fig. 6(a). We adopt the dipole moment in acceleration form in the case of tunneling ionization, as displayed in Fig. 6(b). The acceleration form can provide clearer information for the high order harmonics than the length form, because the differential operator in Eq. (3) acts as a high-pass filter. TF representations for the linear type and quadratic type of transforms are presented in Fig. 7, including the GT (Fig. 7(a)), MWT ( Fig. 7(b)), WVD (Fig. 7(c)), and SPWVD ( Fig. 7(d)). In Fig. 7(c), we observe a periodic repetition of arches, suggesting the chirp-like dynamics of the attosecond radiation, as predicted by the standard semiclassical model [29]. To present such feature in the linear type TF representations, a short window width and a small quality factor τ are needed for the GT and the MWT, respectively. (Applicability of a short window width in the IMT models deserves further investigation.) Here we select σ = 1 (a.u.) for the GT and τ = 6 for the MWT. TF representations in the GT and the MWT illustrate that the arch structure repeats every 0.5T. The inner structure after the 2nd cycle is indistinct because of the broadening caused by the window. While the GT cannot resolve information below the 15th harmonic due to the small window width, the MWT can describe the near and the below threshold harmonics (low order harmonics) by its merit of adaptive resolution. WVD, on the other hand, can characterize the chirp very well without the need of window parameter and the broadening that comes along. However, as described in previous subsection, WVD introduces strong interferences in both temporal and frequency directions. Despite in SPWVD, where we use the parameters σ g = 1 and σ H = 0.13, it does not seem possible to eliminate the interference that mingled with intrinsic structure inside the main arches.
We further apply the RM techniques by reassigning the TF representations by the rule of centroids. In Fig. 8(a) and Fig. 8(b), we show the results of the RM-GT and the RM-SPWVD, representatives for the linear and quadratic type methods, respectively, can address the broadening issue in the original methods and provide distinct inner arches as well as other substructures. Note that in Fig. 8(b), the arch in the duration of 2 − 2.5T and around the 45th harmonic is an example of interference that is difficult to remove. It is because that the RM techniques re-shuffle the original TF representation based on its local center of mass, and hence preserves the intrinsic artifact in the adopted transform.
By a reallocation rule that contains phase information of a linear type TF transform, the SST can enhance the TF representation by linear type TF transforms in a manner similar to that by the RM technique, as shown in Fig. 8(c) and (d). Both the SST-GT and the SST-MWT can depict similar structure with appropriate the STFT and CWT parameters. Note that as mentioned previously, the SST-GT inherited the fixed resolution feature in the STFT, and the SST-MWT maintains the adaptive resolution in the CWT. By comparing Fig. 8(c) with Fig. 8(c), we see that the resolution in the RM-GT is more sharpen than that of the SST-GT. One of the reasons is that while both temporal and frequency reassignments are considered in the conventional reassignment method, temporal reassignment is not taken into account in the SST (both SST-GT and SST-MWT). In addition, the reallocation rule Eq. (20) is only the first order approximation for the IF, resulting in diffusive pattern particularly for fast-varying chirp signal. For the simulation in the multiphoton ionization regime, such diffusive pattern is negligible as the IF changes slowly. By using Eq. (21), a second-order reallocation rule for the SST-GT, the diffusive pattern can be further concentrated, as is shown in Fig. 8(e).
To shed some light on this TF representation for the tunneling ionization mechanism, we follow the standard semiclassical approach suggested independently by Corkum [56] and Kulander et al. [57]. Here the electric-field force corresponding to the applied laser field in a.u. is F z = −E 0 E(t)ẑ , whereẑ is the unit vector in the z−direction. We assume the initial position and initial velocity is both zero [29].
The trajectory of electrons released between 1T and 1.5T by using the extended semiclassical approach is denoted by red circles in Fig. 8(f). This result is superimposed for the sake of comparison with the TF representation of the SST-GT. The trajectory suggests that the largest arch in Fig. 8 is associated with the first return time of the electron wave packet released in between 1T and 1.5T of the laser field, and the smaller arch between 2T and 2.5T , is the second return.
Although the SST depicts the first and multiple returns clearly, the sequence of these returns cannot be known. To retrieve the first and second return quantitatively, Risoud et al. propose the application of the Hilbert transform on the dipole moment with the harmonics lower than the ionization potential filtered [18]. While there are more than one component simultaneously exist at each time, we consider the EMD scheme to decompose the filtered dipole moment. (In this case, harmonics below the 10th harmonic are filtered as they correspond to the bound part of the atomic spectrum.) The stopping criterion in the EMD scheme is that the residue becomes monotonic. The HS for all modes of IMF are displayed in Fig. 9(a). Note that a median filter with a window length of 51 points is applied on all modes as the implementation of Eq. (28) can introduce numerical errors. Despite there is no theoretical base for the IMFs by the EMD and the result is not consistent in the multiphoton ionization case, in this case we see that the first two modes can be associated with the first and second returns. Fig. 9(b) presents the HS for the first two IMFs and Fig. 9(c) compares the IMFs with the representation of the SST-GT. In Fig. 9(b), the IMFs are visually emphasized by additional Gaussian functions in both temporal and frequency domain. Note that the IF revealed by the Hilbert transform is meaningful only if the IMF is consist of a single component, which is not guaranteed in the EMD. For example, for time greater than 2.5 cycles, the HS of the first IMF is distorted because of more wave packet returning from previous emissions. While the Hilbert transform is not theoretically suitable for this kind of oscillatory signal with time-varying AM and IF, we illustrate the IF of these decomposed modes by the SST-GT. We employ the SST-GT on the first two IMFs, as presented in Fig. 9(d) and (e). Results in Fig. 9(d) and (e) suggest that expressing the IF of the IMF by the SST is more stable than using the Hilbert transform. Note that no median filter is necessary for each IMFs in the SST-GT analysis, which reduces the possible artifacts.
A summary for comparison of TF methods in the tunneling ionization regime is provided in Table 3.

Conclusions
This paper provides a benchmark of analyzing time-dependent quantum systems by applying several contemporary TF methods. Within the same physical model and computational scheme, different features of TF representations provided by different TF methods may lead to conflicting interpretation in physics. In the multiphoton ionization regime, linear TF methods with an appropriate window function, can vaguely depict discrete odd harmonics in the HHG process, yet the detail features such as the AC Stark effect cannot be revealed. In addition, in the tunneling ionization regime, the first-return trajectories are roughly illustrated, but the multiple returns are indistinct because of the broadening induced by the window. While present little broadening, TF representations by the WVD suffer from artificial interference patterns and lead to incorrect interpretation of physical mechanisms, such as even harmonics in the multiphoton ionization regime and false trajectories in the tunneling ionization regime. To eliminate the artifacts mentioned above, several modern TF analysis methods are introduced in quantum dynamics for the first time. However, Cohen class distributions such as the SPWVD can only remove portion of interference patterns. Reassignment techniques and the SST can depict accurate dynamic process in the two regimes based on the principle of local centroid and instantaneous frequency, respectively. However, there is no inverse transform in the reassignment techniques. On the other hand, the SST preserves signal causality and allows mode reconstruction. We demonstrate that after removing the broadening, the AC Stark effect in the multiphoton ionization regime and the multiple returns in the tunneling ionization regime are revealed. Finally, we relate the IMFs from the EMD to the multiscattering of the electron energy distribution in the tunneling ionization regime.
We believe that in addition to the atomic hydrogen system, the contemporary TF methods have potential for exploring other complicated quantum systems. In case of analyzing dynamics in an unknown quantum system, this paper suggests the following procedure. First, features obtained by the WVD can be used to determine the parameter for the linear type methods. Then the SST can be applied to obtain TF representations for further interpretation. We hope that this research can serve as a cornerstone in applications of TF analysis for fields including but not limited to, attosecond physics, nuclear magnetic resonance, ultrafast dynamics in atoms and molecules.

Acknowledgments
The authors would like to thank Dr. Elise Y. Li and Yu-lin Sheu for their help in servers.
Appendix A. Adaptive harmonic model Fix constants 0 ≤ 1, 0 < d < 1, c 2 > c 1 > and c 2 > c 3 > . Consider the functional set Q c1,c2,c3 , which consists of functions x(t) = A(t) cos(2πφ(t)) ∈ C 1 (R) ∩ L ∞ (R) where the following conditions hold: We call x = A(t) cos(2πφ(t)) ∈ Q c1,c2,c3 an intrinsic mode type (IMT) function, φ(t) the phase function, the derivative of the phase function φ k (t) is called the instantaneous frequency (IF) and A(t) is called the amplitude modulation. When |φ (t)| φ (t) for all t ∈ R, we call the signal x oscillatory with slowly varying IF; otherwise we call it oscillatory with fast varying IF. Locally an IMT function with slowly varying IF behaves like a harmonic function and an IMT function with fast varying IF behaves like a linear chip function. IMT function serves as a mathematical formula for the IMF considered in the EMD algorithm, i.e., the IMT functions contain the properties defined in the IMF, but the reverse inclusion does not hold.

Appendix B. The Floquet Method
The total Hamiltonian of atomic hydrogen can be expressed as H(t) = H 0 + H 1 (t), where H 0 is the unperturbed hydrogen Hamiltonian and H 1 (t) = −zE(t) describes the laser-atom interaction. In the timedependent generalized pseudospectral simulation, the time period T is fixed but the field amplitude varies. To calculate the shifting of orbital energies, we separate each cycle and take its field amplitude, and then assume the hydrogen in a laser field with the particular constant amplitude. For example, for the third cycle in Fig. 1(a), the maximal amplitude is 4.86 × 10 −4 , and we assume hydrogen in the laser field with the constant amplitude.
When the Hamiltonian is periodic with a time period T , i.e., H(t) = H(t + T ), we can apply the Floquet theory [53] and derive a time-independent infinite-dimensional eigenvalue matrix equation [54] as follows: Here nlmk| H F |n l m k denotes the matrix elements of the time-averaged Floquet Hamiltonian over a period T , where the outer bra-ket notation refers to the inner product over t. We use the orbitals of an unperturbed hydrogen as a basis set since we only consider a weak field regime. As a result, n, l, and m represent the principal quantum number, angular momentum quantum number, and magnetic quantum number of hydrogen, respectively. k is an integer.