Probing and controlling electron behavior in matter at attosecond timescales has become possible with attosecond pulses generated by few-cycle intense lasers. This breakthrough has revolutionized our understanding of the atomic structure and molecular processes^[1–10]. Various experiments have been carried out, claiming exceptional time resolution down to a few tens of attoseconds, using tools such as attosecond photoelectron streaking^[11,12], the reconstruction of attosecond beating by interference of two-photon transitions (RABBIT)^[13,14], high-harmonic spectroscopy^[15], attoclocks^[16,17], and attosecond transient absorption spectroscopy (ATAS)^[18,19]. When observing ultrafast electronic motion, understanding the dynamics of excited states is a fundamental issue in various laser–matter interaction processes^[20–28]. Specifically, the emission time of high-harmonic radiation from excited states is temporally delayed, providing a means to monitor excited-state dynamics^[26,29,30]. The role of laser-dressed excited energy levels in atomic excitation and ionization has been studied using attosecond technology^[4,31,32]. In the process of strong-field tunneling ionization, the transition of electrons to the continuum via multiple bound states has been investigated^[33–36].

With the advancement of pump-probe experiments and the development of attosecond pulses, it has become possible to detect information about the bound states of atoms. De Boer and Muller^[37] obtained the populations of excited states by calculating the area probability of each peak in the photoelectron energy distribution using a nanosecond probe pulse. This method requires long lifetimes of the bound states. More recently, under the high sensitivity of attosecond transient absorption spectroscopy (ATAS), the real-time population of valence states has been experimentally measured using a combination of a few-cycle near-infrared (NIR) pulse and a single XUV attosecond pulse^[38].

Quantum-beat spectroscopy in a pump-probe setting is an alternative interferometric tool for obtaining bound-state dynamics. Building upon earlier research that employed high-harmonic radiation^[39], XUV frequency combs^[40], and synchrotron radiation^[41], the coherence and amplitudes of bound states have been extracted from the emission attosecond quantum beat spectrum in hydrogen^[42], where two optical attosecond pulses are used. Quantum-beat signals carrying bound state information have also been mapped onto the photoelectron spectrum^[43–46]. A wave packet in low-lying excited states of the helium atom, oscillating with an ultrashort period of 2.0 fs, has been mapped out^[43]. Zhang et al.^[45] have reconstructed time-evolved density matrix elements in one-dimensional (1D) atomic hydrogen by exploiting femtosecond XUV pulses and a terahertz (THz) streaking field. Recently, the populations of bound states caused by an XUV pump pulse have been retrieved with high precision from the photoelectron quantum beat spectrum in the 1D hydrogen atom^[46]. In our previous work, the NIR probe pulse used was strong and long to ensure the complete ionization of electrons. However, achieving complete ionization in the experiment is challenging, making it more of a qualitative model in principle.

In this work, we propose a different setting, using strong attosecond XUV pump and weak attosecond XUV probe pulses to reconstruct the populations of bound states from the quantum-beat signals in the photoelectron spectrum. In this approach, the probe pulse is weak and short, which differs from the treatment of the strong-field approximation^[47–49] in Ref. [45]. We examine this scheme for strong attosecond XUV laser-induced ionization of hydrogen atoms by solving a three-dimensional time-dependent Schrödinger equation (TDSE). We demonstrate that the dynamics of bound states can be reconstructed with high precision. This general approach provides a valuable reference for reconstructing the ultrafast dynamics of bound states induced by strong XUV attosecond pulses.

In our scheme, we employ a strong attosecond few-cycle XUV pulse, referred to as the pump pulse in Fig. 1(a), to coherently excite and ionize the atom, as shown in Fig. 1(b). This attosecond pump pulse facilitates a dominant one-photon absorption process. A weaker attosecond XUV pulse with a time delay τd relative to the pump pulse serves as the probe pulse. The probe pulse couples the induced excited states, such as |b⟩ and |j⟩, with the continuum state |q⟩, leading to the transfer of coherent population to the same ionized state |p⟩. Consequently, the state |p⟩ is populated through multiple ionization pathways, as shown in Fig. 1(b). The ionization process carries valuable information about the coherent states. After the complete pump-probe pulse sequence, the photoelectron signal exhibits oscillations, commonly called “quantum beats,” as a function of the time delay τd between the pump and probe pulses. These quantum beats provide insights into the dynamics of the coherent wave packet and the energy spacing involved. By analyzing these quantum beats, we can extract information about the temporal evolution of the electron wave packet and the associated energy dynamics. In the following, we will illustrate this coherent process through the temporal evolution of the electron wave packet.

