Direct reconstruction of isolated XUV or soft x-ray attosecond pulses from high-harmonic generation streaking spectra

To simulate single-atom HHG streaking spectra under a linearly polarized weak XUV (or SXR) IAP and a strong IR (or MIR) laser, we can employ the previously established EQRS model [55]. This model is an improvement of SFA [56]. In the SFA, the single-atom induced dipole along the polarization direction at a fixed time delay τ between the two pulses can be written as

$$\begin{align} x(t,\tau) = i\int_{0}^{t}\mathrm{d}t^{^{\prime}}\left(\frac{\pi }{\varepsilon +i(t-t^{^{\prime}})}\right)^{3/2} \times d^{*}[\,p_{st}(t,t^{^{\prime}}) - A(t)]e^{-iS_{st}(t,t^{^{\prime}})}d[\,p_{st}(t,t^{^{\prime}})-A(t^{^{\prime}})] E(t^{^{\prime}})+ c.c., \end{align} \tag{ 1 }$$

where p_st is the canonical momentum, S_st is the quasi-classical action, $E(t) = E_{\text{IR}}(t) + E_{\text{XUV}}(t-\tau)$ is the combined electric field of IR and XUV pulses, and $A(t) = A_{\text{IR}}(t) + A_{\text{XUV}}(t-\tau)$ is the corresponding vector potential. Note that negative τ means the XUV pulse comes earlier than the peak of IR laser. Here equation (1) is expressed for IR and XUV pulses, but also can be used for MIR and SXR pulses.

Different terms in the time-dependent dipole $x(t,\tau)$ can be separated in the following:

$$\begin{equation}x(t,\tau) = x_1(t)+x_2(t,\tau)+x_3(t,\tau)+x_4(t,\tau). \end{equation} \tag{ 2 }$$

Each term in the right hand of the equation can be explicitly expressed as

$$\begin{align} x_1(t)& = i\int_{0}^{t}\mathrm{d}t^{^{\prime}}\left(\frac{\pi}{\varepsilon+i(t-t^{^{\prime}})}\right)^{3/2} f^{\,\,*}[\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t)]e^{-iS_{st}(t,t^{^{\prime}})}f\,[\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t^{^{\prime}})] \nonumber\\ & \quad \cdot [\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t)][\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t^{^{\prime}})] E_{\text{IR}}(t^{^{\prime}})+c.c., \end{align} \tag{ 3 }$$

$$\begin{align} x_2(t,\tau) = &-i\int_{0}^{t}\mathrm{d}t^{^{\prime}}\left(\frac{\pi}{\varepsilon+i(t-t^{^{\prime}})}\right)^{3/2} f^{\,\,*}[\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t)]e^{-iS_{st}(t,t^{^{\prime}})}f\,[\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t^{^{\prime}})] \nonumber \\ &\cdot [\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t)]A_{\text{XUV}}( t^{^{\prime}}-\tau) E_{\text{IR}}(t^{^{\prime}})+c.c., \end{align} \tag{ 4 }$$

$$\begin{align} x_3(t,\tau) = &-i\int_{0}^{t}\mathrm{d}t^{^{\prime}}\left(\frac{\pi}{\varepsilon+i(t-t^{^{\prime}})}\right)^{3/2} f^{\,\,*}[\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t)]e^{-iS_{st}(t,t^{^{\prime}})}f\,[\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t^{^{\prime}})] \nonumber \\ &\cdot A_{\text{XUV}}(t - \tau)[\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t^{^{\prime}})] E_{\text{IR}}(t^{^{\prime}})+c.c., \end{align} \tag{ 5 }$$

and

$$\begin{align} x_4(t,\tau)& = i\int_{0}^{t}\mathrm{d}t^{^{\prime}}\left(\frac{\pi}{\varepsilon+i(t-t^{^{\prime}})}\right)^{3/2}f^{\,\,*}[\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t)]e^{-iS_{st}(t,t^{^{\prime}})}f\,[\,p_{st}(t,t^{^{\prime}})-A_{\text{IR}}(t^{^{\prime}})] \nonumber\\ & \quad\cdot A_{\text{XUV}}(t-\tau)A_{\text{XUV}}(t^{^{\prime}}-\tau)E_{\text{IR}}(t^{^{\prime}})+c.c., \end{align} \tag{ 6 }$$

where f(p) is related to the dipole transition matrix element d(p) from the ground state to a continuum state as $f(p) = d(p)/p$ with electron momentum p and the stationary action $S_{st}(t,t^{^{\prime}})$ without including the effect of the XUV pulse is

$$\begin{align} S_{st}(t,t^{^{\prime}}) = I_{p}(t-t^{^{\prime}})-\frac{1}{2}p_{st}(t,t^{^{\prime}})(t-t^{^{\prime}})+\frac{1}{2}\int_{t^{^{\prime}}}^{t}\mathrm{d}t^{^{\prime\prime}}A_{IR}^{2}(t^{^{\prime\prime}}). \end{align} \tag{ 7 }$$

Note that the $x_1(t)$ term gives the induced dipole solely by the IR laser, the $x_4(t, \tau)$ term is by the XUV pulse alone, and the $x_2(t,\tau)$ and $x_3(t,\tau)$ terms give the coupling effect between the two pulses.

The first three terms on the right side of equation (2) in frequency domain can be factorized as

$$\begin{equation} x_{1}^{\text{SFA}}(\omega,I) = N^{\text{SFA}}(I)^{1/2}W_{1}(\omega)d^{\text{SFA}}(\omega), \end{equation} \tag{ 8 }$$

$$\begin{equation} x_{2}^{\text{SFA}}(\omega,I,\tau) = N^{\text{SFA}}(I)^{1/8} W_{2}(\omega,\tau)|d^{\text{SFA}}(\omega)|^{3/4}, \end{equation} \tag{ 9 }$$

and

$$\begin{equation} x_{3}^{\text{SFA}}(\omega,I,\tau) = N^{\text{SFA}}(I)^{1/8} W_{3}(\omega,\tau)|d^{\text{SFA}}(\omega)|^{2/5}. \end{equation} \tag{ 10 }$$

where $W_{1}(\omega)$ is the wave packet of the return electron driven solely by the IR laser. $W_{2,3}(\omega,\tau)$ are complex quantities similar to the electron wave packet in $x_{1}(\omega,I)$ and only depend on the external field. Here the peak intensity I of the dressing IR laser is explicitly expressed in the equations. With the same idea of quantitative rescattering (QRS) model [57, 58], ionization probability $N^{\,\text {SFA}}(I)$ and transition dipole moment $d^{\text{SFA}}(\omega)$ under the SFA are replaced by the 'exact' values of $N^{\text {QRS}}(I)$ and $d^{\text {QRS}}(\omega)$ in the QRS and EQRS models. The ionization probability $N^{\,\text {QRS}}(I)$ can be calculated by using the PPT model [59], and $d^{\text {QRS}}(\omega)$ is computed by using the 'exact' wave function for the bound and continuum states within the single-active electron (SAE) approximation. Thus, in the EQRS model, the single-atom induced dipole by the XUV and IR pulses at a time-delay τ can be computed from

