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

EXPLORING LASER-DRIVEN QUANTUM PHENOMENA FR OM A TIME-FREQUENCY ANAL YSIS PERSPECTIVE: A COMPREHENSIVE STUD Y Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU Abstract. Time-frequency (TF) analysis is a pow erful tool for exploring ultrafast dynamics in atoms and molecules. While some TF methods ha v e demonstrated their usefulness and p otential in several of quantum systems, a systematic comparison among these methods is still lac king. T o this end, we compare a series of classical and con temp orary TF metho ds by taking hydrogen atom in a strong laser field as a b enchmark. In addition, sev eral TF methods such as Cohen class distribution other than the Wigner-Ville distribution, reassignment methods, and the empirical mo de decomp osition metho d are first introduced to exploration of ultrafast dynamics. Among these TF metho ds, the synchrosqueezing transform successfully illustrates the physical mechanisms in the multiphoton ionization regime and in the tunneling ionization regime. F urther- more, an empiric al pro cedure to analyze an unknown complicated quantum system is provided, indicating the versatilit y of TF analysis as a new viable ven ue for exploring quantum dynamics. 1. Introduction Time-dep enden t quan tum mechanics is a fundamen tal topic in physics, chemistry , and engineering. T ra- ditionally , dynamics of a quan tum system such as lifetime and energy difference b et w een states can b e rev ealed by spectral lines in the frequency domain by using the F ourier analysis. Ho wev er, the F ourier trans- form presen ts limited c hronological information of a dynamical system. T o explore chronological information in a quantum dynamical system, time-frequency (TF) metho ds are necessary . In recen t experiment adv ances, TF methods such as wa velet transform combined with ultrafast spec- troscop y tec hniques w ere used to probe dynamics in molecules, solids and liquids [1, 2, 3, 4, 5, 6, 7], particular in ligh t harvesting complexes[8, 9, 10]. In attosecond physics, the short-time F ourier transform (STFT), the con tinuous wa velet transform (CWT), Wigner-Ville distribution (WVD) (which b elongs to the Cohen class distribution), were employ ed to uncov er dynamical mechanisms including high-order harmonic generation (HHG) [11, 12, 13, 14, 15, 16, 17]. Lately , Hilb ert transform and the synchrosqueezing transform (SST) were applied to inv estigate multiple scattering [18] and the dynamic origin of near- and b elow-threshold harmonic generation, resp ectiv ely [19, 20]. While these studies sho w 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 quan tum system with unknown complicated dynamics. First, different types of TF metho ds and the c hoice of window functions ma y pro vide extremely differen t TF represen tation, leading to conflicting ph ysical in terpretation. Second, it is difficult to select a particular TF metho d based on previous studies, in which the ph ysical mo dels are different, e.g. , 1D [12, 18], 2D [12] and 3D atoms [12, 15, 16, 19, 20], and solved by differen t n umerical sc hemes [12, 15, 20, 18]. Clearly , TF representations derived from differen t TF metho ds based on differen t physical mo dels cannot provide an impartial comparison. Third, in the past decades, sev eral mo dern TF metho ds, e.g. , the Cohen class distribution, empirical mo de decomp osition with Hilb ert sp ectrum (EMD-HS), reassignment metho ds (RM), and SST hav e b een prop osed and succ essfully applied to classical macroscopic dynamical system s including molecular dynamics [21], cardiopulmonary coupling phe- nomena [22], c hronotaxic systems [23], lamb wa v e propagation [24] and seismic data [25]. Ho wev er, Cohen class distribution other than the WVD, EMD-HS, RM, and different forms of SST ha v e not b een discussed in quantum dynamical systems. As a consequence, a comprehensiv e TF study in the same b enchmark is essen tial. T o this end, w e take 3D h ydrogen atom as a benchmark b ecause h ydrogen is one of the most represen tativ e systems in quantum mechanics. All sim ulations are based on the same numerical metho d (time-dep endent generalized pseudosp ectral metho d), and different contemporary TF analysis, including the STFT, CWT, 1 2 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU Cohen class distributions, RM, SSTs and EMD-HS, are applied to study the sim ulated signals. The system- atic comparison of the TF metho ds enables to interpret physical pro cesses in a quantum system. This article is or ganized as follo ws. In Section 2, w e summarize the ph ysical mo del of the h ydrogen system in a strong laser field and the n umerical simulation. In Section 3, we pro vide an o v erview of the TF methods, as well as a discussion of their pros and cons. In order to show that quantum dynamical phenomena can b e depicted by the TF metho ds, we p erform simulations with tw o different sets of laser parameters, including h ydrogen in the regime of m ultiphoton ionization and tunneling ionization. The results of TF represen tations b y different TF methods and a discussion of their limitations and p otential applications are giv en in Section 4. In Section 5 we conclude the pap er and provide an empirical pro cedure to analyze an unkno wn complicated quan tum system. 2. Theoretical modeling W e simulate hydrogen dynamics in a linear p olarized laser field in the framework of the electric dip ole appro ximation and non-relativistic quan tum mec hanics. Note that for the wa velength range of approximately 500 to 1200 nm, the electric dip ole approximation is v alid [26] only when the laser intensit y I 0 is smaller than 10 15 ∼ 10 16 W / cm 2 . The time-dep enden t Schr¨ odinger equation for atomic hydrogen interacting with a linearly p olarized field along the z -axis in atomic units can b e expressed as i ∂ ψ ( r , t ) ∂ t =  − 1 2 ∇ 2 − 1 r − z E ( t )  ψ ( r , t ) , (1) where ψ ( r , t ) is the wa ve function at p osition r = ( x, y , z ) and at time t , z = r cos( θ ), and E ( t ) is the external laser field. T o numerically solve this equation, w e adopt a time-dep endent generalized pseudosp ectral method [27, 28] whic h consists of tw o essen tial steps: (1) The spatial co ordinates are optimally discretized in a nonuniform fashion by means of the generalized pseudosp ectral technique: the grid is denser near the origin and sparser a wa y from the origin; (2) A second order split op erator technique in the energy representation, which allows the explicit limitation of undesirable fast-oscillating high energy components, is used to obtain an efficien t and accurate time propagation of the wa ve function. According to previous studies [29, 30], dynamical phenomena such as the HHG is asso ciated with the electric dipoles, i.e., the laser-driv en electron oscillating around the stationary n ucleus, which could b e expressed as, resp ectively , the time-dep enden t induced dip ole in the length form, denoted as d L ( t ), and in the acceleration form, denoted as d A ( t ) [28]: (2) d L ( t ) = Z ψ ∗ ( r , t ) z ψ ( r , t )d r . (3) d A ( t ) = Z ψ ∗ ( r , t ) d 2 z d t 2 ψ ( r , t )d r . 3. Time-frequency Anal ysis Supp ose the signal x ( t ) is comp osed of finite oscillatory components, that is, x ( t ) = P L l =1 a l ( t ) cos(2 π φ l ( t )), and each comp onen t has a time-v arying amplitude mo dulation (AM), a l ( t ) > 0, and time-v arying instan ta- neous frequency (IF) φ 0 l ( t ), where 0 stands for the first deriv ativ e, then the TF representation reflects the IF, whic h is a generalization of the notation frequency , and the lo calized phase information could b e extracted. W e call a l ( t ) cos(2 π φ l ( t )) an in trinsic mode type (IMT) function. More details ab out this kind of function could b e found in App endix A. T o disclose the time-v arying nature of this kind of signal, in particular the IF and AM, the TF analysis is a p ow erful approac h. In the follo wing subsections w e provide an o v erview on classical and contemporary TF metho ds. 3 3.1. Linear TF metho ds. The global analysis nature of the F ourier transform is resp onsible for its lim- itation in extracting time-v arying dynamics inside an oscillatory signal. An intuitiv e w a y to resolve this limitation is analyzing the signal lo cally; that is, we could crop a segment of finite length and apply the F ourier transform, and exp ect to observe how a signal oscillates lo cally . This idea leads to the STFT. W e in tro duce a window function g to crop the signal at different time and p erform the F ourier analysis. This analysis results in a time-frequency (TF) representation of the signal [31]: STFT g x ( t, ω ) = Z ∞ −∞ x ( ζ ) g ∗ ( ζ − t ) e − iω ( ζ − t ) d ζ . (4) W e call STFT g x ( t, ω ) a TF representation of the given function x and | STFT g x ( t, ω ) | 2 the sp ectrogram of x . Here the STFT by Eq. (4) is different from the conv entional definition by an additional mo dulation factor e − iω t . Nev ertheless, temporal and frequency resolutions cannot be ac hieved simultaneously b y the STFT, accord- ing to the Heisenberg uncertaint y principal [31]. In other words, a wide windo w provides a go o d frequency estimation in the STFT at the cost of p o or temp oral resolution, while a narrow window has the opp osite trade-off. In this research w e follo w the tradition and c ho ose 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 , whic h leads to the Gab or transform (GT): GT g x ( t, ω ) = Z ∞ −∞ x ( ζ ) e − ( ζ − t ) 2 2 σ 2 e − iω ( ζ − t ) d ζ . (5) Here the GT we apply here differs from the usual one by an additional mo dulation factor e iω t . Based on the same lo cal analysis idea as that of the STFT, we could analyze the momentary b eha vior of the signal x ( t ) by the CWT: CWT g x ( t, a ) = Z ∞ −∞ x ( ζ ) 1 √ a g ∗  ζ − t a  d ζ , (6) where g is the c hosen mother w a v elet 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 represen tation analyzed by the CWT has a go o d time resolution and p o or frequency resolution at the high frequency region and a go o d frequency resolution and p o or time resolution at the low frequency region. One of the commonly applied mother w a velet is the Morlet w av elet, and the CWT based on the Morlet wa v elet is called the Morlet w av elet transform (MWT): MWT g x ( t, ω ) = Z ∞ −∞ x ( ζ ) √ ω g ∗ ( ω ( ζ − t ))d ζ , (7) where ω > 0, g ( ζ ) = 1 √ τ e − ζ 2 2 τ 2 e iζ is the Morlet wa v elet and τ > 0. The standard deviation for the dilated Morlet wa velet is σ = τ /ω . In Eq. (7), w e follow the con ven tion and use ω , whic h is the in v erse of the scale parameter a in Eq. (6). Note the fundamental difference b etw een Eq. (7) and Eq. (5) – in the GT, the Gaussian function is of a fixed width but in the MWT the width v aries. In the MWT the window width σ v aries with ω such that the quality factor τ (the inv erse of the relative bandwidth) is in v ariant on the TF plane [31]. In other w ords, the window width σ b ecomes smaller as ω increases, and vice versa. Since the frequency resolution dep ends on the scale, we say that the TF representation by CWT has an adaptive fr e quency r esolution . Both the STFT and the CWT b elong to the line ar typ e TF analysis , in which the signal is characterized b y their inner pro ducts with a preassigned family of templates with free parameters. Clearly , these TF represen tations dep end on the chosen window, whic h migh t cause artificial patterns on the analysis result. F urther, while in general the underlying structure of the signal under analysis is not kno wn, there is no systematic wa y to design the windo w which faithfully reflects the structure. The ab o ve facts render the linear type TF analysis non-adaptive to the signal under analysis. 4 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU 3.2. Quadratic TF metho ds. T o resolve the non-adaptivity issue of the linear type TF analysis as w ell as finding the higher order structure, the Wigner-Ville distribution (WVD) [31] was prop osed. The WVD is based on the concept of the auto correlation function and is defined as WVD x ( t, ω ) = Z ∞ −∞ x ( t + ζ / 2) x ∗ ( t − ζ / 2) e − iω ζ d ζ . (8) Note that by a slight change of v ariable, we could view the WVD as a v ariant STFT which is free of the windo w choice issue. In this sense, WVD is adaptive to the signal under analysis. Also note that by a direct deriv ation, the sp ectrogram could b e understo o d as follows [31]: | STFT g x ( t, ω ) | 2 = Z Z ∞ −∞ WVD x ( ζ , η )WVD g ( ζ − t, η − ω )d ζ d η . (9) The WVD has s ev eral go o d mathematical prop erties, e.g., the signal energy is preserv ed in the WVD, the WVD is a real-v alued function on the TF-plane and it provides a precise information ab out the chirp signal, lik e x ( t ) = e i 2 π ( β 0 + β 1 t + β 2 t 2 ) , where β 0 ∈ R , β 1 > 0 and β 2 ∈ R . How ev er, choosing the signal itself as the windo w function causes sp ecific artifacts depending on the signal type. F or example, when there are more than one oscillatory component in the signal, the in terference patterns in the TF represen tation is inevitable, whic h migh t lead to mis-in terpretation of the signal. By construction, the WVD is quadratic in the signal which could b e view ed as an energy distribution. A direct generalization of the WVD based on imp osing some constraints on the c ovarianc e structur e [31] of the signal leads to the Cohen’s class. In other words, the Cohen’s class , which comprises all bilinear TF represen tations that are cov ariant to shifts in b oth time and frequency [31] and has the following form: C x ( t, ω ) = Z Z Z x ( s + ζ / 2) x ∗ ( s − ζ / 2) e − iω ζ e − iξ ( s − t ) f ( ξ , ζ )d ξ d s d ζ , (10) where f ( ξ , ζ ) is an arbitrary parameter function. Different parameter functions lead to different TF repre- sen tations. 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 oscil- latory comp onent. T o b e more sp ecific, w e can employ a separable parameter function f ( ξ , ζ ) = G ( ξ ) h ( ζ ) ≡ W h ( ξ , ζ ) in Eq. (10), where suitably c hosen G ( ξ ) and h ( ζ ) permit a contin uous and indep endent control of the in terferences in time and frequency , resp ectiv ely . The corresponding represen tation is called the smo othe d pseudo-Wigner-Vil le distribution (SPWVD) [31] and could b e rewritten as SPWVD g ,h x ( t, ω ) = Z ∞ −∞ g ( t − s ) Z ∞ −∞ h ( ζ ) x ( s + ζ / 2) x ∗ ( s − ζ / 2) e − iω ζ d ζ d s (11) = Z Z ∞ −∞ WVD x ( s, ξ )W h ( s − t, ξ − ω )d s d ξ , (12) where g ( ζ ) is the inv erse F ourier transform of G ( ξ ). Particularly , G ( ξ ) and h ( ζ ) are even functions with g (0) = 1 and H (0) = 1, where H ( ξ ) is the F ourier transform of h ( ζ ). W e mention that we could consider a v ariation 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 co v ariant and the TF representation also has the adaptive resolution feature. T o simplify the discussion, w e do not study the affine class transforms, but we mention that their p erformance is similar to those in the Cohen’s Class. The ab o v e-mentioned transforms are ov erall called the quadr atic TF analysis . A higher-order generalization is also p ossible, and we refer the reader of interest to, for example, [33, 34]. 3.3. Reassignmen t Metho ds. As features of a TF represen tation computed b y a linear t yp e of TF method are smeared by the introduced window functions and those by the quadratic t yp e ones are obscured b y in terferences, the RM was prop osed to sharp en the resolution of the TF representation. Generally sp eaking, the co efficients of the TF represen tation at ( t, ω ) is reallo cated to a different p oint ( ˆ t, ˆ ω ) according to a predefined reallo cation rule [35, 36, 37]. A common choice of the reallo cation rule is to assign v alues of a TF represen tation to the lo cal centroids. Note that this is different from the a veraging idea b ehind the STFT 5 (Eq. (9)). In this study we apply the RM to the STFT and SPWVD, and called the metho ds RM-STFT and RM-SPWVD, resp ectively . The reallo cation rule for the RM-STFT is derived by estimating the lo cal centers of gravit y , denoted as: ˆ t x ( t, ω ) =    <  STFT tg x ( t, η )(STFT g x ( t, η )) ∗  | STFT g x ( t, η ) | 2 , when STFT g x ( t, η ) 6 = 0 ∞ , when STFT g x ( t, η ) = 0 . (13) ˆ ω x ( t, ω ) =      −= n STFT d g x ( t, η )(STFT g x ( t, η )) ∗ o | STFT g x ( t, η ) | 2 , when STFT g x ( t, η ) 6 = 0 ∞ , when STFT g x ( t, η ) = 0 . (14) The notations tg and d g stand for the windo w function tg ( t ) and the first deriv ativ e of g ( t ) in the STFT, resp ectiv ely . The RM-STFT is thus defined as RM-STFT g x ( ˆ t, ˆ ω ) = Z Z ∞ −∞ | STFT g x ( t, ω ) | 2 δ ( ˆ t x ( t, ω ) − ˆ t, ˆ ω x ( t, ω ) − ˆ ω )d t d ω . (15) Essen tially , this formula reassign the energetic conten ts of the sp ectrogram to the new lo cation ( ˆ t x , ˆ ω x ). As a consequence, the RM leads to a substantially improv ed resolution in the TF representation. Note that the reassignmen t 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 un wrapp ed phase of the STFT g x ( t, ω ). Ho wev er these physical meaningful expressions are numerically inefficien t. Other forms of reassignment op erators can b e found in [31]. Next we consider SPWVD. Although the SPWVD can smo oth out the interferences in the WVD, the smo othing functions introduce artificial broadening. Applying the reassignment technique on the TF repre- sen tation of the SPWVD can reduce such artifacts. Similar to Eq. (15), the reassigned representation of the SPWVD is RM-SPWVD x ( ˆ t, ˆ ω ) = Z Z ∞ −∞ SPWVD g x ( t, ω ) δ ( ˆ t x ( t, ω ) − ˆ t, ˆ ω x ( t, ω ) − ˆ ω )d t d ω . (16) where the reassignment rule ( ˆ t x , ˆ ω x ) is determined by the concept of exp ectation: ˆ t x ( t, ω ) = Z Z ∞ −∞ s WVD x ( s, ξ )W h ( s − t, ξ − ω )d s d ξ (17) ˆ ω x ( t, ω ) = Z Z ∞ −∞ ξ WVD x ( s, ξ )W h ( s − t, ξ − ω )d s d ξ . (18) Despite RMs are intuitiv e techniques to sharpen the linear type and quadratic t ype TF representations, in verse routines as well as mo de reconstruction are not av ailable. 3.4. Sync hrosqueezing transform. In this section we describ e the SST, which is a sp ecial RM aiming to address the intrinsic blurring issue in the linear TF metho ds. T o b e more precise, the SST manifests IF c haracteristic according to the reallo cation rule that consists of solely the frequency information, rather than the cen troid of the TF represen tation [38, 39]. In addition to sharpening the TF represen tation, the causality prop ert y of the signal is preserved in the SST, which allows the decomp osition of osc illatory comp onen ts when the signal is comp osed of several oscillatory comp onen ts. W e refer the reader who has interest in the theoretical analysis and other details, like reconstructing eac h oscillatory comp onent and robust to noise, of SST to [38, 39, 40]. In this section we describ e the synchrosqueezed STFT (SST-STFT) and the sync hrosqueezed CWT (SST-CWT). 6 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU 3.4.1. SST-STFT. Consider a STFT g x ( t, ω ) of x ( t ) with a window function g ( t ) suc h that supp( G ( ω )) ⊂  − d 2 , d 2  , where G ( ω ) is the F ourier transform of g and d is the smallest gap of IFs of any t w o consecutive oscillatory comp onents. The SST-STFT is given as SST-STFT g x ( t, ω ) = Z ∞ −∞ STFT g x ( t, η ) 1 α √ π e − | ω − ˆ ω x ( t,η ) | 2 α d η , (19) where ˆ ω x ( t, η ) is the reallo cation rule and α > 0 is a controllable smo othing parameter for the resolution, whic h in practice is chosen to b e small. The reallo cation rule given by the following equation utilizes the phase information hidden inside the smeared TF representation: ˆ ω x ( t, η ) =    − i∂ t STFT g x ( t, η ) STFT g x ( t, η ) when STFT g x ( t, η ) 6 = 0 ∞ when STFT g x ( t, η ) = 0 . (20) Note that only the frequency reassignment operator is considered in the SST-STFT, so that the causality prop ert y of the signal can b e preserv ed. A sligh t modification of ˆ ω ( t, η ) based on the second-order information of the IF [41] migh t further impro v e the TF representation. Indeed, when the interested oscillatory component could be approximated b y a chirp function [41, 36], although the TF represen tation is sharpened by the SST- STFT, a mild spreading is inevitable. W e call this mild spreading the diffusive p attern . T o cop e with the diffusiv e pattern, w e could consider the following reassignment rule: ˇ ω x ( t, η ) =  ˆ ω x ( t, η ) + c ( t, η )( t − ˆ t x ( t, η )) when ∂ ω ˆ t x ( t, η ) 6 = 0 ˆ ω x ( t, η ) otherwise , (21) where ˆ t x ( t, η ) is defined as ˆ t x ( t, η ) = t + i ∂ ω STFT g x ( t, η ) STFT g x ( t, η ) and c ( t, η ) = ∂ t ˆ ω x ( t, η ) ∂ ω ˆ t x ( t, η ) . (22) The reconstruction formulas [39] for each IMT function x k ( t ) from SST g , STFT x ( t, ω ) is x k ( t ) = < ( C − 1 g Z ω k ( t )+  2 ω k ( t ) −  2 SST-STFT g x ( t, ω )d ω ) (23) 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. F or a CWT with a mother wa v elet g such that supp( G ( ω )) ⊂ [1 − ∆ , 1 + ∆], with ∆ < d 1+ d , the SST-CWT is given by SST-CWT g x ( t, ω ) = Z ∞ −∞ η − 3 / 2 CWT g x ( t, η ) 1 α √ π e − | ω − ˆ ω x ( t,η ) | 2 α d η , (24) where α > 0 is a smo othing parameter and ˆ ω x ( t, η ) is the reallo cation rule defined by ˆ ω x ( t, η ) =    − i∂ t CWT g x ( t, η ) CWT g x ( t, η ) when CWT g x ( t, η ) 6 = 0 ∞ when CWT g x ( t, η ) = 0 . (25) When describing an oscillatory comp onen t with a fast v arying IF, the TF representation in SST-CWT ma y sho w a slight spreading. While the second-order information of the IF has been discussed in the SST-STFT[41], it is not yet b een studied. The reconstruction formula [39, 40] for each IMT function x k ( t ) under SST-CWT is x k ( t ) = < ( R − 1 g Z ω k ( t )+  2 ω k ( t ) −  2 a − 3 / 2 SST-CWT g x ( t, a )d a ) (26) where   1, and R g := R ∞ −∞ G ( η ) η d η . The results of the SST are “adaptive” to the signal in the sense that the error in the estimation dep ends only on the first three moments of the mother wa v elet instead of the profile of the mother wa v elet [38, 39, 40]. In other words, the influence of the c hosen window on the asso ciated TF representation is minimized 7 compared with the linear TF metho ds. In addition to the ab ov e prop erties, the SST is robust to several differen t kinds of noise, which might b e slightly non-stationary [40]. An imp ortan t prop erty shared by the SST-STFT, SST-CWT and RM-STFT is that by taking a short window, the fast v arying IF could b e w ell captured. A t the first glance, the results of the reassignmen t technique and SST seem to break the well-kno wn Heisen b erg uncertain principle, which sa ys that the temp oral and frequency resolution cannot b e achiev ed sim ultaneously . Ho w ev er, in the reassignment technique and SST, we hav e shown that we could obtain a TF represen tation with an almost p erfect time and frequency resolution. The main reason is that w e focus on the oscillatory signals, in particular the adaptive harmonic mo del but not the whole L 2 space. W e mention that the b ehavior of the RM and SST on the general L 2 function is an op en question. 3.5. EMD. One p opular wa y to define the IF of a given signal x ( t ) is via finding the analytic representation of x ( t ) with the aid of the Hilb ert transform: ˜ x ( t ) = x ( t ) + i H ( x ( t )) = a ( t ) e i Φ( t ) , (27) where H ( ˜ x ( t )) denotes the Hilb ert transform of x ( t ). In this equation, a ( t ) and Φ( t ) are the mo dulus and phase of ˜ x ( t ), resp ectiv ely [42, 43]. The IF of x ( t ) is thus defined as the rate of the v arying phase: ω ( t ) = 1 2 π dΦ(t) d t . (28) Ho wev er, the IF pro vided by Eq. (28) is not alw ays meaningful. T o deal with this problem, the empirical mode decomp osition algorithm (EMD) [44] was prop osed to decomp ose x ( t ) into a series of functions called the in trinsic mo de function (IMF), x = P K k =1 x k ( t ), on which the Hilb ert transform can be applied subsequently , via an algorithm called the sifting pro cess. IMFs must satisfy the following tw o c onditions: (i) In the whole data set of a signal, the num ber of extrema and the num ber of zero crossings must either b e equal or differ at most by one; (narro w band) (ii) A t any p oint, the mean v alue of the env elop e defined b y the lo cal maxima and the en velope defined by the lo cal minima is zero (adoption of lo cal prop erties). Note that in general an IMF is different from an IMT function. Details of the EMD algorithm could b e found in, for example, [44]. F or each IMF x k ( t ), the corresp onding IF and AM amplitudes can b e estimated by Eq. (27) via the Hilb ert transform and Eq. (28). By assigning the IF and amplitudes of all IMFs on a TF plane, we obtain the Hilb ert sp ectrum (HS) HS x (t , ω ) (HS), which is a TF representation for x ( t ) determined by HS x ( t, ω ) = X k a k ( t ) δ ( ω − ω k ( t )) , (29) where δ ( ω ) denotes the Dirac delta measure. Although the EMD along with HS has b een applied to several fields, due to a num b er of heuristic and ad ho c elemen ts in the EMD algorithm, it is difficult to analyze its accuracy and limitation. In addition, despite the solid mathematical supp ort, the Hilb ert transform might b e limited when applied to analyze the momentary dynamics of an oscillatory signal. First, due to the slow deca y nature of the kernel in the Hilb ert transform, keeping the causality of the signal structure might b e difficult. Second, when the signal has time-v arying amplitude and frequency , in general there is no guarantee to get the correct analytic signal from a given real oscillatory signal by the Hilb ert transform in Eq. (27). W e mention that the part of using the Hilb ert transform to obtain the AM and IF of each IMF could b e replaced b y using the SST-STFT or SST-CWT. In this case, the metho d is regarded as the EMD-SST-STFT and EMD-SST-CWT. A comparison of different TF metho ds discussed in this study is summarized in T able 1. 4. Resul ts and Discussions In this section we emplo y the aforementioned TF metho ds on the time-dep endent dipole computed by the time-dep enden t generalized pseudosp ectral metho d at the ab initio lev el. In a strong laser field, a v ariety of dynamical pro cesses can o ccur in atoms and molecules, such as p onderomotive effect, the AC Stark effect, HHG [46], multiphoton ionization[47, 48, 49], ab o ve threshold ionization, tunneling ionization, rescattering of electron w av epack et [50, 51], etc . Up to date, the STFT [11, 12, 13], CWT [14, 16, 15], WVD [17], and Hilb ert transform [18] hav e b een engaged in inv estigation of estimation of emission time and multiple 8 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU T able 1 Summary of TF metho ds in This Study T yp e Multiresolution Choice of In verse Artifacts in TF Representation P arameters T ransform STFT No Y es Y es broadening by the window function CWT Y es Y es Y es broadening b y the windo w function WVD No No Y es in terference patterns SPWVD No Y es Y es broadening by the filter function RM-STFT No Y es No causalit y is not preserved RM-SPWVD No Y es No causalit y is not preserved SST-STFT No Y es Y es Diffusiv e pattern in the fast-v arying IF 2nd order SST-STFT No Y es Y es SST-CWT Y es Y es Y es Diffusiv e pattern in the fast-v arying IF EMD-HS Y es? † Y es ‡ Y es/No  mo de mixing, etc. ∗ EMD-SST-STFT Y es? † Y es ‡ Y es/No  mo de mixing † : The multiresolution-lik e b ehavior of the sifting pro cess was studied in [45]. ‡ : The stopping criteria of the sifting pro cess dep ends on tuning several parameters.  : The inv ersion from the HS is not p ossible for the definition (29). It is p ossible if we define the HS as P k a k ( t ) e i Φ k ( t ) δ ( ω − ω k ( t )). ∗ : It is sensitive to noise; a shortly existing oscillatory comp onent destroys the whole analysis. scattering. How ev er, to the b est of our understanding there is no comprehensiv e study comparing these TF metho ds in a quantum dynamical system. By taking the w ell-studied atomic hydrogen as a benchmark, here we fo cus on the laser-driv en h ydrogen in the m ultiphoton ionization regime and in the tunneling ionization regime. When the laser photon energy ~ ω is m uch smaller than the ionization p oten tial I p , the m ultiphoton ionization and the tunneling ionization can be describ ed by a dimensionless Keldysh parameter γ K = p I p / 2 U P , where U P = E 2 0 / 4 ω 2 0 is the p onderomotive p oten tial [29, 48]. Similar ph ysical dynamics can b e achiev ed for any atom-field interaction given a fixed γ K [52]. Generally sp eaking, multiphoton ionization of atoms b ecome dominant when γ K  1, while tunneling ionization is predominant when γ K  1 [29]. In the following we p erform simulations with laser parameters in these tw o regimes. 4.1. Multiphoton Ionization Regime. In the first sim ulation, the laser field parameters are arranged suc h that the Keldysh parameter is γ K = 3 . 07, indicating that multiphoton ionization is the ma jor mechanism. The laser wa v elength is 880 nm, which corresp onds to ω 0 ≈ 0 . 05178273 in the atomic unit (a.u.), and the laser intensit y of I 0 = 10 13 W/cm 2 , which corresp onds to a field amplitude E 0 ≈ 0 . 0169 (a.u.). Note that for this γ K v alue, b oth tunneling and multiphoton ionization can o ccur, 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 dip ole in the length form d L ( t ) is presented in Fig. 1(b). The p ow er sp ectrum [28], computed by the squared F ourier sp ectrum, of the laser field and d L ( t ) are shown in Fig. 1(c)and 1(d), resp ectiv ely . Both the length and the acceleration forms of dip ole moment presen t the same detail structures in their p ow er sp ectra [28]. Here w e present the analysis using d L ( t ) for the multiphoton ionization case. Although the profiles of d L ( t ) and the applied laser field lo ok similar in the time domain, they are differen t in the frequency domain. F or example, while the p o w er sp ectrum of the laser field has only one p eak lo cated at ω 0 , that of d L ( t ) reveals o dd harmonics due to the parit y symmetry [29]. Ho w ever, the meaning of the substructures within the o dd harmonics and their corresp onding dynamics are unclear. T o unv eil the dynamics of d L ( t ), first we apply the con ven tional linear and quadratic TF metho ds, as sho wn 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 display ed in Fig. 2(a) and Fig. 2(b), resp ectiv ely . Both the GT and MWT depict separate broad lines regarding the o dd harmonics in the HHG pro cess. While frequencies of the extracted comp onen ts are consisten t with the information shown in the F ourier sp ectrum in Fig. 1(d), the GT and MWT further capture the momentary b eha vior of each comp onen t. F or a fixed resolution in the GT, low frequency comp onents could b e lost given an insufficient window width. In the TF representation 9 0 10 20 30 40 50 60 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 Time (Cycle) Laser Field (atomic unit) 0 10 20 30 40 50 60 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 Time (Cycle) d L (t) (atomic unit) (a) (b) 0 5 10 15 20 25 30 −25 −20 −15 −10 −5 Harmonics log 10 P l ( ω ) 0 5 10 15 20 25 30 −25 −20 −15 −10 −5 Harmonics log 10 P( ω ) (c) (d) Figure 1. The simulation of laser-driven hydrogen in the multiphoton ionization regime. The laser wa v elength is 880 nm, and the laser intensit y of I 0 = 10 13 W/cm 2 , corresp onding to the Keldysh parameter of γ K = 3 . 07. (a) The laser profile. (b) The induced dip ole momen t d L ( t ). (c) The p ow er sp ectrum of the laser profile. (d) The p ow er sp ectrum of d L ( t ), resp ectively . Note that the laser profile and d L ( t ) are very similar, yet very different in their sp ectral comp onents. of MWT, the same comp onen ts are exp ected to be captured but the representation should b e 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 b elow the 7th harmonic is impro v ed compared with Fig. 2(b), while harmonics on the upp er TF plane remain broaden. In addition, the frequency-dep enden t w eighting 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 b etw een neigh b oring harmonics in the upp er TF plane. The crescent shap e distribution near the b oundaries of the TF plane is caused by the b oundary effect in CWT – for a low frequency , the corresp onding scale is large and the influence of the window cut-off near the b oundary b ecomes apparen t. The b oundary effect is also inevitable in the GT, like the vertical artificial pattern in the very b eginning and end. How ever it is less dominan t. As discussed ab o ve, while the window is av oided in the WVD, a strong interference pattern is inevitable [31]. The TF representation of the WVD, sho wn in Fig. 2(c), reveals components b et w een the odd harmonics, 10 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU whic h is inconsistent with the features in the F ourier sp ectrum and violates the parity symmetry in physics. Moreo ver, harmonics ab o ve 15, of which intensities are weak, are not observed within the range of colorbar. In terference along the time-direction is also observed. As there are more oscillatory components and more o verlaps caused by x ( s + ζ / 2) x ∗ ( s − ζ / 2) in Eq. (8) as time increased, interference b ecomes stronger, as sho wn on the righ t-half TF plane in Fig. 2(c). Despite generating inaccurate information, note that in this sim ulation the WVD suggest the features of the TF represen tation, i.e., slow-v arying comp onents, regardless of a window parameter. This suggest that the window width in the linear type TF m ethods should b e large enough such that patterns depict the harmonics as that p erformed by the WVD. In order to remov e the interference, we emplo y 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, temp oral interference cannot b e fully remov ed. The in troduced filters also w orsen both temp oral and frequency resolutions, as well as generating unexp ected features. Note that the dynamic range (colorbar) is increased to rev eal weak details. By reassigning the energy distribution of the linear t yp e TF methods and the SPWVD to their lo cal cen troids, the TF representations are sharp ened. 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 significan tly impro ved the resolution, rev ealing clear and distinct o dd harmonics, as well as their subtle v ariation. 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 w ith filtering functions are remov ed and the TF repre- sen tation depicts o dd harmonics similar to those in Fig. 3(a), the interference artifact cannot b e eradicated. W e further employ the SST on the GT (SST-GT) and MWT (SST-MWT), in which the allo cation rule is a go od approximation of the IF of eac h IMT. The results sho wn 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, weigh ted frequencies, and the b oundary effect of the CWT. In Fig. 3(a), (c) and (d), we show that the IMT functions can b e decomposed without am biguit y 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 lo cal centroids, there is no routine for reconstructing the harmonics. SST, on the other hand, comes equipped with an inv erse metho d (See Eq. (23) and Eq. (26)). T o explore the physical meaning of the shifting in Fig. 4(a)-(c), we analyze it using the Flo quet metho d [53], which has b een extensively used in chemical physics [54, 55]. Details of the Flo quet metho d can b e found in App endix 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 , resp ectively , computed by the Flo quet metho d, in the unit of ω 0 . Due to the breaking of the spherical symmetry by the laser field, the energy lev els of 2s, 2p x , 2p y , and 2p z shift and split. The quasi-energies derived b y the Flo quet computation and the curve by the SST-GT are sup erimposed for comparison. W e found that the sp ectral 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 b eginning cycles of the 7th harmonic, which could b e the AC Stark effect. In the following con text, w e discuss ab out whether the IF comp onen ts obtained by the reassigned v alues hav e a ph ysical meaning, or they are simply clamped by artificial pro cessing. 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 differen t from from that in Fig. 4(b) and (c), since in the RM-GT the squared TF represen tation 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, whic h corresp onds to the AC Stark effect. The 7 . 241-th harmonic at the b eginning cycles corresp onds to the energies for 1s-2p transition ( 1 2 (1 − 1 2 2 ) /ω 0 = 7 . 241 in a.u.), manifested as a small p eak ov erlapping the H7 of the p ow er sp ectrum in Fig. 1(d). The intensit y of such sp ectral line is small b ecause it arises from the near resonance absorption (not resonance absorption). 11 Time (Cycle) Harmonic Order 0 20 40 60 0 5 10 15 20 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 Time (Cycle) Harmonic Order 0 20 40 60 0 5 10 15 20 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 (a) (b) Time (Cycle) Harmonic Order 0 20 40 60 0 5 10 15 20 −9 −8 −7 −6 −5 −4 −3 −2 −1 0 Time (Cycle) Harmonic Order 0 20 40 60 0 5 10 15 20 −14 −12 −10 −8 −6 −4 −2 0 (c) (d) Figure 2. TF represen tations in the multiphoton ionization regime by the (a) GT, (b) MWT, (c) WVD, and (d) SPWVD. The lines indicating odd harmonics in the GT and MWT are sub ject to broadening issues arising from a window. Although the broadening artifact can b e alleviated by the WVD, the evok ed interference artifacts in the TF representation result in incomprehensible analysis. By applying additional filters, the SPWVD can mo derate the in terference patterns, at the price of broadening of the features in the TF representation. Fig. 4(b) and (c) suggest that the IFs shown b y different SST metho ds are indep enden t of the c hosen linear TF metho d. Based on the indep endence of the c hosen SST metho ds and the in timate matching of the analysis result and the theoretical prediction shown in Fig. 4(d), we conclude that the decomp osed IMTs are physically meaningful. The ionization pro cess for the laser field in Fig. 1(a) is summarized in Fig. 5. In Fig. 5(a), when the laser field is small, the sp ectral line 7 . 241 aroused from the 1s-2p near resonance absorption caused by the atomic structure. As the laser intensit y gradually increases, the energy levels of 2s, 2p x , 2p y , and 2p z shift and split due to the breaking of the symmetry b y the electric field, whic h is regarded as the AC Stark effect. When the laser intensit y increases furthermore, the high-order harmonic pro cess is induced and the dynamics is dominated by the transition b etw een dressed states | M , N i formed b y the electron state M and the photon state N , as illustrated in Fig. 5(b). 12 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU Time (Cycle) Harmonic Order 0 20 40 60 0 5 10 15 20 −18 −16 −14 −12 −10 −8 −6 −4 −2 0 Time (Cycle) Harmonic Order 0 20 40 60 0 5 10 15 20 −14 −12 −10 −8 −6 −4 −2 0 (a) (b) Time (Cycle) Harmonic Order 0 20 40 60 0 5 10 15 20 −9 −8 −7 −6 −5 −4 −3 −2 −1 0 Time (Cycle) Harmonic Order 0 20 40 60 0 5 10 15 20 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 (c) (d) Figure 3. TF representations in the multiphoton ionization regime by the (a) RM-GT and (b) RM-SPWVD. The broadening caused by the window in the GT and filter functions in SPWVD is remo ved. Note that the frequency shift in the beginning 10 cycles corresponds to the AC Stark effect. TF representations in the multiphoton ionization regime by the (c) SST-GT and (d) SST-MWT. The broadening issue caused by the window in GT and MWT are remov ed. Note that the frequency shift in the b eginning 10 cycles corresp onds the AC Stark effect. A summary for comparison of TF metho ds in the multiphoton ionization regime is provided in T able 2. Note that in this table we mention that the o dd harmonics and the AC Stark effect cannot b e obtained by EMD with the HS. In fact, we could not yet find ph ysical meaning for the series of IMFs decomp osed by the EMD in the multiphoton ionization regime. 4.2. T unneling Ionization Regime. In the second sim ulation, we discuss the TF represen tations for signals in the tunneling ionization regime. The laser wa velength in this subsection is ω 0 ≈ 0 . 05696100 in atomic units (a.u.), corresp onding to 800 nm, and the laser intensit y 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 13 Time (Cycle) Harmonic Order 2 4 6 8 10 7.19 7.2 7.21 7.22 7.23 7.24 7.25 7.26 −14.5 −14 −13.5 −13 −12.5 −12 −11.5 −11 −10.5 −10 Time (Cycle) Harmonic Order 2 4 6 8 10 7.19 7.2 7.21 7.22 7.23 7.24 7.25 7.26 −6.5 −6 −5.5 −5 −4.5 −4 −3.5 −3 −2.5 −2 (a) (b) Time (Cycle) Harmonic Order 2 4 6 8 10 7.19 7.2 7.21 7.22 7.23 7.24 7.25 7.26 −6.5 −6 −5.5 −5 −4.5 −4 −3.5 −3 −2.5 −2 (c) (d) Figure 4. The A C Stark effect revealed by the (a) RM-GT, (b) SST-GT and (c) SST- MWT. (d) Comparison of the frequency shift caused by the A C Stark effect computed by the SST-GT (gray scale) and the Flo quet metho d. T able 2 Comparison of TF metho ds in the m ultiphoton ionization regime Phenomenon Discrete Odd Harmonics The AC Stark Effect GT Y es (blurred) Y es (blurred) MWT Y es (blurred) Y es (blurred) WVD No No SPWVD Y es, for Y es (blurred) harmonics > 5 (blurred) RM-GT Y es Y es RM-SPWVD Y es Y es (blurred) SST-GT Y es Y es SST-MWT Y es Y es EMD-HS No No In this study , the GT and MWT are examples for the STFT and CWT, resp ectively . 14 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU (a) (b) Figure 5. Illustration for the (a) AC-Stark effect and (b) the high-order harmonic pro cess. and its profile is ramp ed on according to E 0 ( t ) =    sin 2 ( π t 10 T ) when 0 ≤ t ≤ T 1 when T ≤ t ≤ 9 T sin 2 ( π t 10 T ) when 9 T ≤ t ≤ 10 T , (30) whic h leads to the laser field E 0 F ( t ) sin( ω 0 t ), as shown in Fig. 6(a). W e adopt the dip ole momen t in acceleration form in the case of tunneling ionization, as display ed in Fig. 6(b). The acceleration form can pro vide clearer information for the high order harmonics than the length form, b ecause the differen tial op erator 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 p eriodic rep etition of arches, suggesting the chirp-lik e dynamics of the attosecond radiation, as predicted by the standard semiclassical mo del [29]. T o present such feature in the linear t yp e TF representations, a short windo w width and a small quality factor τ are needed for the GT and the MWT, resp ectively . (Applicability of a short window width in the IMT mo dels deserves further inv estigation.) 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 rep eats every 0 . 5T. The inner structure after the 2nd cycle is indistinct b ecause of the broadening caused by the window. While the GT cannot resolve information b elo w the 15th harmonic due to the small windo w width, the MWT can describ e the near and the b elow threshold harmonics (low order harmonics) b y 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. Ho wev er, as describ ed in previous subsection, WVD introduces strong in terferences in b oth temp oral and frequency directions. Despite in SPWVD, where we use the parameters σ g = 1 and σ H = 0 . 13, it do es not seem p ossible to eliminate the in terference that mingled with intrinsic structure inside the main arches. W e further apply the RM techniques b y 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, representativ es for the linear and quadratic type metho ds, resp ectively , can address the broadening issue in the original metho ds and provide distinct inner arc hes 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 in terference that is difficult to remov e. It is b ecause that the RM techniques re-shuffle the original TF representation based on its lo cal center of mass, and hence preserves the intrinsic artifact in the adopted transform. By a reallo cation rule that contains phase information of a linear t yp e TF transform, the SST can enhance the TF represen tation b y linear t yp e TF transforms in a manner similar to that by the RM tec hnique, 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 15 0 2 4 6 8 10 −0.1 −0.05 0 0.05 0.1 Time (Cycle) Laser Field (atomic unit) 0 2 4 6 8 10 −0.04 −0.03 −0.02 −0.01 0 0.01 0.02 0.03 0.04 Time (Cycle) d A (t) (atomic unit) (a) (b) Figure 6. The simulation of laser-driv en hydrogen in the tunneling ionization regime. The laser wa velength is 800 nm, and the laser intensit y of 3 . 5 × 10 14 W/cm 2 , corresp onding to the Keldysh parameter of γ K = 0 . 57. (a) The laser profile. (b) The induced dip ole momen t d A ( t ). fixed resolution feature in the STFT, and the SST-MWT maintains the adaptiv e resolution in the CWT. By comparing Fig. 8(c) with Fig. 8(c), we see that the resolution in the RM-GT is more sharp en than that of the SST-GT. One of the reasons is that while b oth temp oral and frequency reassignments are considered in the conv entional reassignmen t metho d, temp oral reassignment is not taken in to account in the SST (b oth SST-GT and SST-MWT). In addition, the reallo cation rule Eq. (20) is only the first order approximation for the IF, resulting in diffusiv e pattern particularly for fast-v arying chirp signal. F or the simulation in the multiphoton ionization regime, such diffusive pattern is negligible as the IF changes slowly . By using Eq. (21), a second-order reallo cation rule for the SST-GT, the diffusive pattern can b e further concentrated, as is shown in Fig. 8(e). T o shed some light on this TF represe n tation for the tunneling ionization mec hanism, we follow the standard semiclassical approach suggested indep endently b y Corkum [56] and Kulander et al. [57]. Here the electric-field force corresp onding to the applied laser field in a.u. is F z = − E 0 E ( t ) ˆ z , where ˆ z is the unit v ector in the z − direction. W e assume the initial p osition and initial velocity is b oth zero [29]. The tra jectory of electrons released b etw een 1 T and 1 . 5 T by using the extended semiclassical approach is denoted by red circles in Fig. 8(f ). This result is sup erimposed for the sake of comparison with the TF represen tation of the SST-GT. The tra jectory suggests that the largest arch in Fig. 8 is asso ciated with the first return time of the electron w a v e pack et released in betw een 1 T and 1 . 5 T of the laser field, and the smaller arch b etw een 2 T and 2 . 5 T , is the second return. Although the SST depicts the first and multiple returns clearly , the sequence of these returns cannot b e kno wn. T o retrieve the first and second return quantitativ ely , Risoud et al. prop ose the application of the Hilb ert transform on the dip ole moment with the harmonics low er than the ionization p otential filtered [18]. While there are more than one comp onent simultaneously exist at each time, we consider the EMD scheme to decomp ose the filtered dip ole moment. (In this case, harmonics below the 10th harmonic are filtered as they corresp ond to the b ound part of the atomic sp ectrum.) The stopping criterion in the EMD scheme is that the residue b ecomes monotonic. The HS for all mo des of IMF are display ed in Fig. 9(a). Note that a median filter with a windo w length of 51 p oints is applied on all mo des as the implementation of Eq. (28) can introduce n umerical errors. Despite there is no theoretical base for the IMFs by the EMD and the result is not consistent in the m ultiphoton ionization case, in this case we see that the first tw o mo des can be asso ciated with the first and second returns. Fig. 9(b) presen ts the HS for the first tw o IMFs and Fig. 9(c) compares the IMFs with the representation of the SST-GT. In Fig. 9(b), the IMFs are visually emphasized 16 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −5 −4.5 −4 −3.5 −3 −2.5 −2 Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −4 −3.5 −3 −2.5 −2 −1.5 −1 (a) (b) Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −9 −8 −7 −6 −5 −4 −3 −2 Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −7 −6.5 −6 −5.5 −5 −4.5 −4 −3.5 −3 (c) (d) Figure 7. TF represen tations in the tunneling ionization regime by the (a) GT, (b) MWT, (c) WVD, and (d) SPWVD. The rep eating arches suggest excitation of the harmonics is strongly correlated with the laser field. In the representations of the GT and MWT, the arc hes are broaden by the window, leading to the difficulty to differen tiate the in tricate structures inside the arches. The representation of the WVD is obscured b y the interference pattern, and it is difficult to remov e the interference pattern with filters in the SPWVD. b y additional Gaussian functions in b oth temp oral and frequency domain. Note that the IF rev ealed by the Hilb ert transform is meaningful only if the IMF is consist of a single comp onent, which is not guaranteed in the EMD. F or example, for time greater than 2 . 5 cycles, the HS of the first IMF is distorted because of more wa v e pac ket returning from previous emissions. While the Hilbert transform is not theoretically suitable for this kind of oscillatory signal with time-v arying AM and IF, w e illustrate the IF of these decomp osed mo des b y the SST-GT. W e employ the SST-GT on the first tw o 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 Hilb ert transform. Note that no median filter is necess ary for each IMFs in the SST-GT analysis, which reduces the p ossible artifacts. A summary for comparison of TF metho ds in the tunneling ionization regime is pro vided in T able 3. 17 Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −6 −5 −4 −3 −2 −1 0 Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −6 −5.5 −5 −4.5 −4 −3.5 −3 (a) (b) Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −3 −2.5 −2 −1.5 −1 Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −3 −2.5 −2 −1.5 −1 (c) (d) Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −3 −2.5 −2 −1.5 −1 1st return 2nd return (e) (f ) Figure 8. TF represen tations in the tunneling ionization regime by the (a) RM-GT and (b) RM-SPWVD. The broadening artifacts caused by the window in b oth metho ds are eliminated, but it is difficult to remo ve the in terference pattern in the SPWVD. TF represen tations in the tunneling ionization regime by the (c) SST-GT and (d) SST-MWT. The broadening artifacts caused by the window in b oth metho ds are alleviated. (e) TF represen tation in the tunneling ionization regime by the second order SST-GT. Note that the diffusive pattern in the SST-GT is mo dified for the arc hes. (f ) Comparison of the TF represen tation of the SST-GT and a tra jectory of the electron released b et w een the 1 T and 1 . 5 T (denoted by red circles) computed by the standard semiclassical approach suggested indep enden tly by Corkum [56] and Kulander et al. [57]. 18 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −7 −6 −5 −4 −3 −2 Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −6 −5 −4 −3 −2 −1 0 (a) (b) (c) Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −4 −3 −2 −1 0 Time (Cycle) Harmonic Order 1.5 2 2.5 3 3.5 0 10 20 30 40 50 60 −4 −3 −2 −1 0 (d) (e) Figure 9. (a) TF representation in the tunneling ionization regime of the EMD-HS algorithm, whic h con tains all IMF modes extracted b y EMD. (b) TF represen tation in the tunneling ionization regime of the first t wo IMFs extracted by EMD. T o enhance the visualization, the curves are thick ened. (c) Comparison of the first t w o IMFs with the TF represen tation by the SST-GT. According to the semiclassical tra jectory , the first IMF is the first return and the second IMF is the second return. (d) The first IMF and (e) the second IMF analyzed by the SST-GT. 19 T able 3 Comparison of TF metho ds in the tunneling ionization regime Phenomenon First Return Second Return Sequencing of Near and Below Multiple Returns Threshold Harmonics GT Y es (blurred) Y es (blurred) No No MWT Y es (blurred) Y es (blurred) No Y es (blurred) WVD No No No No SPWVD Y es (blurred) No No No RM-GT Y es Y es No Y es RM-SPWVD Y es No No No SST-GT Y es (diffused) Y es No No 2nd order SST-STFT Y es Y es No No SST-MWT Y es (diffused) Y es No Y es EMD-HS Y es (mo de mixing) Y es (mo de mixing) Y es No EMD-SST Y es (mo de mixing) Y es (mo de mixing) Y es No 5. Conclusions This pap er provides a b enc hmark of analyzing time-dep endent quan tum systems by applying several con- temp orary TF methods. Within the same physical mo del and computational sc heme, different features of TF represen tations pro vided by different TF metho ds ma y lead to conflicting interpretation in physics. In the m ultiphoton ionization regime, linear TF metho ds with an appropriate window function, can v aguely depict discrete o dd harmonics in the HHG pro cess, yet the detail features such as the AC Stark effect cannot b e rev ealed. In addition, in the tunneling ionization regime, the first-return tra jectories are roughly illustrated, but the multiple returns are indistinct because of the broadening induced by the windo w. While present little broadening, TF representations b y the WVD suffer from artificial interference patterns and lead to incorrect in terpretation of physical mec hanisms, suc h as even harmonics in the m ultiphoton ionization regime and false tra jectories in the tunneling ionization regime. T o eliminate the artifacts mentioned ab o v e, several mo dern TF analysis metho ds are introduced in quantum dynamics for the first time. How ev er, Cohen class distribu- tions suc h as the SPWVD can only remo v e p ortion of in terference patterns. Reassignment techniques and the SST can depict accurate dynamic pro cess in the tw o regimes based on the principle of lo cal centroid and instan taneous frequency , resp ectively . How ev er, there is no inv erse transform in the reassignmen t techniques. On the other hand, the SST preserv es signal causality and allows mo de reconstruction. W e demonstrate that after removing the broadening, the AC Stark effect in the multiphoton ionization regime and the m ultiple returns in the tunneling ionization regime are rev ealed. Finally , we relate the IMFs from the EMD to the m ultiscattering of the electron energy distribution in the tunneling ionization regime. W e b elieve that in addition to the atomic hydrogen system, the contemporary TF metho ds hav e p oten tial for exploring other complicated quantum systems. In case of analyzing dynamics in an unknown quantum system, this pap er suggests the following pro cedure. First, features obtained by the WVD can b e used to determine the parameter for the linear type metho ds. Then the SST can b e applied to obtain TF represen- tations for further in terpretation. W e hop e that this research can serve as a cornerstone in applications of TF analysis for fields including but not limited to, attosecond ph ysics, n uclear magnetic resonance, ultrafast dynamics in atoms and molecules. 6. A cknowledgments The authors would like to thank Dr. Elise Y. Li and Y u-lin Sheu for their help in servers. 20 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU Appendix A. Adaptive harmonic model Fix constan ts 0 ≤   1, 0 < d < 1, c 2 > c 1 >  and c 2 > c 3 >  . Consider the functional set Q c 1 ,c 2 ,c 3  , whic h consists of functions x ( t ) = A ( t ) cos(2 π φ ( t )) ∈ C 1 ( R ) ∩ L ∞ ( R ) where the following conditions hold:              A ∈ C 1 ( R ) ∩ L ∞ ( R ) , φ ∈ C 3 ( R ) , inf t ∈ R A ( t ) ≥ c 1 , inf t ∈ R φ 0 ( t ) ≥ c 1 , sup t ∈ R A ( t ) ≤ c 2 , sup t ∈ R φ 0 ( t ) ≤ c 2 , sup t ∈ R | φ 00 ( t ) | ≤ c 3 , | A 0 ( t ) | ≤ φ 0 ( t ) , | φ 000 ( t ) | ≤ φ 0 ( t ) for all t ∈ R , W e call x = A ( t ) cos(2 π φ ( t )) ∈ Q c 1 ,c 2 ,c 3  an intrinsic mo de typ e (IMT) function, φ ( t ) the phase function, the deriv ative of the phase function φ 0 k ( t ) is called the instantane ous fr e quency (IF) and A ( t ) is called the amplitude mo dulation . When | φ 00 ( t ) |  φ 0 ( t ) for all t ∈ R , we call the signal x oscillatory with slow ly varying IF ; otherwise we call it oscillatory with fast varying IF . Lo cally an IMT function with slowly v arying IF b ehav es like a harmonic function and an IMT function with fast v arying IF b eha ves 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 prop erties defined in the IMF, but the reverse inclusion do es not hold. T o describ e an oscillatory signal, the follo wing adaptive harmonic model is commonly considered in the time-frequency analysis literature [38, 39, 40]. Fix constants 0 ≤   1, d > 0 and c 2 > c 1 > 0. Consider the functional set Q c 1 ,c 2 ,c 3 ,d , which consists of functions in C 1 ( R ) ∩ L ∞ ( R ) with the following format: x ( t ) = P K ` =1 x k ( t ), where K is finite and x k ( t ) = A ` ( t ) cos(2 π φ ` ( t )) ∈ Q c 1 ,c 2 ,c 3  . Let φ 0 k ( t ) > φ 0 k − 1 ( t ). In the STFT-SST, success iv e IMTs for all time t hav e to b e separated by at least d , i.e., φ 0 k ( t ) − φ 0 k − 1 ( t ) > d . Appendix B. The Floquet Method The total Hamiltonian of atomic hydrogen can b e expressed as H ( t ) = H 0 + H 1 ( t ), where H 0 is the unp erturbed hydrogen Hamiltonian and H 1 ( t ) = − z E ( t ) describ es the laser-atom interaction. In the time- dep enden t generalized pseudosp ectral simulation, the time p erio d T is fixed but the field amplitude v aries. T o calculate the shifting of orbital energies, we separate each cycle and take its field amplitude, and then assume the h ydrogen in a laser field with the particular constant amplitude. F or example, for the third cycle in Fig. 1(a), the maximal amplitude is 4 . 86 × 10 − 4 , and w e assume hydrogen in the laser field with the constan t amplitude. When the Hamiltonian is p erio dic with a time p erio d T , i.e., H ( t ) = H ( t + T ), we can apply the Flo quet theory [53] and derive a time-indep endent infinite-dimensional eigenv alue matrix equation [54] as follows: X nlmk hh nlmk | b H F | n 0 l 0 m 0 k 0 ii φ n 0 l 0 m 0 k 0 λξ = ε λξ φ nlmk λξ , (31) Here hh nl mk | b H F | n 0 l 0 m 0 k 0 ii denotes the matrix elements of the time-av eraged Flo quet Hamiltonian ov er a p eriod T , where the outer bra-ket notation refers to the inner pro duct o v er t . W e use the orbitals of an unp erturbed hydrogen as a basis set since we only consider a weak field regime. As a result, n , l , and m represent the principal quan tum num b er, angular momentum quantum n umber, and magnetic quantum n umber of h ydrogen, resp ectiv ely . k is an integer. The matrix element in Eq. (31) can b e ev aluated as hh nlmk | b H F | n 0 l 0 m 0 k 0 ii = E nlm δ nn 0 δ ll 0 δ mm 0 δ kk 0 − h nlm | z | n 0 l 0 m 0 i eE 0 2 ( δ k,k 0 +1 + δ k,k 0 − 1 ) + k ~ ω δ nn 0 δ ll 0 δ mm 0 δ kk 0 (32) The v alues of h nlm | z | n 0 l 0 m 0 i can b e computed analytically . F or examples, in atomic units, h 210 | z | 1 0 0 0 0 0 i = 0 . 7449, h 310 | z | 1 0 0 0 0 0 i = 0 . 2983, h 210 | z | 3 0 0 0 0 0 i = 0 . 5418, and h 310 | z | 2 0 0 0 0 0 i = 1 . 7695. 21 References [1] T. F uji, T. Saito, and T. Kobay ashi, “Dynamical observ ation of Duschinsky rotation by sub-5-fs real-time sp ectroscopy ,” Chem. Phys. Lett. 11 (2), 99–116 (2011). [2] T. Kobay ashi and A. Y abushita, “T ransition-state sp ectroscopy using ultrashort laser pulses,” Chem Rec. 11 (2), 99–116 (2011). [3] M. J. J. V rakking, D. M. Villeneuve, and A. Stolo w “Observation of fractional revivals of a molecular wa ve pac ket,” Ph ys. Rev. A 54 (1), R37– (1996). [4] L. W ang, W. Liu and C. F ang “Elucidating low-frequency vibrational dynamics in calcite and water with time-resolved third-harmonic generation spectroscopy ,” Ph ys. Chem. Chem. Phys. 17 (26), 17034 (2015). [5] A. V olpato and E. Collini, “Time-frequency methods for coherent spectroscopy ,” Opt. Express 23 (15), 20040–20050 (2015). [6] J. Prior, E. Castro, A. W. Chin, J. Almeida, S. F. Huelga, and M. B. Plenio, “W av elet analysis of molecular dynamics: Efficient extraction of time-frequency information in ultrafast optical pro cesses,” J. Chem. Phys 139 (22), 224103 (2013). [7] Y. Deng, Q. Sun, F. Liu, C. W ang, and Q. Xing “T erahertz Time-resolved Spectroscopy with W a velet-transform,” in International Congress on Image and Signal Pro c essing (IEEE, 2010), pp. 3462–3464. [8] C. Y. W ong, R. M. Alvey , D. B. T urner, K. E. Wilk, D. A. Bryan t, P . M. G. Curmi, R. J. Silb ey , and G. D. Scholes, “Electronic coherence l ineshapes rev eal hidden excitonic correlations in photosyn thetic light harv esting,” Nat. Chem. 4 (5), 396–404 (2012). [9] S. D. McClure , D. B. T urner , P . C. Arpin , T. Mirko vic , and G. D. Scholes, “Coherent oscillations in the PC577 Cryptophyte antenna occur in the excited electronic state,” J. Phys. Chem. 118 (5), 1296–1308 (2014). [10] D. B. T urner, R. Dinshaw, K.-K. Lee, M. S. Belsley , K. E. Wilk, P . M. G. Curmic and G. D. Scholes, “Quantitativ e inv estigations of quantum coherence for a light-harv esting protein at conditions simulating photosynthesis,” Ph ys. Chem. Chem. Phys. 14 , 4857–4874 (2012). [11] F. Krausz and M. Iv anov, “Attosecond physics,” Rev. Mod. Phys. 81 (1), 163–234 (2009). [12] C. C. Chirila, I. Dreissigack er, E. V. v an der Zwan, and M. Lein, “Emission times in high-order harmonic generation,” Phys. Rev. A 81 (3), 033412 (2010). [13] C. Figueira de Morisson F aria, M. Dorr, W. Beck er, and W. Sandner, “Time-frequency analysis of t wo-color high-harmonic generation,” Phys. Rev. A 60 (2), 1377–1384 (1999). [14] P . An toine, B. Piraux, D. B. Milosevic, and M. Ga jda, “Generation of ultrashort pulses of harmonics,” Ph ys. Rev. A 54 (3), R1761–R1764 (1996). [15] X. M. T ong and Shih I. Chu, “Probing the sp ectral and temp oral structures of high-order harmonic generation i n in tense laser pulses,” Phys. Rev. A 61 (2), 021802(R) (2000). [16] X. Ch u, S. I. Chu, and C. Laughlin, “Sp ectral and temporal structures of high-order harmonic generation of Na in intense mid-ir laser fields,” Phys. Rev. A 64 (1), 013406 (2001). [17] J. Chen, S. I. Chu, and J. Liu, “Time-frequency analysis of molecular high-harmonic generation sp ectrum by means of wa velet transform and Wigner distribution tec hniques,” J. Phys. B: At. Mol. Opt. Phys. 39 (22), 4747–4758 (2006). [18] F. Risoud, J. Caillat, A. Maquet, R. T aieb, and C. Leveque, “Quantitativ e extraction of the emission times of high-order harmonics via the determination of instantaneous frequencies,” Phys. Rev. A 88 (4), 043415 (2013). [19] P .-C. Li, Y.-L. Sheu, C. Laughlin, and S.-I Chu, “Role of laser-driven electron-multirescattering in resonance-enhanced below-threshold harmonic generation in He atoms,” Phys. Rev. A 90 (4), 041401(R) (2014). [20] P .-C. Li, Y.-L. Sheu, C. Laughlin, and S.-I Chu, “Dynamical origin of near- and b elow-threshold harmonic generation of Cs in an intense mid-infrared laser field,” Nat. Comm. 6 , 7178 (2015). [21] S. Phillips, R. J. Gledhill, J. W. Essex, and C. M. Edge, “Application of the Hilb ert-Huang T ransform to the Analysis of Molecular Dynamics Simulations,” J. Phys. Chem. A 107 (24), 4869–4876 (2003). [22] M. Orini, R. Bailon, L. T. Mainardi, P . Laguna, and P . Flandrin, “Characterization of Dynamic Interactions Betw een Cardiov ascular Signals by Time-F requency Coherence,” IEEE T rans. Biomed. Eng 59 (3), 663–673 (2012). [23] Y. F. Suprunenko, P . T. Clemson, and A. Stefanovsk a, “Chronotaxic Systems: A New Class of Self-Sustained Nonau- tonomous Oscillators,” Phys. Rev. Lett. 111 (2), 024101 (2013). [24] M. Niethammer, L. J. Jacobs, J. Qu, and J. Jarzynski, “Time-frequency representations of Lamb w aves,” J. Acoust. So c. Am. 109 (5 Pt 1), 1841–1847 (2001). [25] P . W ang, J. Gao, and Z. W ang, “Time-F requency Analysis of Seismic Data Using Synchrosqueezing T ransform,” IEEE Geosci. Remote Sens. Lett. 11 (12), 2042–2044 (2014). [26] A. Ludwig, J. Maurer, B. M. May er, C. R. Phillips, L. Gallmann, U. Keller, “Breakdown of the Dipole Approximation in Strong-Field Ionization,” Phys. Rev. Lett., 113 (24), 243001 (2014). [27] C. Y ao and S. I. Chu, “Generalized pseudosp ectral methods with mappings for b ound and resonance state problems,” Chem. Phys. Lett., 204 (3-4), 381–388 (1993). [28] X. M. T ong and S. I. Chu, “Theoretical study of m ultiple high-order harmonic generation by in tense ultrashort pulsed laser fields: A new generalized pseudosp ectral time-dependent method,” Chem. Phys., 217 (2-3), 119–130 (1997). [29] Z. Chang, F undamentals of Attose c ond Optics (CRC Press, 2011). [30] M. Lewenstein, P . Balcou, M. Y. Ivano v, A. L’Huillier, and P . B. Corkum, “Theory of high-harmonic generation by low-frequency laser fields,” Phys. Rev. A, 49 (3), 2117–2132 (1994). [31] P . Flandrin, Time-fr e quency/time-sc ale Analysis, Wavelet Analysis and Its Applications (Academic Press Inc., 1999). [32] I. Daub ec hies, T en L e ctur es on Wavelets (Philadelphia: So ciet y for Industrial and Applied Mathematics, 1992). 22 Y AE-LIN SHEU, HAU-TIENG WU, AND LIANG-Y AN HSU [33] J. C. O’Neill and P . Flandrin, “Virtues and Vices of Quartic TimeF requency Distributions,” IEEE T rans. Signal Process., 48 (9), 2641–2650 (2000). [34] B. Boashash and B. Ristich , “Polynomial Wigner-Ville distributions and time-varying higher-order sp ectra,” in Time- F r e quency and Time-Sc ale Analysis , (Proc. IEEE-SP Int. Symp., 1992), pp. 31–34. [35] P . Flandrin, F. Auger, and E. Chassande-Mottin, Time-F r e quency Signal Pr o c essing (CRC Press, 2003). [36] F. Auger, and P . Flandrin, “Improving the Readability of Time-F requency and Time-Scale Representations by the Reas- signment Metho d,” IEEE T rans. Signal Pro cess., 43 (5), 1068–1089 (1995). [37] F. Auger, and P . Flandrin, Y. T. Lin, S. McLaughlin, S. Meignen, T. Ob erlin, and H. T. W u, “Time-F requency Reassignment and Synchrosqueezing,” IEEE Signal Pro cessing Mag., 30 (6), 32–41 (2013). [38] I. Daubechies, J. Lu, H. T. W u, “Synchrosqueezed W av elet T ransforms: An empirical mode decomp osition-like to ol,” Appl. Comput. Harmon. Anal., 30 (1), 243–261 (2011). [39] H. T. W u, A daptive analysis of c omplex data sets Do ctoral dissertation (Princeton Universit y , 2011). [40] Y.-C. Chen, M. Y. Cheng, H. T. W u, “Nonparametric and adaptive mo deling of dynamic p erio dicit y and trend with heteroscedastic and dependent errors,” J. R. Stat. Soc. B, 76 (3), 651–682 (2014). [41] T. Oberlin, S. Meignen, V. P errier, “Second-Order Sync hrosqueezing T ransform or In vertible Reassignmen t? T ow ards Ideal Time-F requency Representations,” IEEE T rans. Signal Pro cess., 63 (5), 1335–1344 (2015). [42] L. Cohen, Time-fr e quency Analysis, Signal Pr o c essing Series (Prentice-Hall, 1995). [43] D. Gab or, “Theory of comm unication,” Pro c. IEE, 93 , 429–457 (1946). [44] N. E. Huang, Z. Shen, S. R. Long, M. C. W u, H. H. Shih, Q. Zheng, N. C. Y en, C. C. T ung, H. H. Liu, “The Empirical Mode Decomp osition and the Hilb ert Spectrum for Nonlinear and Nonstationary Time Series Analysis,” Pro c. R. Soc. A, 454 (1971), 903–995 (1998). [45] G. Rilling and P . Flandrin, “One or Two F requencies? The Empirical Mo de Decomp osition Answ ers,” IEEE T rans. Signal Process., 56 (1), 85–95 (2008). [46] S. H. Lin, A dvanc es in Multi-Photon Pr o c esses and Spe ctr osc opy, V ol. 19. (W orld Scientific Publishing, 2010). [47] K. Mishima, K. Nagay a, M. Hay ashi and S. H. Lin, “T ow ards the realization of the quantum chemistry approach to tunneling photoionization processes in strong laser fields,” J. Chem. Phys., 122 (2), 024104 (2005). [48] H. Mineo, K. Nagay a, M. Hay ashi and S. H. Lin, “Theoretical studies of high-harmonic generation based on the Keldysh- F aisalRei ss theory ,” J. Phys. B: A t. Mol. Opt. Phys., 40 (12), 2435–2451 (2007). [49] Y. Nomura, Y. F ujimura and H. Kono, “Theory of quantum b eats in timeresolved multiphoton ionization of molecules,” J. Chem. Phys., 88 (3), 1501–1510 (1988). [50] M. Gavrila, A toms in Intense L aser Fields. (Academic Press, 1992). [51] P . Gibbon, Short Pulse Laser Inter actions with Matter (Imp erial College Press, 2005). [52] E. P . Po wer, A. M. March, F. Catoire, E. Sistrunk, K. Krushelnick, P . Agostini and L. F. DiMauro, “XFROG phase measurement of threshold harmonics in a Keldysh-scaled system,” Nature Photon., 4 (6), 352–356 (2010). [53] J. H. Shirley , “Solution of the sc hrdinger equation with a hamiltonian perio dic in time,” Phys. Rev. B, 138 (4B), B979–B987 (1965). [54] L. Y. Hsu, D. Xie, and H. Rabitz, “Light-driv en electron transp ort through a molecular junction based on cross-conjugated systems,” J. Chem. Phys., 141 (12), 124703 (2014). [55] L. Y. Hsu and H. Rabitz, “Coherent light-driv en electron transp ort through polycyclic aromatic hydrocarbon: laser fre- quency , field intensit y , and p olarization angle dependence,” Phys. Chem. Chem. Ph ys., 17 (32), 20617–20629 (2015). [56] P . B. Corkum, “Plasma persp ective on strong field m ultiphoton ionization,” Ph ys. Rev. Lett., 71 (13), 1994–1997 (1993). [57] K. C. Kulander, K. J. Schafer, and J. L. Krause, in Pro ceedings of the W orkshop on Super-Intense Laser Atom Physics (SILAP) II I. P . Piraux et al. , ed. (Plenum Press, 1993). 1215 Hunters Glen Dr., Plainsboro, New Jersey 08536, USA, Dep ar tment of Ma thema tics, University of Toronto, Toronto ON M5S 2E4, Canad a. Email: hauwu@ma th.toronto.edu Dep ar tment of Chemistr y, Princeton University, Princeton, New Jersey 08544, USA. Email: lianghsu@princeton.edu