The electron wave packet for the atom initially in its ground state |1⟩ induced by the pump pulse can be written as $$\begin{gathered} |ψ1⟩=∑j⟨j|U1|1⟩|j⟩+∫dp⟨p|U1|1⟩|p⟩ \\ ≡∑jU1j1|j⟩+∫dpU1p1|p⟩, \end{gathered}$$(1)where |j⟩ denotes the jth bound state with energy εj. The state |p⟩ represents the continuum state with momentum p. U1 stands for the evolution operator in the presence of the pump pulse. Next, the electrons evolve freely over a duration, td [td=τd−(T1+T2)/2], which induces a stationary evolution phase to the state. This can be expressed as |ψ1(td)⟩=∑jU1j1e−iεjtd|j⟩+∫dpU1p1e−ip22td|p⟩,(2)where τd is sufficiently long to ensure that the second pulse does not overlap with the first one, as shown in Fig. 1(a). After the second weak and ultrashort probe pulse, the wave packet of electrons becomes $$\begin{gathered} |ψ2(td)⟩=∑b(∑jΛjb1e−iεjtd+∫dppΛb1e−ip22td)|b⟩ \\ +∫dp(∑jΛjp1e−iεjtd+∫dqqΛp1e−iq22td)|p⟩, \end{gathered}$$(3)where Λγαβ≡U2αγU1γβ, U2bj=⟨b|U2|j⟩, U2bp=⟨b|U2|p⟩, and |b⟩ represents the bth bare bound state. U2pj=⟨p|U2|j⟩, U2pq=⟨p|U2|q⟩, and |q⟩ represents the continuum state with momentum q. U2 represents the evolution operator perturbed by the second pulse. For numerical calculations, we can project Eq. (3) onto the scattering state of the atom to detect the photoelectron yield, which is given by $$\begin{gathered} P(p,td)=∑b∑jΛbp1Λjp1*e−i(εb−εj)td+|∫dqqΛp1e−iq22td|2 \\ +∑j[∫dqΛp1jΛqp1*e−i(εj−q22)td+c.c.]. \end{gathered}$$(4)

The time-resolved photoelectron signal exhibits quantum beats due to interfering pathways. Fourier analysis of the quantum beat signal P(p,td) can analyze the quantum beat pattern. However, in the realistic numerical calculations, the range of the time delay τd is finite [compared to the analytical derivation where the range is (−∞,∞)], which introduces abrupt changes in the photoelectron signals. To suppress this background noise, we apply a Gaussian envelope given by g(τd)=exp[−2ln2(τd−t′)2τF2], where τF represents the full-width at half-maximum of the laser pulse. Here, t′=t0+T/2, with t0 and T denoting the starting point and duration of the time delay, respectively. The resulting quantum beat spectrum is given by $$\begin{gathered} G(p,Ω)=F[P(p,td)·g(τd)], \\ =F[P(p,td)]*F[g(τd)]≡G˜(p,Ω)*g˜(Ω), \end{gathered}$$(5)where $$\begin{gathered} G˜(p,Ω)=∑b∑jΛbp1Λjp1*δ[Ω−(εj−εb)] \\ +∬dqdq′Λqp1Λq′p1*δ[Ω−(q′22−q22)] \\ +∑j{∫dqΛjp1Λqp1*δ[Ω−(q22−εj)] \\ +∫dqΛqp1Λjp1*δ[Ω−(−q22+εj)]}, \end{gathered}$$(6)and g˜(Ω)=12πae−iΩt′−Ω24a,a=2ln2τF2.(7)