$$\begin{align} x^{\text {EQRS}}(\omega,I,\tau) = x_1^{\text {QRS}}(\omega,I)+ x_2^{\text {EQRS}}(\omega,I,\tau)+ x_3^{\text {EQRS}}(\omega,I,\tau). \end{align} \tag{ 11 }$$

In the above equation, the $x_4(\omega,\tau)$ term is neglected because of its much smaller amplitude. Details of the derivation of EQRS model can be found in [55].

Note that HHG spectra obtained from equation (11) modulate as a function of τ. We have checked that such modulation appears also from TDSE-based harmonic spectra, and thus should be also from the experiment. Detail about the TDSE simulations can be seen in section 3.5.

When retrieving the spectral phase of XUV pulse from the HHG streaking spectra, we assume that the IR laser pulse and the spectral distribution of the XUV pulse are known. We do not need to know the exact intensity of the XUV pulse, and it is assumed weak, but strong enough to modulate the HHG spectra. The XUV pulse in the frequency domain without the time delay can be expressed as

$$\begin{equation} E_{\text{XUV}}(\omega) = U_{\text{XUV}}(\omega)e^{i[\phi_{\text{XUV}}(\omega)]}. \end{equation} \tag{ 12 }$$

Here $U_{\text{XUV}}(\omega)$ is the spectral amplitude, and $\phi_{\text{XUV}}(\omega)$ is the spectral phase to be retrieved. We assume that the XUV pulse is composed of multiple sub-pulses in the frequency domain. Thus, it can be expanded as in the following:

$$\begin{equation} E_{\text{XUV}}(\omega) = \sum_{i = 0}^{\infty}E_{i}(\omega), \end{equation} \tag{ 13 }$$

in which

$$\begin{equation} E_{i}(\omega) = u_{i}(\omega)e^{i[\phi_{\text{XUV}}(\omega)]}, \end{equation} \tag{ 14 }$$

$$\begin{equation} u_{i}(\omega) = \varepsilon_{i}e^{-2\ln 2 \frac{(\omega-\omega_{i})^{2}}{D_{i}^2}}, \end{equation} \tag{ 15 }$$

where $u_{i}(\omega)$ is a Gaussian envelope centered at ω_i with the full-width-at-half-maximum (FWHM) value of D_i and the amplitude of ε_i . In equation (14), the XUV spectral phase $\phi_{\text{XUV}}(\omega)$ is maintained, and it can be expanded into terms of a Taylor series about the central frequency ω_i . For the i-th component, it can be written as

$$\begin{equation} \phi_{\text{XUV}}(\omega) = \Phi_{i} + \alpha_{i}(\omega-\omega_{i}). \end{equation} \tag{ 16 }$$

Here the constant term $\Phi_{i}$ is the carrier-envelope phase (CEP), and the linear spectral phase term $\alpha_{i} = \frac{\mathrm{d} \phi_{\text{XUV}}(\omega)}{\mathrm{d}\omega}|_{\omega_{i}}$ is a constant group delay (GD) time.

By taking the Fourier transform of the i-th XUV field $E_{i}(\omega)$ in the frequency domain, it becomes

$$\begin{equation} E_{i}(t) = \frac{\varepsilon_{i}}{2\pi}\int_{-\infty}^{+\infty}e^{-2\ln2\frac{(\omega-\omega_i)^2}{D_{i}^{2}}}e^{i[\Phi_i+\alpha_i(\omega-\omega_i)]}e^{-i\omega t}\mathrm{d}\omega. \end{equation} \tag{ 17 }$$

And its temporal expression can be written as

$$\begin{equation} E_{i}(t) = \frac{\varepsilon_{i}D_{i}}{2\sqrt{2\pi\ln2 }}e^{-\frac{D_{i}^{2}(t-\alpha_{i})^2}{8\ln2}}\cos(\omega_{i}t-\Phi_{i}). \end{equation} \tag{ 18 }$$

The derivation of equation (18) is presented in appendix 'Derivation of the XUV temporal expression in equation (18)'. At a time delay τ, the i-th temporal XUV field can be written as

$$\begin{equation} E_{i}(t,\tau) = \frac{\varepsilon_{i}D_{i}}{2\sqrt{2\pi\ln2 }}e^{-\frac{D_{i}^{2}(t+\tau-\alpha_{i})^2}{8\ln2}}\cos[\omega_{i}(t+\tau)-\Phi_{i}]. \end{equation} \tag{ 19 }$$

Note that temporal fields in equations (18) and (19) are expanded using Gaussian envelope basis functions, the temporal form of IAP pulse can be arbitrary.

We have checked that for a spectral region overlapped with the XUV pulse, $|x_3^{\text {EQRS}}(\omega,I,\tau)|^2$ is about one order of magnitude smaller than $|x_2^{\text {EQRS}}(\omega,I,\tau)|^2$, thus the modulated HHG spectra are mainly caused by the interference of x₁ and x₂ terms. We have checked that this conclusion does not depend on the parameters of the IR laser, the XUV pulse, nor on the atomic target. In equation (4), we assume that all terms involving the IR laser in the integral change slowly over the integration interval in comparison with the vector potential of $A_{\text{XUV}}(t^{^{\prime}},\tau)$ such that the IR electric field can be pulled out from the integral. Thus, under the adiabatic approximation, equation (4) is simplified to

$$\begin{equation} x_{2}(t,\tau) \approx -iC_{\text{IR}}E_{\text{IR}}(t)\int_{0}^{t}A_{\text{XUV}}(t^{^{\prime}},\tau)\mathrm{d}t^{^{\prime}}+c.c., \end{equation} \tag{ 20 }$$

where C_IR is a constant related to the IR laser. With the integral in the following

$$\begin{equation} \int_{0}^{t}-A_{\text{XUV}}(t^{^{\prime}},\tau)\mathrm{d}t^{^{\prime}}\propto E_{\text{XUV}}(t,\tau), \end{equation} \tag{ 21 }$$

we can obtain

$$\begin{equation} x_{2}(t,\tau) \propto E_{\text{XUV}}(t,\tau) E_{\text{IR}}(t). \end{equation} \tag{ 22 }$$