We can observe that the information regarding the energy gap is encoded in the quantum beat spectrum. After performing the convolution, as shown in Eq. (5), we obtain $$\begin{gathered} G(p,Ω)=12πa×{∑b∑je−i[Ω−(εj−εb)]t′e−[Ω−(εj−εb)]24aΛbpgΛjpg* \\ +∬dqdq′e−i[Ω−(q′22−q22)]t′e−[Ω−(q′22−q22)]24aΛqpgΛq′pg* \\ +∑j∫dqe−i[Ω−(q22−εj)]t′e−[Ω−(q22−εj)]24aΛjpgΛqpg* \\ +∑j∫dqe−i[Ω−(−q22+εj)]t′e−[Ω−(−q22+εj)]24aΛqpgΛjpg*}. \end{gathered}$$(8)

The Gaussian envelope is used solely for noise reduction, which does not alter the underlying physical results. As a result, we can extract the interference information exclusively from G˜(p,Ω) as given by Eq. (6).

Quantum beats provide valuable information about impulsively induced coherence. In the next step, we aim to reconstruct the populations of bound states induced by the pump pulse by isolating the different interference processes present in the quantum beat spectrum. The population of the bth bound state, denoted as |U1b1|2, is related to the quantity Λbp1, which is defined as Λbp1≡U2pbU1b1. As a result, the population of the bound state |b⟩ can be expressed as |U1b1|2=|Λbp1|2|U2pb|2.(9)

In our previous scheme^[46], the probe pulse was designed to be extremely long so that electrons are completely ionized at the end of the pulse, i.e., ∫|U2pb|2dp≈1. In this case, once Λbp1 is known, the population of the bound state |b⟩ is given by |U1b1|2=∫|Λbp1|2dp.(10)

By referring to Eq. (6), |Λbp1| can be obtained from the first term in the quantum beat spectrum. However, obtaining the quantum beat spectrum in full-dimensional space is challenging in practical experiments. It requires detecting the three-dimensional photoelectron momentum distribution. Second, for every polarized angle θ, every azimuthal angle φ, and every momentum p, one needs to collect signals at different time delays and then perform a Fourier transform of the time delay to obtain the 3D quantum beat spectrum. The system exhibits cylindrical symmetry for a linearly polarized laser pulse, and a single arbitrary azimuthal angle φ is sufficient. Furthermore, it should be noted that the electrons need to be completely ionized.

We now propose an improvement to the previous reconstruction scheme to make it more feasible experimentally. In our scheme, we only detect the photoelectron momentum distribution along the polarization direction of the laser field, i.e., the z-axis. Instead of using an extremely long and strong probe pulse, we suggest using a weak and ultrashort attosecond pulse for the probing measurement, which only perturbs the system. In this case, the population of the bound state |b⟩ is given by |U1b1|2=|Λbp1|2|U2pb|2.(11)

First, we calculate |Λbp1| from the quantum beat spectrum. According to Eq. (6), the first term only contains information about the bound states characterized by a delta function. The finite range of τd in numerical calculations results in a non-negligible spacing dΩ. As a result, the characteristic signals appearing in the quantum beat spectrum will span a certain width along Ω. To obtain Λbp1, we integrate the beat frequency Ω near δεbj on both sides of Eq. (8), where δεbj=εb−εj, and b≠j (with b,j=1,2,3,…). This integration leads to the following expression, $$\begin{gathered} ∫δεbj−dεδεbj+dεG(p,Ω)dΩ \\ ≈12πa∫δεbj−dεδεbj+dεe−i(Ω−δεbj)t′e−(Ω−δεbj)24aΛbp1Λjp1*dΩ. \end{gathered}$$(12)

After further calculations, we obtain |∫δεbj−dεδεbj+dεG(p,Ω)dΩ|=|Λbp1||Λjp1|,(13)where 2δε represents the approximate distribution width, which observes that the artificial envelope g(τd) no longer appears in Eq. (13), confirming that the envelope does not impact the physical results. For each pair of distinct bound states, one equation involving Λbp1 can be obtained. When Cn2≥n (where n is the number of bound states), specifically for n≥3, it is possible to solve for Λbp1.

Second, we calculate the photoionization cross section |U2pb|. An analytical expression for the photoionization cross section of the ground state can be found in Refs. [50,51]. Experimental detection of the photoionization cross sections for excited states can be found in Ref. [52]. In our numerical simulations, we calculate |U2pb| using the time-dependent perturbation theory. In the weak and short probing measurement regime, the one-photon transition is the dominant process. The transition amplitude in the velocity gauge can be expressed as^[53,54]U2pb=−∫−T2/2T2/2⟨p|A(t)·∇|b⟩ei(E−εb)tdt,(14)where E=p2/2. This equation indicates that we can reconstruct the population of the bound state induced by the pump pulse by calculating Eq. (11).

We now examine this reconstruction scheme for realistic atomic hydrogen by solving the time-dependent Schrödinger equation (TDSE) in the velocity gauge and dipole approximation, i∂tψ(r,t)=[H0−iA(t)·∇]ψ(r,t),(15)where H0=−∇2/2−1/r is the field-free Hamiltonian. The vector potential is given by A(t)=A(t)z^, with A(t,τd)=A1f1(t)sin(ω1t)+A2f2(t−τd)sin[ω2(t−τd)],(16)where a sin-squared envelope is used, fi(t)=sin2(πtTi),(17)with Ti being the duration of the ith laser pulse. Ai=Ii/ωi represents the vector potential amplitude, where Ii is the peak intensity, and ωi is the photon energy for the ith pulse. The time interval τd corresponds to the interval between the envelope maxima of the first and second pulses, as depicted in Fig. 1(a).

To propagate the wave function in time, we employ the split-Lanczos propagator^[55–57]. In our calculations, the wave function is expressed as a unique sum of spherical harmonics, and the radial part of the wave function is discretized using the finite-element discrete variable representation (FE-DVR) method^[58–60]. To solve the TDSE, we split the solution into inner and outer regions at specific time points ti (with i=1,2,…,Nsplit) using an absorbing function s(r), ψ(r,t)=s(r)ψ(r,t)+[1−s(r)]ψ(r,t).(18)

The absorbing function s(r) is defined as follows: s(r)=1−11+e(r−Rc)/Δ,(19)where Rc represents the size of the inner box. The inner wave function follows the TDSE strictly, while the outer part is propagated using the Coulomb–Volkov propagator^[61,62].

Our numerical simulations set the maximal angular momentum as lmax=10. The radial box size is rmax=350 atomic units (a.u.), and the inner box size is Rc=200a.u. The time step used in the simulations is Δt=0.01a.u. The pump pulse used in our simulations has an intensity of 1015W/cm2 and a wavelength of 100 nm. The pulse duration, denoted as T1, is set to 2.5 optical cycles, which is approximately equal to 834.4 as. The pump pulse prepares a superposition of bound states and continuum. The photon energy of the first pulse, ω1, is blue detuned with respect to the |3p⟩ state and red detuned with respect to the |4p⟩ state. Under these parameters, the primary contribution to electron dynamics arises from the one-photon transition. After interacting with the pump pulse, electrons predominantly populate the bound states |1s⟩, |2p⟩, |3p⟩, and |4p⟩, as illustrated in Fig. 2. The probe pulse, on the other hand, is extremely short, lasting only one optical cycle (approximately 275.7 as). Like the pump pulse, the probe process is primarily governed by a one-photon transition.

After the interaction is completed, we obtain the final momentum distribution P(p,τd) along the polarized direction of the laser field for different time delays τd, as shown in Fig. 3(a). The photoelectron signals exhibit quantum beats and demonstrate the presence of attosecond coherence in both bound states and continuum. Interestingly, they display a distribution of light–dark cycles as a function of τd. In our calculations, we consider a range of time delays from τd(start)=58a.u. to τd(end)=3598a.u., with a spacing of 4 a.u., resulting in a total of 886 interval points. This number is significantly smaller than the 2048 time delay samples used in our previous scheme. The choice of a large time delay (τd≥58a.u.) ensures that the coherent state prepared by the first pulse remains unaffected by the probe pulse, indicating no overlap between the two laser pulse sequences. To extract the interference information encoded in the momentum distribution[Fig. 3(a)], we perform a Fourier transform of the quantum beat signal P(p,τd), yielding a momentum-frequency correlation spectrum G(p,Ω). To mitigate the effect of the sudden change in the photoelectron signals induced by the time delay, we apply a Gaussian window function g(τd) with a width equal to one-fourth of the range [τd(end)−τd(start)]. The resulting quantum beat pattern is given by G(p,Ω)=(2π)−1∫P(p,τd)g(τd)exp(−iΩτd)dτd.