The wave packet $W_{2}(\omega,\tau)$ in equation (9) is the same in the SFA and EQRS [55], thus equation (22) is valid in the EQRS model. Multiplying by $e^{i\omega_{j}t}$ and then integrating both sides of equation (22) with time t, we get

$$\begin{equation} \int_{-\infty}^{+\infty}x_{2}(t,\tau)e^{i\omega_{j}t}\mathrm{d}t \propto \int_{-\infty}^{+\infty}\sum_{i = 0}^{\infty}E_{i}(t,\tau)E_{\text{IR}}(t)e^{i\omega_{j}t}\mathrm{d}t. \end{equation} \tag{ 23 }$$

Only when $\omega_i = \omega_j$, the right hand of above equation has a non-zero value. We thus get the following relation:

$$\begin{equation} x_{2}(\omega_{i},\tau) \propto \int_{-\infty}^{+\infty}E_{i}(t,\tau)E_{\text{IR}}(t)e^{i\omega_{i}t}\mathrm{d}t. \end{equation} \tag{ 24 }$$

The derivation of equation (24) can be found in appendix 'Derivation of equation (24)'. By substituting equation (18) into equation (24), and adopting the adiabatic approximation, one has

$$\begin{equation} x_{2}(\omega_{i},\tau) \propto \frac{\varepsilon_{i}}{2}e^{i(\Phi_{i}-\omega_{i}\tau)}E_{\text{IR}}(\alpha_{i}-\tau). \end{equation} \tag{ 25 }$$

The derivation of equation (25) is given in appendix 'Derivation of equation (25)'. This formula can be further written as

$$\begin{equation} |x_{2}(\omega_i,\tau)|^{2} \propto \varepsilon_{i}^{2}E_{\text{IR}}^{2}(\alpha_{i}-\tau). \end{equation} \tag{ 26 }$$

We assume that the CEP of IR laser pulse is zero. The variation of the CEP only leads to a linear term difference in the reconstructed phase, which does not alter the temporal shape and the duration of IAP pulse. For a given ω_i , the maximum of $|x_{2}(\omega_i,\tau)|^{2}$ is reached by varying τ when $\alpha_i - \tau$ is zero such that the $E_{\text{IR}}^{2}(\alpha_i - \tau)$ reaches its maximum, and we label the satisfied delay time as τ_i . Thus, we have

$$\begin{equation} \tau_i = \alpha_{i} = \frac{\mathrm{d}\phi_{\text{XUV}}(\omega)}{\mathrm{d}\omega}|_{\omega_i}. \end{equation} \tag{ 27 }$$

If the $\Delta \omega = \omega_{i+1} -\omega_i$ is small enough, we replace the discrete ω_i and τ_i by variables ω and $\tau(\omega)$ by omitting the subscript i. Equation (27) is rewritten as

$$\begin{equation} \tau(\omega) = \frac{\mathrm{d}\phi_{\text{XUV}}(\omega)}{\mathrm{d}\omega}. \end{equation} \tag{ 28 }$$

Thus the spectra phase of XUV pulse can be obtained in the following:

$$\begin{equation} \phi_{\text{XUV}}(\omega) = \int_{\omega_s}^{\omega}\tau(\omega)\mathrm{d}\omega + \phi_s. \end{equation} \tag{ 29 }$$

For convenience, $\phi_{\text{XUV}}(\omega_s)$ at the boundary frequency is set as an arbitrary constant φ_s so that $\phi_{\text{XUV}}(\omega)$ at the central frequency of XUV pulse is equal to zero.

Note that the spectral phase of XUV pulse can be retrieved from the time-delayed spectra of $|x_{2}(\omega,\tau)|^2$. The HHG streaking spectra in the interested spectral region are mainly resulting from the interference of x₁ and x₂ terms. The former one does not depend on τ, thus for a given ω, the maximum of $|x_2|^2$ along τ is in coincidence with the maximum in the envelope of HHG streaking spectra. This forms the basis for retrieving the spectral phase of XUV pulse from the time-delayed HHG streaking spectra.

We first check the validity of equation (26) through two numerical examples. In the simulations, the IR pulse has a central wavelength of 800 nm, a FWHM duration of 5 fs, and a CEP of zero. The IR peak intensity of 2.5 × 10¹⁴ W cm⁻² provides a sufficient spectral range to cover the XUV pulse with a central photon energy of 71.3 eV and a spectral width of 9 eV. We choose two different phases of XUV pulse as shown in figure 1(a) together with their spectral intensities (assuming a Gaussian distribution and the same for both cases). They correspond to the transform-limited (TL) and chirped XUV pulses, respectively. The peak intensity of the two XUV pulses is 2.5 × 10¹¹ W cm⁻², and their temporal profiles are plotted in figure 1(b).

Zoom In Zoom Out Reset image size

We show the spectra of $|x_{2}(\omega,\tau)|^2$ calculated by the EQRS model in figures 1(c) and (d). The dotted (black) closed curves are the contour lines indicating the time delays where the spectral peak intensity decrease to its 40% for each given ω. The dashed (white) lines label the peak positions of spectra in terms of $\tau(\omega)$ in equation (28). One can see that the shapes of spectra are greatly modified by the spectral phase of the XUV pulse. We choose three photon energies: 67.5 (blue), 71.5 (black), and 75.5 eV (red), and plot the spectra of $|x_{2}(\omega,\tau)|^2$ as a function of the time delay in figures 1(e) and (f). We also show the profiles of $E_{\text{IR}}^{2}(\alpha_i - \tau)$ (dashed lines) in these figures, where α_i is obtained by equation (27). It is clearly shown that the peak of $|x_{2}(\omega,\tau)|^2$ always agrees with that of the intensity envelope of the shifted IR laser by varying the photon energy and the XUV phase. Thus, equation (26) is verified. Next, we show the modulated HHG spectra resulting from the interference of x₁ and x₂ terms in figures 1(g) and (h). They have similar envelope shapes in comparison with those of $|x_{2}(\omega,\tau)|^2$, and their peaks appear at the same positions as $E_{\text{IR}}^{2}(\alpha_i - \tau)$. Therefore, the relation of the GD time α_i and the time delay τ of peak intensity can be set through the HHG streaking spectra. This fact will be used to retrieve the spectral phase of XUV pulse. Note that in the retrieval the complete HHG streaking spectra are used, not from the $|x_{2}(\omega,\tau)|^2$ term only.

We next demonstrate how to retrieve the spectral phase (or temporal profile) of XUV isolated pulse with a narrow bandwidth from the HHG streaking spectra. Such spectra are first simulated by using the EQRS model with known parameters. The parameters of IR laser and XUV pulses are the same as those in figure 1 except for a variety of XUV phase. We assume that the spectral phase of the XUV pulse can be written as