Figure 3(b) shows the frequency-resolved photon-electronic spectroscopy, i.e., the quantum-beat spectrum G(p,Ω). By applying the Gaussian window, we effectively reduce the background noise compared to Fig. 1(c) in Ref. [46]. The spectrum exhibits two prominent interference structures, for which horizontal lines, labeled as L1,L2,…, are observed at a constant frequency Ω, and upward-opening parabolic curves are also observed. These characteristics observed in the full-dimensional system are consistent with those obtained in the 1D case^[46]. Notably, the characteristic signals in the spectrum have a certain width, which is wider than the 1D case, despite the use of the same spacing dΩ=0.00178a.u.

As can be seen from Eq. (6), we can identify the physical origins of the two prominent interference structures observed in Fig. 3(b). The horizontal lines in the spectrum correspond to the first term of Eq. (6), specifically when Ω=εj−εb in the δ function. This can be visualized as interference path 1, as depicted in Fig. 1(b). Initially, electrons populated in the ground state |1⟩ transition to different excited states |b⟩ and |j⟩ under the influence of the first laser pulse U1. During the subsequent free evolution time td, a dynamical quantum beat phase is accumulated, which is equal to the red shaded area in Fig. 1(b) and given by |εj−εb|td. With the second pulse U2, the excited electrons from states |j⟩ and |b⟩ are ultimately ionized into the same continuum state |p⟩, resulting in interference. This accumulated phase manifests as the horizontal lines in the quantum beat spectrum, indicating interference between electron wave packets originating from two different bound states.

Similarly, the parabolic curves in the spectrum result from the third term of Eq. (6), occurring when Ω=q22−εj. This corresponds to interference path 2 in Fig. 1(b). Electrons are ionized into the continuum state |q⟩ and excited to the state |j⟩ by the first pulse U1. During the subsequent free evolution, these electrons accumulate a relative phase given by the purple shaded area in Fig. 1(b) and equal to |q22−εj|td. After the second pulse U2, which transitions the electrons from |q⟩ and |j⟩ to the same |p⟩ continuum state, interference occurs. This interference between electron wave packets originating from the continuum and the bound state gives rise to the observed parabolic curves in the quantum beat spectrum.

It is important to note that these two interference structures also appear in the range Ω<0, but the opening of the parabolic curves is inverted, induced by the fourth term of Eq. (6). Additionally, it should be mentioned that due to the degeneracy of the angular momentum l, contributions from different l states within the same n manifold cannot be distinguished.

The horizontal lines correspond to the energy differences of the bound states of the hydrogen atom. In our reconstruction, we specifically focus on the bound states that electrons primarily populate at the end of the first pulse. That is, our goal is to reconstruct the populations of the first four bound states. By referring to the first term of Eq. (6) and the property of the δ-function, we can determine the following values: Ω(L1)≡δε21=ε2−ε1=0.375a.u., Ω(L2)≡δε31=ε3−ε1=0.444a.u., Ω(L3)≡δε32=ε3−ε2=0.0694a.u., and Ω(L4)≡δε42=ε4−ε2=0.0938a.u. Here, L1,L2,L3,L4 mark the different horizontal lines, as shown in Fig. 3(b). However, due to the finite range of the time delay τd in the numerical calculations, the resulting horizontal lines in the quantum beat spectrum have a certain width along Ω. Based on the treatment described in Section 2.2, we can derive a set of equations from Eq. (13), |Λ1p1||Λ2p1|=|∫0.375−dε0.375+dεG(p,Ω)dΩ|≡a1,(20)|Λ1p1||Λ3p1|=|∫0.4444−dε0.4444+dεG(p,Ω)dΩ|≡a2,(21)|Λ2p1||Λ3p1|=|∫0.0694−dε0.0694+dεG(p,Ω)dΩ|≡a3,(22)|Λ2p1||Λ4p1|=|∫0.0938−dε0.0938+dεG(p,Ω)dΩ|≡a4.(23)