$$\begin{equation} \phi_{\text{XUV}}(\omega) = \phi(\omega_0) + \alpha(\omega-\omega_0) + \frac{\beta}{2}(\omega-\omega_0)^2 + \frac{\gamma}{6}(\omega-\omega_0)^3, \end{equation} \tag{ 30 }$$

where ω₀ is the central frequency, $\phi(\omega_0)$ is a constant, α describes the group delay, β is the coefficient of the group delay dispersion (GDD), and γ accounts for the third-order dispersion (TOD). By choosing different values of these parameters, we can construct five different forms of the spectral phase as shown in figure 2(f). When we only keep the term $\phi(\omega_0)$, the spectral phase is a constant (solid line), leading to the transform-limited (TL) XUV pulse. If only the third term in equation (30) is used, the spectral phase is labeled as 'GDD' (dotted lines). The spectral phase is noted as 'TOD' (dashed-dotted lines) by only using the forth term in equation (30).

Zoom In Zoom Out Reset image size

Figures 2(a)–(e) show the simulated HHG streaking spectra of $|x^{\text {EQRS}}(\omega,\tau)|^2$ from equation (11) by using different XUV phases. The spectra are plotted in a photon-energy region overlapped with the spectral region of XUV pulse only. In figure 2(a), the range of time delay is set at about three half-optical cycles (1.33 fs) of the IR laser. One can clearly see the half-optical cycle modulation structure. This structure also appears if other phases are applied (not shown). To see clearly the change of the modulated spectra, we set the range of time delay as about one half-optical cycle in figures 2(b)–(e). The dashed (white) lines are the contour lines in which the intensity is 40% of the peak intensity for each photon energy. With the help of these lines, we can visually identify the change of the streaking spectra by varying the spectral phase of the XUV pulse, especially the change of time delay for the peak intensity. Since the HHG streaking spectra are sensitive to the XUV phase, it would be possible to precisely retrieve the XUV phase from such spectra.

Staring from HHG streaking spectra in figure 2, we reconstruct the spectral phase of XUV pulse as in the following. We take figure 2(b) for example. First, for a specific frequency ω_i , i.e. photon energy of 71.3 eV, the original data describing the variation of spectral intensity with the time delay is shown in figure 3(a). Here the step in time delay is assumed as 10 a.s. To obtain the accurate peak position in each oscillation cycle, we increase the number points in time delay and get the spectral intensity by the Newton interpolation. The results are shown in dashed (red) line in figure 3(a). We then use these peak positions to construct a smooth envelope by using the least square fitting method as shown by a dotted (blue) line in figure 3(a). The peak of this envelope determines the value of time delay when the modulated spectra reach the maximum at a given photon energy. This forms one point labeled as 'cross' in the solid (black) line in figure 3(b). This procedure can be repeated for other photon energies to fully obtain the data of $\tau(\omega)$ (solid black line). Then we smooth these data (dashed red line) to get rid of some peculiar points due to errors in the fitting procedure for the smooth envelope of HHG modulation intensity. Finally, according to equation (29), the spectral phase can be reconstructed from HHG streaking spectra.

Zoom In Zoom Out Reset image size

For the three selected HHG streaking spectra shown in figures 2(a), (b), and (d), the reconstructed XUV spectral phases are shown in the left column of figure 4. Clearly the retrieved phases are very accurate compared to the 'input' data, some discrepancies appear near the wings where the spectral amplitudes are already small. To further check the accuracy of the retrieval algorithm, we show the reconstructed temporal profiles of XUV pulses on the right column of figure 4. The spectral amplitudes of XUV pulse are assumed known as Gaussian. One can see that the retrieved XUV pulses overlap almost perfectly with the 'input' ones. The error in the FWHM duration of XUV pulse is only a few attoseconds, and the biggest relative error is less than 3%.

Zoom In Zoom Out Reset image size

We check whether our retrieval algorithm can be employed to reconstruct the spectral phase of IAPs with a broad bandwidth in the SXR region. We choose three attosecond pulses with the central photon energy of 300 eV and the spectral FWHM bandwidth of 45 eV, which would result in a 40 as pulse under the transform-limited condition. The peak intensity of the SXR pulse is fixed at 5.0 × 10¹¹ W cm⁻², and three different input spectral phases are shown as dashed (red) lines in figures 5(d)–(f). The distribution of spectral intensity is also shown as a dotted (blue) line in figure 5(d). We employ a MIR laser pulse with pulse duration of 20 fs, wavelength of 1600 nm, peak intensity of 5.0 × 10¹⁴ W cm⁻², and CEP of zero. The simulated HHG streaking spectra by the MIR laser pulse and the time-delayed isolated SXR pulses are plotted in figures 5(a)–(c). We cannot visually observe the interference fringes in these spectra because the optical period of the SXR pulse is much smaller than the time-delay range. One can see that the prominent modulation structures are dramatically changed by varying the SXR spectral phase. Using the same procedure as discussed in section 3.2, we reconstruct the spectral phases shown as solid (black) lines in figures 5(d)–(f). And the retrieved temporal profiles of isolated SXR pulses are plotted in figures 5(g)–(i). The retrieved pulse durations are 480.1 as, 70.8 as, and 41.4 as, compared with the 'input' ones of 487.4 as, 70.3 as, and 41.7 as, respectively. The biggest relative error in the duration of SXR pulse is less than 1.5%.

Zoom In Zoom Out Reset image size

In the examples above, the assumed form of the phase over a broad spectral region is relatively simple. Can the spectral phase be accurately retrieved if it has a more complicated behavior in a more broader SXR region. We consider a chirped SXR pulse with the central photon energy of 280 eV and the spectral bandwidth of 90 eV. Its spectral phase is plotted as a dashed (red) curve in figure 6(b). We employ the same MIR laser pulse, and keep the same peak intensity of SXR pulse as those in figure 5. With the simulated HHG streaking spectra in figure 6(a), we obtain the retrieved spectral phase of SXR pulse shown as a solid (black) line in figure 6(b). It shows some derivations with the 'input' phase. The reconstructed temporal profile of SXR pulse, as shown in figure 6(c), presents a good overlap with the 'input' one. The retrieved pulse duration is 27.7 as while the 'input' value is 24.3 a.s. This example verifies that our method is capable of accurately retrieving ultrashort isolated SXR pulses with complicated phase.

Zoom In Zoom Out Reset image size