The right-hand side of the equations is obtained numerically, where a1(2,3,4) corresponds to the integration of the horizontal lines L1(2,3,4). By solving these four equations, we can obtain the values of |Λ1p1|, |Λ2p1|, |Λ3p1|, and |Λ4p1|. Using time-dependent perturbation theory, we calculate the photoionization cross section [Eq. (14)] for the first four bound states of the hydrogen atom. Combining this with Eq. (11), we can retrieve the populations of the ground state |1⟩, the first excited state |2⟩, the second excited state |3⟩, and the third excited state |4⟩, which are given by |U111|2, |U121|2, |U131|2, and |U141|2, respectively.

During the reconstruction process, we obtain different values of p and find that a smaller value of p yields better results (p=0.04455a.u.). This originates from the fact that the signals are stronger at smaller momenta, which decreases the numerical calculation errors. To further validate the accuracy of our reconstruction method, we introduce changes to the wavelength and optical cycle of the first pulse while keeping the second pulse unchanged (I2=1014W/cm2, ω2=15eV, T2=1 optical cycle). In Fig. 4, we compare the reconstructed populations with the results of solving the TDSE. It can be observed that the reconstructed populations closely match the results from the TDSE calculations, demonstrating the effectiveness and accuracy of our reconstruction technique. Compared to the other reconstructed populations, more disagreements always exist for the reconstructed state |1⟩. In addition, the reconstruction error increases as the population of the ground state increases. This can be attributed to the depletion of the ground state, which makes the transition amplitude U2p1, calculated by the perturbation theory [Eq. (14)], deviate from the actual value slightly.

To investigate the effect of different dΩ values on the accuracy of the reconstructed populations, we consider a scenario where the first pulse has a wavelength of λ1=100 nm and a duration of T1=2.5 optical cycles. The results are summarized in Table 1. For smaller values of dΩ=0.00102a.u., the accuracy of the reconstruction is comparable to that of dΩ=0.00178a.u. When we increase dΩ to 0.00205 a.u., the reconstructed population of state |4⟩ deviates from the result of solving the TDSE. To understand the dependence of the reconstructed populations on dΩ, we analyze the quantum beat at a specific momentum p=0.04455a.u. in Fig. 5, where L1(2,3,4) corresponds to the horizontal lines shown in Fig. 3. According to Eqs. (20)–(23), the distributions at L1, L2, L3, and L4 should be separate and isolated peaks to ensure the integration range and accuracy of the reconstructed results. From Fig. 5(c), we can observe that when dΩ=0.00205a.u., the distribution at L3 becomes mixed with other peaks, significantly affecting the accuracy of the reconstructed results. Thus, it is crucial to carefully choose an appropriate value of dΩ to ensure well-separated peaks in the distributions of the quantum beat spectrum G(p,Ω) along the horizontal lines, improving the accuracy of the reconstructed populations.Table 1.

The Constructed Populations for Different Frequency Intervals dΩ (in a.u.)^a

In conclusion, we have developed a weak-field-induced quantum beat method that enables the determination of bound state populations, providing in-depth insights into the dynamics of excited states in strong-field phenomena. By applying this method to the strong-field induced ionization of atomic hydrogen, we have effectively validated the accuracy of the reconstructed populations by comparing them with results obtained from the time-dependent Schrödinger equation (TDSE). This approach exhibits the potential for further expansion in investigating complex systems extending beyond hydrogen atom ionization.

Our method focused exclusively on electrons that predominantly occupy the |1s⟩ and |np⟩ states through one-photon transitions induced by the strong pump pulse. To determine the transition amplitude caused by the probe pulse, we included the contributions from these dominant angular momentum states while disregarding those from others. Specifically, we approximated the reconstructed population in the |2p⟩ (|3p⟩, |4p⟩) state as the population of the |2⟩ (|3⟩, |4⟩) state. Moreover, the photoelectron quantum-beat spectrum is a valuable tool for identifying quantum interference pathways involving bound and continuum states, primarily driven by one-photon transitions. Considering the third term of Eq. (6), extracting information regarding transitions between these continuum states becomes feasible. For future work, we intend to refine our method further to address these limitations, enabling us to accurately interpret the results and understand the information obtained from the beat spectra.