In the above numerical investigations, we assume that the IR (or MIR) laser pulse can produce the high harmonic spectrum, whose spectral region can well overlap with that of the XUV (or SXR) pulse. However, the IR (or MIR) laser pulse cannot be known exactly in the experiment. We next investigate how the change of IR (or MIR) laser can influence the accuracy of our retrieval method. We take an example of the isolated SXR pulse in figures 5(e) and (h), and the parameters of the MIR laser are kept the same except for the peak intensity. For three peak intensities, the simulated HHG spectra by the MIR laser pulse alone are varied in both the spectral structure and spectral intensity, see figure 7(d). As shown in figures 7(a)–(c), with the increase of the MIR intensity, the modulated HHG spectra maintain their shapes, but become much brighter compared to the background signals. This does not affect the accuracy of the retrieved spectral phase. As shown in figure 7(e), three retrieved spectral phases at different MIR intensities are all in good agreement with the 'input' value. The agreement of retrieved temporal profiles of SXR pulses with the 'input' one in figure 7(f) further shows that the accuracy of retrieval method does not change by varying the peak intensity of the MIR laser. We have also checked the CEP effect of IR (or MIR) laser, which contributes to a linear term in the spectral phase, i.e. α_i differs from τ_i in equation (27) by a constant, leading to an overall time shift of XUV (or SXR) pulse only. Thus the accuracy of retrieval method is not changed by changing the CEP of the IR (or MIR) laser.

Zoom In Zoom Out Reset image size

We next check how the intensity of the IAP pulse affects the accuracy of our retrieval method. For the examples of the isolated SXR pulse shown in figures 5(d) and (g), all parameters are kept the same except for the peak intensity of the IAP. We show the simulated HHG spectra in figures 8(a)–(c) by varying the intensity of the SXR pulse by an order of magnitude. The retrieved spectral phases and the temporal pulse profiles are shown in figures 8(d) and (e), respectively. One can see the relative error in the retrieved duration of the SXR pulse is less than 4%. Note that with the decrease of intensity of the IAP, the modulation structure in the HHG streaking spectra becomes faint which might not be desirable in such experiments.

Zoom In Zoom Out Reset image size

We also check another factor of the time-delay step for obtaining the 'input' HHG streaking spectra. In the experiment, this is varied due to the capability of precisely controlling the two-color laser beams. We use the same SXR pulse and MIR laser as those in figure 7, and the peak intensity of MIR laser is set at 5.0 × 10¹⁴ W cm⁻². The pulse duration of the MIR laser is changed to 10 fs. The central photon energy of the SXR pulse is 300 eV, and its spectral bandwidth is 45 eV. We set the step of time delay between the SXR pulse and the MIR laser as 20 as, 30 as, 40 as, and 50 as, respectively. These values are all larger than the modulation period of HHG spectra along the time delay, which is close to one optical cycle of the central SXR frequency, i.e. 15 a.s. Thus the resulted streaking spectra in figures 9(a)–(d) show different patterns visually. The retrieved spectral phases of SXR pulse are plotted in figure 9(e). We have checked that the result obtained by the step of 20 as is almost identical with the 'input' one, and it can be considered as a reference here. One can see that the retrieved phase changes a little bit randomly by increasing the step of time delay. The reconstructed temporal profiles of SXR pulses are plotted in figure 9(f). The peak of retrieved pulse is slightly shifted with the time-delay step, but its shape remains and the retrieved duration still has the high accuracy. From the duration values labeled in the figure, the biggest relative error in the SXR pulse duration is about 4%. So the accuracy of retrieval maintains if the step of time delay in measuring the HHG streaking spectra is much bigger than the spectra modulation period. In a real experiment, the time delay cannot be perfectly stabilized and a jitter with the scale of few tens attoseconds should be considered. We have checked that retrieved results are accurate and stable even a random and reasonable shift in time delay is added for simulating the HHG streaking spectra.

Zoom In Zoom Out Reset image size

The HHG streaking spectra used for retrieving the spectral phase of XUV (or SXR) pulse so far are all simulated by using the EQRS model. Here we test if the spectral phase can be accurately reconstructed from the 'experimental' spectra, which is taken to be the HHG streaking spectra obtained by numerically solving the TDSE [57]. The TDSE solution is a three dimensional simulation obtained under the dipole and the single-active electron approximations with the effective potential taken from [60]. Based on the TDSE simulations, we define that $D_{\text{IR,XUV}}(t,\tau)$ is the dipole acceleration under the IR laser and the time-delayed XUV pulse, $D_{\text{XUV}}(t,\tau)$ and $D_{\text{IR}}(t)$ are the dipole accelerations solely caused by the XUV pulse at a time delay τ and by the IR pulse, respectively. We then use these quantities in the frequency domain, and let

$$\begin{align} |D(\omega,\tau)|^{2} = |D_{\text{IR,XUV}}(\omega,\tau)|^2-|D_{\text{IR}}(\omega)+D_{\text{XUV}}(\omega,\tau+nT_{\text{IR}})|^{2}, \end{align} \tag{ 31 }$$

where n is an integer and it meets the condition of $nT_{\text{IR}}\geqslant \frac{1}{2}T_{\text{total}}$ with the IR period T_IR, and the IR total pulse duration T_total. $D_{\text{IR,XUV}}(\omega,\tau)$ gives the dipole response when the IR laser and the XUV pulse are overlapped while $D_{\text{IR}}(\omega)+D_{\text{XUV}}(\omega,\tau+nT_{\text{IR}})$ is the response when the two pulses are well separated in time. Both the spectra from the two terms can be directly measured in the experiment at different time delays between the IR laser and the XUV pulse. Thus, through equation (31), modulation of the HHG streaking spectra with time delay becomes clearly visible and is consistent with that simulated by the EQRS model. We take the example of retrieving a narrow-bandwidth XUV pulse. The parameters of IR laser and XUV pulse are the same as those in section 3.2. For two different spectral phases, the TDSE-based modulation spectra are modified, see figures 10(a) and (b). Using our phase retrieval method, the spectral phases can be well reconstructed, see figures 10(c) and (d). The retrieved temporal profiles of the XUV pulse also show good agreement with the 'input' ones, see figures 10(e) and (f). The retrieved pulse durations of 215.9 as and 445.9 as can be well compared with the 'input' ones of 200.5 as and 486.3 as, respectively. The biggest relative error in the retrieved pulse duration is about 8% in these two cases.

Zoom In Zoom Out Reset image size

By taking the Fourier transform of the i-th electric field $E_{i}(\omega)$ in the frequency domain in equation (14), its temporal expression can be written as

$$\begin{align} E_{i}(t)& = \frac{1}{2\pi}\int_{-\infty}^{+\infty}\varepsilon_{i}e^{\frac{-2\ln2(\omega-\omega_i)^2}{D_{i}^{2}}}e^{i[\Phi_i+\alpha_i(\omega-\omega_i)]}e^{-i\omega t}\mathrm{d}\omega\nonumber\\ & = \frac{\varepsilon_{i}}{2\pi}e^{i(\Phi_i-\omega_{i}t)}\int_{-\infty}^{+\infty}e^{\frac{-2\ln2(\omega-\omega_i)^2}{D_{i}^{2}}}e^{i(\alpha_i-t)(\omega-\omega_i)}\mathrm{d}\omega. \end{align} \tag{ 32 }$$

We define

$$\begin{align} a = \frac{\sqrt{2\ln2}}{D_i};b = -i(\alpha_i-t);x = \omega-\omega_i, \end{align} \tag{ 33 }$$

and the field expression of $E_i(t)$ can be simplified as

$$\begin{align} E_i(t)& = \frac{\varepsilon_{i}}{2\pi}e^{i(\Phi_i-\omega_{i}t)}\int_{-\infty}^{+\infty}e^{-(a^{2}x^{2}+bx)}\mathrm{d}x. \end{align} \tag{ 34 }$$

The above equation can be calculated as

$$\begin{align} E_i(t)& = \frac{\varepsilon_{i}}{2\pi}e^{i(\Phi_i-\omega_{i}t)}\int_{-\infty}^{+\infty}e^{-(ax+\frac{b}{2a})^{2}+\frac{b^2}{4a^2}}\mathrm{d}x = \frac{\varepsilon_{i}}{2\pi}e^{i(\Phi_i-\omega_{i}t)}e^{\frac{b^2}{4a^2}}\int_{-\infty}^{+\infty}e^{-(ax+\frac{b}{2a})^{2}}\mathrm{d}x\nonumber\\ & = \frac{\varepsilon_{i}}{2a\pi}e^{i(\Phi_i-\omega_{i}t)}e^{\frac{b^2}{4a^2}}\int_{-\infty}^{+\infty}e^{-(ax+\frac{b}{2a})^{2}}\mathrm{d}\left(ax+\frac{b}{2a}\right) = \frac{\varepsilon_{i}}{2a\sqrt{\pi}}e^{i(\Phi_i-\omega_{i}t)}e^{\frac{b^2}{4a^2}}. \end{align} \tag{ 35 }$$

The definite integral of $\int_{-\infty}^{+\infty}e^{-x^{2}}\mathrm{d}x = \sqrt{\pi}$ is used.

Therefore, the complex expression of $E_i(t)$ is

$$\begin{align} E_i(t) = \frac{\varepsilon_{i}D_{i}}{2\sqrt{2\pi\ln2}}e^{-\frac{D_{i}^{2}(t-\alpha_i)^{2}}{8\ln2}}e^{i(\Phi_i-\omega_{i}t)}. \end{align} \tag{ 36 }$$

If we only keep the real part, $E_i(t)$ can be expressed as

$$\begin{align} E_i(t) = \frac{\varepsilon_{i}D_{i}}{2\sqrt{2\pi\ln2}}e^{-\frac{D_{i}^{2}(t-\alpha_i)^{2}}{8\ln2}}\cos(\omega_{i}t-\Phi_i). \end{align} \tag{ 37 }$$

The equation (23) is expressed as

$$\begin{equation} \int_{-\infty}^{+\infty}x_{2}(t,\tau)e^{i\omega_{j}t}\mathrm{d}t \propto \int_{-\infty}^{+\infty}\sum_{i = 0}^{\infty}E_{i}(t,\tau)E_{\text{IR}}(t)e^{i\omega_{j}t}\mathrm{d}t. \end{equation} \tag{ 38 }$$

Its right hand side can be written as

$$\begin{equation} \int_{-\infty}^{+\infty}\sum_{i = 0}^{\infty}E_{i}(t,\tau)E_{\text{IR}}(t)e^{i\omega_{j}t}\mathrm{d}t = \sum_{i = 0}^{\infty} F_{ij}(\tau), \end{equation} \tag{ 39 }$$

where

$$\begin{equation} F_{ij}(\tau) = \int_{-\infty}^{+\infty}E_{i}(t,\tau)E_{\text{IR}}(t)e^{i\omega_{j}t}\mathrm{d}t. \end{equation} \tag{ 40 }$$

By substituting the explicit form of $E_{i}(t,\tau)$ in equation (18) into the above equation, one obtains

$$\begin{align} F_{ij}(\tau) = & \frac{\varepsilon_{i}D_{i}}{2\sqrt{2\pi\ln2}}\int_{-\infty}^{+\infty}e^{-\frac{D_{i}^{2}(t+\tau-\alpha_i)^{2}}{8\ln2}}E_{\text{IR}}(t) \cos[\omega_{i}(t+\tau)-\Phi_i]e^{i\omega_{j}t}\mathrm{d}t. \end{align} \tag{ 41 }$$

We define

$$\begin{align} C_{i} = \frac{\varepsilon_{i}D_{i}}{2\sqrt{2\pi\ln2}}; B_{i} = \frac{D_{i}^{2}}{8\ln2}; t^{^{\prime}} = t + \tau - \alpha_i, \end{align} \tag{ 42 }$$

and then we have

$$\begin{align} F_{ij}(\tau)& = C_{i}e^{i\omega_{j}(\alpha_{i}-\tau)}\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau) \cos[\omega_{i}(t^{^{\prime}}+\alpha_{i})-\Phi_i]e^{i\omega_{j}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}\nonumber\\ & = \frac{C_{i}}{2}e^{i\omega_{j}(\alpha_{i}-\tau)}\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}e^{i(\omega_{i}t^{^{\prime}} + \omega_{i}\alpha_{i}-\Phi_{i})} E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)e^{i\omega_{j}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}\nonumber\\ & \quad +\frac{C_{i}}{2}e^{i\omega_{j}(\alpha_{i}-\tau)}\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}e^{-i(\omega_{i}t^{^{\prime}} + \omega_{i}\alpha_{i} - \Phi_{i})} E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)e^{i\omega_{j}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}\nonumber\\ & = \frac{C_{i}}{2}e^{i[(\omega_{j}+\omega_{i})\alpha_{i}-\Phi_{i}-\omega_{j}\tau]}\!\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}e^{i (\omega_{j}+\omega_{i})t^{^{\prime}}} E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)\mathrm{d}t^{^{\prime}}\nonumber\\ &\quad + \frac{C_{i}}{2}e^{i[(\omega_{j}-\omega_{i})\alpha_{i}+\Phi_{i}-\omega_{j}\tau]}\!\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}e^{i (\omega_{j}-\omega_{i})t^{^{\prime}}} E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)\mathrm{d}t^{^{\prime}}. \end{align} \tag{ 43 }$$

Then the common integral term in the above equation is

$$\begin{align} \int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)e^{i\omega^{^{\prime}}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}& = E_{\text{IR}}(\alpha_{i}-\tau)\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}+i\omega^{^{\prime}}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}\nonumber\\ & = E_{\text{IR}}(\alpha_{i}-\tau)\sqrt{\frac{\pi}{B_{i}}}e^{-\frac{\omega^{^{\prime} 2}}{4B_{i}}}, \end{align} \tag{ 44 }$$

where $\omega^{^{\prime}}$ = $\omega_{j}+\omega_{i}$ or $\omega_{j}-\omega_{i}$. The integral with the same form as equation (44) also appears in equation (48), below which the approximated approach to get its value is given in equation (54). When $\omega^{^{\prime}} = \omega_{j}+\omega_{i}$, the exponential term in equation (44) is approximately equal to zero since $D_i \lt \lt \omega^{^{\prime}}$. Therefore, equation (43) can be simplified as

$$\begin{align} F_{ij}(\tau) = \frac{\varepsilon_{i}}{2}e^{i[(\omega_{j}-\omega_{i})\alpha_{i}+\Phi_{i}-\omega_{j}\tau]}E_{\text{IR}}(\alpha_{i}-\tau)e^{-2\ln2\frac{(\omega_{j}-\omega_{i})^2}{D_{i}^{2}}}. \end{align} \tag{ 45 }$$

In equation (45), only when $\omega_i = \omega_j$, $F_{ij} (\tau)$ has a considerable value in the summation in equation (39), thus equation (38) can be reduced as

$$\begin{equation} x_{2}(\omega_{i},\tau) \propto \int_{-\infty}^{+\infty}E_{i}(t,\tau)E_{\text{IR}}(t)e^{i\omega_{i}t}\mathrm{d}t. \end{equation} \tag{ 46 }$$

This is equation (24) in section 2.2.

Let

$$\begin{align} F_{i}(\tau) = &\int_{-\infty}^{+\infty}E_{i}(t,\tau)E_{\text{IR}}(t)e^{i\omega_{i}t}\mathrm{d}t = \frac{\varepsilon_{i}D_{i}}{2\sqrt{2\ln2\pi}}\int_{-\infty}^{+\infty}e^{-\frac{D_{i}^{2}(t+\tau-\alpha_i)^{2}}{8\ln2}}E_{\text{IR}}(t)\nonumber\\ &\times \cos[\omega_{i}(t+\tau)-\Phi_i]e^{i\omega_{i}t}\mathrm{d}t. \end{align} \tag{ 47 }$$

Using defined constants in equation (42), we have

$$\begin{align} F_{i}(\tau)& = C_{i}e^{i\omega_{i}(\alpha_{i}-\tau)}\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)\cos[\omega_{i}(t^{^{\prime}}+\alpha_{i})-\Phi_i]e^{i\omega_{i}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}\nonumber\\ & = \frac{C_{i}}{2}e^{i\omega_{i}(\alpha_{i}-\tau)}\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}e^{i(\omega_{i}t^{^{\prime}} + \omega_{i}\alpha_{i}-\Phi_{i})}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)e^{i\omega_{i}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}\nonumber\\ & \quad +\frac{C_{i}}{2}e^{i\omega_{i}(\alpha_{i}-\tau)}\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}e^{-i(\omega_{i}t^{^{\prime}} + \omega_{i}\alpha_{i} - \Phi_{i})}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)e^{i\omega_{i}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}\nonumber\\ & = \frac{C_{i}}{2}e^{i[\omega_{i}(2\alpha_{i}-\tau)-\Phi_{i}]}\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}e^{i (2\omega_{i}t^{^{\prime}})}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)\mathrm{d}t^{^{\prime}}\nonumber\\ &\quad +\frac{C_{i}}{2}e^{i(\Phi_{i}-\omega_{i}\tau)}\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)\mathrm{d}t^{^{\prime}}. \end{align} \tag{ 48 }$$

In the following, we calculate the two integrals in the above equation separately. First we expand the electric field $E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)$ by Taylor series at the time $t^{^{\prime}} = 0$ as

$$\begin{align} E_{\text{IR}}&(t^{^{\prime}}+\alpha_{i}-\tau) = E_{\text{IR}}(\alpha_i-\tau)+E^{^{\prime}}_{\text{IR}}(\alpha_i-\tau)t^{^{\prime}}+\frac{E^{^{\prime\prime}}_{\text{IR}}(\alpha_i-\tau)}{2}t^{^{\prime} 2} + \ldots+ \frac{E_{\text{IR}}^{(n)}(\alpha_i-\tau)}{n!}t^{^{\prime} n}. \end{align} \tag{ 49 }$$

The $E_{\text{IR}}^{(n)}(\alpha_i-\tau)$ is the nth derivative of the electric field at time $t^{^{\prime}} = 0$. We keep the first three terms in the expansion, and substitute the expansion form of $E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)$ into the second integral of equation (48) as

$$\begin{align} \int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)dt& = E_{\text{IR}}(\alpha_i-\tau)\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}\mathrm{d}t^{^{\prime}}\nonumber\\ & \quad +E^{^{\prime}}_{\text{IR}}(\alpha_i-\tau)\int_{-\infty}^{+\infty}t^{^{\prime}}e^{-B_{i} t^{^{\prime} 2}}\mathrm{d}t^{^{\prime}}+\frac{E^{^{\prime\prime}}_{\text{IR}}(\alpha_i-\tau)}{2}\int_{-\infty}^{+\infty}t^{^{\prime} 2}e^{-B_{i} t^{^{\prime} 2}}\mathrm{d}t^{^{\prime}}. \end{align} \tag{ 50 }$$

Using the definite integral of the Gamma function in the following

$$\begin{align} \int_{-\infty}^{+\infty}x^{n}e^{-x^2}\mathrm{d}x = \left\{\begin{matrix} &0 && \mathrm{n\ is\ an\ odd\ number}\\ &\frac{(n-1)!!}{\sqrt{2}^n}\sqrt{\pi} && \mathrm{n\ is\ an\ even\ number}, \end{matrix}\right. \end{align} \tag{ 51 }$$

the second integral becomes

$$\begin{align} \int_{-\infty}^{+\infty}&e^{-B_{i} t^{^{\prime} 2}}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)\mathrm{d}t^{^{\prime}} = E_{\text{IR}}(\alpha_i-\tau)\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}\mathrm{d}t^{^{\prime}}+\frac{E^{^{\prime\prime}}_{\text{IR}}(\alpha_i-\tau)\sqrt{\pi}}{4B_{i}^{3}}. \end{align} \tag{ 52 }$$

Since the FWHM duration of the i-th XUV pulse $E_i(t)$ is shorter than the optical cycle of the IR field, the IR electric field can be considered slowly varying, meaning that $E^{^{\prime\prime}}_{\text{IR}}(\alpha_i-\tau) \approx 0$. Thus we have

$$\begin{align} &\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)\mathrm{d}t^{^{\prime}} = E_{\text{IR}}(\alpha_i-\tau)\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}\mathrm{d}t^{^{\prime}}. \end{align} \tag{ 53 }$$

We use the same approach to calculate the first integral of equation (48) as

$$\begin{align} \int_{-\infty}^{+\infty} e^{-B_{i} t^{^{\prime} 2}}e^{i2\omega_{i}t^{^{\prime}}}E_{\text{IR}}(t^{^{\prime}}+\alpha_{i}-\tau)\mathrm{d}t^{^{\prime}} & = E_{\text{IR}}(\alpha_i-\tau)\int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}e^{i2\omega_{i}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}\nonumber\\ & \quad +\sum_{n = 1}^{\infty}\frac{E_{\text{IR}}^{(n)}(\alpha_i-\tau)}{n!}\int_{-\infty}^{+\infty}t^{^{\prime} n}e^{-B_{i} t^{^{\prime} 2}}e^{i2\omega_{i}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}. \end{align} \tag{ 54 }$$

We expand the function of $e^{i2\omega_{i}t^{^{\prime}}}$ as

$$\begin{align} e^{i2\omega_{i}t^{^{\prime}}} = \sum_{m = 0}^{\infty}\frac{(i2\omega_{i})^m}{m!}t^{^{\prime} m}. \end{align} \tag{ 55 }$$

Thus we have

$$\begin{align} \sum_{n}\frac{E_{\text{IR}}^{(n)}(\alpha_i-\tau)}{n!}\int_{-\infty}^{+\infty}t^{^{\prime} n}e^{-B_{i} t^{^{\prime} 2}}e^{i2\omega_{i}t^{^{\prime}}}\mathrm{d}t^{^{\prime}} & = \sum_{n = 1}^{\infty}\frac{E_{\text{IR}}^{(n)}(\alpha_i-\tau)}{n!}\int_{-\infty}^{+\infty}t^{^{\prime} n}e^{-B_{i} t^{^{\prime} 2}}\sum_{m = 0}^{\infty}\frac{(i2\omega_{i})^m}{m!}t^{^{\prime} m}\mathrm{d}t^{^{\prime}}\nonumber\\ & = \sum_{n,m}\frac{(i2\omega_{i})^m}{m!}\frac{E_{\text{IR}}^{(n)}(\alpha_i-\tau)}{n!}\int_{-\infty}^{+\infty}t^{^{\prime} n+m}e^{-B_{i} t^{^{\prime} 2}}\mathrm{d}t^{^{\prime}}\nonumber\\ & = \sum_{n,m}\frac{(i2\omega_{i})^m}{m!}\frac{E_{\text{IR}}^{(n)}(\alpha_i-\tau)}{B_{i}^{n+m+1}n!}\int_{-\infty}^{+\infty}t^{^{\prime\prime} n+m}e^{-t^{^{\prime\prime} 2}}\mathrm{d}t^{^{\prime\prime}}. \end{align} \tag{ 56 }$$

In the above equation, $t^{^{\prime\prime}} = \sqrt{B_{i}}t^{^{\prime}}$. By considering the integral of the Gamma function, we get

$$\begin{align} \sum_{n}\frac{E_{\text{IR}}^{(n)}(\alpha_i-\tau)}{n!}\int_{-\infty}^{+\infty}t^{^{\prime} n}e^{-B_{i} t^{^{\prime} 2}}e^{i2\omega_{i}t^{^{\prime}}}\mathrm{d}t^{^{\prime}} = \sum_{n,m}\frac{(i2\omega_{i})^m}{m!}\frac{E_{\text{IR}}^{(n)}(\alpha_i-\tau)}{B_{i}^{n+m+1}n!}\frac{\sqrt{\pi}(n+m-1)!!}{\sqrt{2}^{n+m}}, \end{align} \tag{ 57 }$$

where n + m is set as even number. Using the same argument for the second order derivative of $E_{\text{IR}}(\alpha_i-\tau)$ previously, $E^{(n)}_{\text{IR}}(\alpha_i-\tau) \approx 0$ with $n \geqslant 2$ (even number), so the above equation equals to zero. We define

$$\begin{align} a = \sqrt{B_i}; b = 2i\omega_i, \end{align} \tag{ 58 }$$

and we have the integral in the first term of equation (54) as

$$\begin{align} \int_{-\infty}^{+\infty}e^{-B_{i} t^{^{\prime} 2}}e^{i2\omega_{i}t^{^{\prime}}}\mathrm{d}t^{^{\prime}}& = \int_{-\infty}^{+\infty}e^{-a^2t^{^{\prime} 2}+bt^{^{\prime}}}\mathrm{d}t^{^{\prime}}\nonumber\\ & = e^{\frac{b^2}{4a^2}}\int_{-\infty}^{+\infty}e^{-(at^{^{\prime}}-\frac{b}{2a})^2}\mathrm{d}t^{^{\prime}} = \frac{\sqrt{\pi}}{a}e^{\frac{b^2}{4a^2}}. \end{align} \tag{ 59 }$$

Thus we have

$$\begin{align} \int_{-\infty}^{+\infty}&e^{-B_{i} t^{2}}e^{i2\omega_{i}t}\mathrm{d}t = \sqrt{\frac{\pi}{a}}e^{-\frac{8\ln2\omega_i^2}{D_i^2}}. \end{align} \tag{ 60 }$$

Since $\omega_i \gg D_i$, the above equation is approximately equal to zero. The first integral term in equation (48) can thus be ignored. Therefore, equation (48) can be written as

$$\begin{align} F_{i}(\tau) = &\frac{\varepsilon_{i}D_{i}}{\sqrt{32\pi\ln2}}e^{i(\Phi_{i}-\omega_{i}\tau)}E_{\text{IR}}(\alpha_{i}-\tau)\int_{-\infty}^{+\infty}e^{-\frac{D_{i}^{2}t^{2}}{8\ln2}}\mathrm{d}t\nonumber\\ = &\frac{\varepsilon_{i}}{2}e^{i(\Phi_{i}-\omega_{i}\tau)}E_{\text{IR}}(\alpha_{i}-\tau). \end{align} \tag{ 61 }$$

Thus equation (24) is expressed as

$$\begin{align} x_2(\omega_i,\tau) \propto \frac{\varepsilon_{i}}{2}e^{i(\Phi_{i}-\omega_{i}\tau)}E_{\text{IR}}(\alpha_{i}-\tau). \end{align} \tag{ 62 }$$

This is equation (25) in section 2.2.