Retrieval of photoionization group delay in chirped-attosecond pulses photoelectron streaking experiments

Photoionization time delays have been investigated in many streaking experiments in which an extreme ultraviolet (XUV) attosecond pulse is used to ionize the target in the presence of a dressing infrared laser field. The discrepancies between the photoionization time delays thus experiment measured and those from many sophisticated theoretical simulations have generated a great deal of controversy in recent years. The difficulty of achieving an accuracy of the retrieved time delays comes from two facts: a so-called wavepacket approximation is introduced to construct the photoionization electron wavepacket, this approximation is invalid if atto-chirp of XUV is non-zero; the other one is that the lower sensitivity of the streaking spectra to the phase of the photoionization transition dipole. Here we present a time delay retrieval method born from our recently proposed ‘phase retrieval of broadband pulses auto correlation’ (PROBP-AC) technology to overcome above limitations. We carefully exam the validity of our method and make a few compare with some other common used retrieval codes, and the simulations demonstrate that more accurate results can be retrieved using PROBP AC. Based on the present method, the angular dependent photoionization time delays can also be retrieved. Our investigation casts doubts on the measured group time delays in previous streaking experiments. We also mention here that a single photoionization group time delay at the XUV peak energy is not enough to represent a complete photoemission process; instead, a fully characterization of the photoionization group time delay over the whole bandwidth of the wave packet is required.


Introduction
The advent of attosecond technology since 2001 [1][2][3][4][5][6][7][8][9][10][11][12][13] paves the way to study the nature of electron dynamics in atomic [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30], molecules [31][32][33][34][35] and condensed materials [36][37][38][39][40][41]. In a typical attosecond streaking experiment, two-dimension photoelectron spectra are generated by extreme ultraviolet (XUV) attosecond pulses in the presence of a moderate intense infrared (IR) laser, with tunable time delay between the two pulses. Such spectra are called a spectrogram or a spectral trace. Traditionally, these streaking spectrogram are widely used to characterize the time domain attosecond pulses, by employing the frequency-resolved optical gating for complete reconstruction of attosecond bursts (FROG-CRAB) method [42] for the characterization of isolated attosecond pulses, and the 'reconstruction of attosecond beating by interference of two-photon transition' [43] method for the characterization of attosecond pulse trains, or the other published methods, such as PROBP [44], PROBO-AC [45], Volkoff transform generalized projections algorithm [46] and phase retrieval by omega oscillation filtering [47]. The spectrogram also contains the information of the so-called photoelectron wave packet, which consists of the XUV field and the complex photoionization transition dipole moment (TDM) of the target. Thus, the streaking trace also has been used to characterize the TDM phase. In fact, what have been interested usually is not the phase of the dipole, but rather the first derivative of the phase with respect to the electron energy. Such a derivative has been widely called photoemission time delay or 'photon-ionization time delay' in the attosecond science. The concept of photon-ionization time delay was an extension of the Eisenbud-Wigner-Smith delay in scattering process [48][49][50], and was generally reported in past decade. We should mention here that, in physics, the group delay is defined as the derivative of the spectral phase. Similarly, the photoionization time delay used here is the derivative of the photoionization TDM phase respect to energy so that it should be called a Wigner group delay (to avoid confusion, we still call it time delay or group time delay through out this paper). The topic of photonionization time delay has generated a great deal of controversy since the first experiment, by Schultze et al [23] in 2010, where a time delay of 21 as was reported between the ionization from the 2p and that from the 2s sub shells of Ne. After this experiment reported, many sophisticated theoretical calculations have all been devoted to reproduce this quantity, but few (or at least, only part of them) can get the 'right' number. Therefore, accurately extracting the information of photon-ionization time delay from these attosecond streaking experiments is still a widely open question and rather challenging test.
In order to characterize the photon-ionization time delay between two separate ionization channels which can be in the same target or in different two targets, usually one would extract two electron wave packets from corresponding two photoelectron spectrograms from a typical streaking experiment. The TDM phase difference between the two channels is expected to be the phase difference of the two electron wave packets, whereas the unknown XUV phase is assumed to be canceled after taking difference, this is the so-called wavepacket approximation (WPA), which was first derived in [51]. However, there will be a prominent error in the photon-ionization time delay if the atto-chirp of XUV pulse is large [52][53][54][55], one of the reason is that WPA is not valid when a large chirped phase XUV used to generate streaking trace. The relationship between TDM and XUV phases in the photoionization wavepacket is more complex than WPA (detail can be found in figure 3 in [52]), therefore, the unknown XUV phase cannot be directly removed by such a simple operation. Another reason is that the chirped XUV phase is much larger than that of TDM, leading to the sensitivity from XUV phase is much larger than that from TDM phase, this will make the TDM phase or photon-ionization time delay retrieval more difficult.
To address the above questions, in this work, we focus on accurately reconstructing the photoionization time delay (difference) of different ionization channels from photoelectron spectra generated by large chirped attosecond pulses. Our model is rely on the strong field approximation (SFA), which is a simple theory and calculation cost is small. Instead of employing the whole set of the streaking spectra as in the previous phase retrieval algorithms, we calculate the autocorrelation (AC) of the streaking spectrogram [45,56], and use it to extract the photoionization time delay. AC pattern has been applied to retrieve the broadband single attosecond pulses in two recent experiments in [45]. It was found that applying AC patterns to retrieve the attosecond pulses can achieve a much faster convergence by at least one or two orders of magnitude, and a more accurate spectral phases can be retrieved. In this paper, we show that the AC of the photoelectron spectra also reveals more clear variations with TDM phases. We still call the method in this paper PROBP AC for simplicity. By separating TDM and XUV from the photo-ionization wavepacket, we remove the WPA in the retrieval optimization calculation, and thus the XUV phase, TDM phase, as well as the IR laser field can all be accurately retrieved, even the atto-chirp of XUV pulse is large. With the help of this method, the whole photoionization time delay over all the electron-emit angles can also be determined.
The article is organized as the following: section 2 briefly discusses the principle of the theory; in section 3, we show how the AC patterns evolves with the TDM phases; in section 4 based on SFA model, we run various simulations using different targets and XUV pulses, and compare the retrieved photoionization group delays from method present here with that from PROBP and FROG-CRAB; section 5 shows the result for the results retrieved from TDSE input AC patterns in both high and low energy domain; in section 6, we extend our method to retrieve the angular dependent photoionization time delay. Section 7 discusses the error effect; section 8 gives the conclusion. Atomic units are used unless otherwise stated.

Model and theory
The measured photoelectron spectrogram can be modeled by SFA [42,57]: Here the polarization of the XUV, the IR, and the photoelectrons are all taken along the +z direction. p is the asymptotic momentum of the photoelectron, then the energy of the electron E = p 2 /2. E XUV (t) is the XUV field in the time domain. In the energy domaiñ where U(Ω) and Φ(Ω) are the XUV spectral amplitude and phase, respectively. In equation (1), is the vector potential of the IR field. t d is the time delay between XUV and IR. A positive t d means the XUV comes after the peak of the IR field. I p is the ionization potential energy of the target. The phase function ϕ(p, t) reads Equation (1) includes a single-photon TDM element d(E) = Eẑ|z|i , where |i is the initial bound state, and |Eẑ is the continuum state in which the electron has an asymptotic momentum p = √ 2E toward the +z direction. The bound and continuum states are eigenstates of the field-free Hamiltonian of the target atom. The transition dipole is a complex quantity d(E) = |d(E)|e iη(E) which can be calculated within the single-active-electron approximation. The details of evaluating the transition dipole can be found in reference [54]. The photoionization group time delay is then defined as A single-color, multi-cycle IR field with a frequency ω L and carrier-envelope phase φ CEP in the time domain is expressed as From the spectrogram S(E, t d ) an AC is defined by the following: The integration is performed typically cover the whole photoelectron spectrum.
In this work, we assume that the frequency ω L , carrier-envelope phase φ CEP , XUV spectral amplitude U(Ω) and TDM amplitude d(Ω) are known. Then the unknown functions are IR envelop f(t), XUV spectral phase Φ(Ω) and TDM phase η(Ω). In addition, there is a drawback of our model is that it cannot determine the absolute time t, that is to say, the output of our algorithm could be d(t − t 0 ), f(t − t 0 ) and E XUV (t − t 0 ), here t 0 is arbitrary. Equivalently, such uncertainty would add a linear term Ωt 0 to the TDM phase η(Ω) and t 0 to the photoionization time delay. Thus if one applies our method to the α and β spectrograms individually, the same temporal axis for the two extracted TDM cannot be guaranteed, as a result, the obtained time delay is uncertain. Therefore, the electron spectrogram have to be expressed as the total 'sum' of the traces from the α and β channels: Then the photoionization time delay between α and β channels is: where η α(β) is the TDM phase from channel α(β). Thus, the uncertainty t 0 to the photoionization time delay can be canceled. As a result, we get totally four unknown functions needed to be retrieved: f(t), Φ(Ω), η α and η β . In our method, the unknown functions are constructed by B-spline basis functions [58]:  [44].
In our method, the optimal parameters {a i , b i , c i , d i } are found by genetic algorithm (GA). From guessed parameters one can reconstruct η α , η β , Φ(Ω) and f(t) and then reconstruct AC Q r (t d1 , t d2 ) using equation (6) for the two ionization channels. Compared with the input AC Q i (t d1 , t d2 ), the goal of the algorithm is to minimize the merit function Here we discretize S(E, t d ) at two sets of grid points E(i), (i = 1, . . . , NE) and t d (j), (j = 1, . . . , Nt d ), respectively. Typically we chose Nt d = 100 to 500 points in time delay and NE = 500 to 600 points in electron energy. The GA runs a large number of generations (typically 2000 to 10 000 generations) until convergence is achieved. The efficiency of the GA depends sensitively to the choice of B-spline function parameters, the order k and the number of B-spline functions n. Experimenting on the choice of the input parameters for the GA is necessary. Our strategy is first to run over different combinations of (nj, k) for a few generations (about 300 in this work) and pick the basis set with the fastest change in the fitness function. Based on our observations, running k from 4 to 7 and nj from 5 to 7 would often guide us to select the best nj and k for our cases. In the optimization, once we get the best B-spline parameters for IR in the first retrieval calculation, we use these sets of parameters and do not scan them in the remaining retrieval calculations to save time. Details of GA and B-spline functions can be found in our previous publications [44,45,53,56].
Different with the other photon-ionization time delay retrieval technology (details of these methods can be found in [42,52]) so far we know, WPA is removed in our method thus the TDM phase and XUV phase are identified.

Streaking spectra and AC pattern
For simplicity, consider an XUV pulse with a Gaussian spectral amplitude and a quadratic spectral phase Φ(Ω) = a 2 (Ω−Ω 0 ) 2 (ΔΩ/2) 2 . The coefficient a 2 is a measure of the attochirp, or equivalently the duration of the chirped pulse compared with the TL duration [56]. Here we use Ω 0 = 70 eV as the XUV central energy and the full width at half maximum (FWHM) bandwidth ΔΩ = 11.5 eV, which is corresponding a TL pulse of about 160 as FWHM (if a 2 = 0). The targets are chosen Ar 3p and Ne 2p channels, the TDMs are calculated from solving schrödinger equation, details can be found in [54]. Figures 1(a), (c), (e) and (g) show the spectrograms and AC patterns for Ar 3p channel (figures 1(a) and (e) are from the TL pulse, figures 1(c) and (g) are from chirped pulse). Meanwhile, figures 1(b), (d), (f) and (h) are the spectrograms and AC patterns for Ne 2p channel (figures 1(b) and (f) are from the TL pulse, figures 1(d) and (h) are from chirped pulse). The spectrograms are simulated using SFA according to equation (1). The IR field is 800 nm in wavelength, cosine-squared envelope, 6.2 fs in FWHM duration, and 10 12 W cm −2 in peak intensity. The AC patterns are calculated based on equation (6) with t d1 , t d2 over one optical period. In this example, while one can see that the streaking spectra display hardly fine differences (except an energy shift from the different ionization energies) as the target changed, the change is much more pronounced in the AC pattern, even in the chirped XUV cases. The sensitivity of AC patterns respect to TDM phases makes our method easier to retrieve TDM phases. Also, we can easily distinguish the chirped pulse from the TL pulse from both the trace and AC patterns.

Test with spectrogram generated from SFA
The aim of this section is to test the validity of our method and compare the retrieved photoionization time delay with the most widely used technology, FROG-CRAB and with our previous published method PROBP. The input traces and corresponding AC patterns are calculated using SFA equation (1). Since both the input and output data are simulated by SFA, the error generated in a specific method is only from the approximations and model made by this method itself. The input TDMs are calculated using Tong's model potential, which can be found in reference [60].

Group delay between photoionization from 2p and 2s states of Ne
We first focus on testing the validity of our method by retrieving of photoionization time delay (difference) Δτ 2p-2s between 2p and 2s channels of Ne. In this example, two XUV pulses are used in the simulation, they differ only by their spectral phases. One is TL with a FWHM of 160 as and the other is chirped with a FWHM of 280 as. Both pulses have the same U(Ω) which takes a Gaussian form with a central photon energy Ω 0 = 105 eV and a FWHM bandwidth ΔΩ = 11.5 eV. The IR field is 800 nm in wavelength, cosine-squared envelope, 6.2 fs in FWHM duration, 10 12 W cm −2 in peak intensity, and φ CEP = 0. Numerically the input spectrograms S(E, t d ) are formed by 501 × 301 matrices with energy step dE = 0.2 eV and delay step dt d = 53 as. Figures 2(a) and (b) show the input dipole amplitude |d(Ω)| and phase η(Ω) for the two ionization channels, respectively. The dipole phase of 2p and 2s channels vary smoothly in the photon energy range of the XUV pulse. The comparison between the input and retrieved group delay Δτ 2p-2s are shown in figure 1(c) for the TL case and in figure 1(d) for the chirped case. We have retrieved the photoionization time delay using PROBP and FROG-CRAB presented in reference [52], and now we add the result from the present method to the comparison. The input photoionization time delay (difference) calculated from Tong's model potential is almost constant (4-5 as). We see clearly these quantities can be accurately retrieved by FROG-CRAB, PROBP and PROBP AC if the XUV pulse is TL. However, in the case of 280 as chirped XUV, the retrieved photoionization time delay (difference) from both PROBP and FROG-CRAB vary significantly with the photon energy (−3 to 12 as for PROBP and −10 to 20 as for FROG-CRAB), and the result from PROBP AC is still stable and matches well with the input data. Because FROG-CRAB and PROBP actually retrieve the photon-ionization wave packet instead of the TDM, and the trace are approximated as: where

then the photoionization time delay (difference) is
this is the so-called WPA. However, the energy domain wave packet does not always follow equation (13), this approximation is not valid when the atto-chirp of XUV is large, we have illustrated this in figure 3 in [52]. This problem happens in all FROG-CRAB type codes where the wave packets are constructed by equation (13). Since a certain amount of residual attochirp will always be present in a real experiment, the capability of PROBP AC to retrieve photoionization group delays with chirped XUV pulses is needed. Besides, as the TDM phase changes smoothly in this spectral range, the streaking trace thus become insensitive to the TDM phase, which lead to difficulties in extracting the TDM phase for FROG CRAB and PROBP, as shown in figure 1. We note that the streaking trace in figure 2(A) of Schultze et al [23] appears to be generated from a chirped XUV pulse since the trace is more like that in figures 1(c) and (d). Thus the retrieved 21 as time delay might include error due to the chirp of their XUV pulse. In general, one expected XUV pulse obtained from harmonic generation or free-electron-laser is chirped, and we do not know how large this atto-chirp is, such an uncertain might easily lead to substantial error in the characterization of photoionization time delay or TDM phase. In addition, in order to compare the rate of convergence, the computation times from the three methods are given: for a typical photoionizaiton time delay retrieval like in figure 2, PROBP AC method would take about 500-1000 iterations to reach the convergence with calculation time about 5-10 min, while for PROBP method would take about 2500-5000 iterations to reach the convergence with calculation time about 30-60 min. FROG CRAB is the slowest, the calculation would take a whole day to get the results in figure 2. All the calculations were done on a single Skylake-SP Intel Xeon CPU core. Clearly the PROBP-AC method proposed here is much faster. The performance of PROBP AC in retrieving chirped XUV pulse and IR laser field is presented in figure 3. As shown in this figure, we have successfully retrieved the spectral phases and temporal intensity for chirped XUV pulse, and the IR laser field is also accurately retrieved. A note is made here that when people measure XUV phase by other algorithms, the TDM phase is approximated from solving time independent Schrödinger equation and treated as a known function, this approximation helps the program to get the converge faster, however, it may also introduce an extra error to XUV phase if the TDM phase is not accurately calculated. Our method gives a way to cancel this error.

Group delay between photoionization from Ar(3p) and Ne(2p)
Then, we retrieve the photoionization group delay (difference) Δτ Ar-Ne between Ar(3p) and Ne(2p) with PROBP AC, and we compare these results with those from PROBP and FROG-CRAB. The input photoionization TDM amplitude and phase for Ar and Ne atoms are plotted in figures 4(a) and (b), respectively. The input spectrograms and AC patterns are computed using equation (1). Both the XUV pulses used here have a central photon energy Ω 0 = 60 eV with the FWHM bandwidth ΔΩ = 23 eV,  supporting a TL pulse with FWHM of 80 as. The IR field is 800 nm in wavelength, cosine-squared envelope, 8.8 fs in FWHM duration, and 10 12 W cm −2 in peak intensity. The input and retrieved group delays Δτ Ar-Ne are compared in figure 4(c) for the case of TL XUV pulse, and in figure 4(d) for the case of chirped XUV pulse with FWHM about 280 as. Details of the FROG-CRAB method can be found in [54]. The retrieved photo-ionization time delays from PROBP AC are more accurate than those from PROBP and FROG-CRAB, especially in the chirped case in figure 4(d). The error of FROG-CRAB comes from the well-known central momentum approximation, the remaining error of PROBP and FROG-CRAB come from the WPA, which have already been discussed detailed in section 4.1 in [52]. We can conclude that by getting rid of the WPA, PROBP AC is more stable than PROBP and FROG-CRAB.
We also note that the input Δτ Ar-Ne is from −100 to 20 as in the photon energy range 40 to 80 eV, this large photoionization time delay change indicates that a single photoionization group time delay at the XUV peak energy is not enough to represent a complete photoemission process; instead, a fully

Test with spectrogram generated from TDSE
In this section, we generate the input streaking traces and corresponding AC patterns by solving TDSE as the 'experiment input data' for both Ar and Ne targets. As we know, SFA is the only theory available today for retrieving the spectral phase and photoionization time delay, but this model is not accurate enough for describing the streaking spectra when the photoelectron energy is less than 30 or 40 eV since SFA does not take into account the interaction between the continuum electron and the ionic core. Here, we address this problem and try to see how much error would SFA introduce both in high photoelectron energy and low photoelectron energy domains. We use PCTDSE code [59] numerically to solve the three-dimension (3D) TDSE with the one-electron model potential given in [60]. The box size and number of grid points are chosen to ensure convergence.

High energy
Here, we try to retrieve photoionization time delay (difference) Δτ 2p-2s between 2p and 2s channels of Ne in case of TL and chirped XUV pulses. The two XUV pulses and IR laser field used here to generate the streaking traces are totally the same with those in section 4.1. The retrieved Δτ 2p-2s are plotted in figure 5(a) for TL XUV and (b) for chirped XUV to compare with the input ones, respectively. Even using the TDSE data as the input, the error of Δτ 2p-2s retrieval is quite small for both PROBP AC and PROBP in TL case shown in figure 5(a). This shows that both PROBP AC and PROBP can accurately retrieve Δτ 2p-2s when using TL XUV pulse with central photon energy of about 105 eV for the Ne target investigated here. Meanwhile, in the chirped case shown in figure 5(b), the retrieved Δτ 2p-2s from PROBP AC are 4 to 5 as, which still agree very well with the input value. The retrieved Δτ 2p-2s from PROBP are from −5 to 25 as, note that these numbers from TDSE input data are quite similar with that in SFA case (−7.5 to 20 as in section 4.1). This result indicates that the photo-ionization spectrograms calculated by SFA model is accurate enough in the high photoelectron energy region, where the Coulomb effect on the photoelectron can be neglected.

Low energy
Yet in this example, we use two XUV pulses with a same Gaussian form amplitude U(Ω) with a central energy Ω = 30 eV, the spectral FWHM is set to ΔΩ = 11.5 eV which supports an 160 as TL FWHM. The two XUV pulses only different in the atto-chirp, one is TL and the other has a larger chirp supporting the FWHM of about 300 as. The IR field is 800 nm in wavelength, cosine-squared envelope, 8.8 fs in FWHM duration, and 10 12 W cm −2 in peak intensity. From this new TDSE-simulated input IR-dressed spectrogram, we repeat our method to retrieve the photoionization group delay (difference) Δτ Ar-Ne between Ar(3p) and Ne(2p) with PROBP AC and PROBP. The results are shown in figure 6(a) for TL XUV case and (b) for chirped case, respectively. In this example, the input Δτ Ar-Ne is from −104 to 5 as, however, all the retrieved results prominently deviate from this number. This is not surprising since it has been well established that the SFA model does not accurately describe low-energy photoelectrons (16 to 34 eV in this example). Thus one has to conclude that it is not possible to accurately retrieve the partial wave phase shifts, or equivalently, the phase of low energy photoelectrons with SFA model.

Angular group delay retrieval
In most of the previous photoionization time delay experiments and theories the amplitude and phase of the photoelectron wave packet is studied along photoelectron emission angle θ k = 0 direction only, corresponding that the polarization of the XUV, the IR, as well as the photoelectrons are all taken along the same direction. By measuring the streaking trace in this direction, one can reconstruct the photoionization time delay or the TDM at θ k = 0. Meanwhile, the photoionization also generates electrons at other angles. In principle one would like to determine the photoionization time delay (difference) over all the angles. In this section, we apply our method in retrieving photoionization time delays (difference) Δτ Ar-Ne between Ar(3p) and Ne(2p) in arbitrary photoelectron emission angle θ k . It should be mentioned here that, although our example is based on the calculations of Ar(3p) and Ne(2p) target, the method is directly applicable to other targets or ionization channels. Specifically, we only simulate a large atto-chirped XUV case since the TL case have already been investigated in our previous work [53]. All the input data in this section are generated by TDSE. The XUV pulse has a central energy Ω = 60 eV, the spectral FWHM is set to ΔΩ = 11.5 eV, the temporal FWHM of this XUV is about 700 as, and IR laser field used to generated the input streaking trace is totally the same with that in section 4.2.
For simplicity, we directly give the final complex angular dependent TDM that we are interested: Where R nl (r) is the radial wave function, and Y lm is the spherical harmonic, η L (E) = − Lπ 2 + σ L (E) + δ L (E) is the phase shift for the partial wave with angular quantum number L, σ L = arg Γ L + 1 − iZ c /k is the Coulomb phase shift, and δ L is the phase shift due to the short-range potential. Details of the derivation can be found in our previous paper [53].
We take the Ar(3p) atom as example where n i = 3, l i = 1. Denote the radial matrix elements as W L (E) = R EL |r|R 31 . By working out the angular part we have the angular TDM: W 0 (E)e iη 0 (E) + W 2 (E)e iη 2 (E) 3 cos 2 θ k − 1 (17) for the case where the initial state has m i = 0. For the initial state with m i = ±1, we have From the above equation, once we obtain W 0 , W 2 , η 0 and η 2 the TDM amplitude and phase in any direction can be calculated. This can be done in two steps. First, we retrieve W 0 and G(E) = W 2 (E)cos[η 0 (E) − η 2 (E)] from IR-free photoionization differential cross section: and the β parameter here we only give the final expression of equations (19) and (20), the derivation can be found in our previous paper [53]. From equations (19) and (20) we have The left-hand side of equation (21) is a known function, whereas the four unknown functions on the right-hand side can be combined into two, W 0 (E) and G(E) = W 2 (E)cos[η 0 (E) − η 2 (E)]. We expand the two unknowns by B-spline functions: The unknown coefficients {a i , b j } are optimized by GA in order to obtain the best fit to the left-hand side of equation (21). From the retrieved W 0 (E), one can calculate W 2 (E) from equation (19). Then from the retrieved G(E) one can deduce the phase shift difference η 0 (E) − η 2 (E). Figures 7(a) and (b) are the input σ total and β parameter for Ar atoms. Mathematically, equation (21) is a simple quadratic functions of the unknowns, therefore it only takes a few seconds for our GA solver to converge. After about 10 000 generations we obtained W 0 (E) and G(E), then from equation (19), we also obtain the phase difference η 0 − η 2 . Figures 7(c) and (d) compare the input and retrieved radial matrix elements W 0 and W 2 . The input and retrieved data match very well. The retrieved phase difference η 0 − η 2 also agrees well with the input data. Similar calculation can be repeated for Ne(2p) channel, the retrieved W 0 and W 2 of Ne(2p) channel are shown in figures 9(a) and (b), which are consistent with the input one.
Then in the second step we retrieve the partial wave phase shift η 0 of Ar and Ne from AC pattern for any measured direction θ k : where the angular dependent trace is: θ in equation (25) is the angle between vector potential and the electron momentum. Note that although our goal is to retrieve the photoinoization time delay in the whole angular range, we actually do not need to measure the streaking trace at each of the angles, only one specific angle is enough for the retrieval. The four unknown functions, η 0 of Ar and Ne, as well as XUV phase and IR envelop are all expanded by the B-spline functions with fitting parameters, then we do the optimization calculation with GA until the code converges. Once this phase η 0 of Ar and Ne are retrieved, using the known W 0 , W 2 , and η 0 − η 2 , all the four parameters W 0 , W 2 , η 0 and η 2 as functions of E are obtained. Figures 8(a) and (b) and 9(c) and (d) are the input and output partial wave phase shifts η 0 and η 2 for Ar(3p) (figures 8(a) and (b)) and Ne(2p) (figures 9(c) and (d)) atoms. We can see from the figures that the The phase difference between the two partial waves vs photon energy from the input and from the retrieval also agree well (not shown).

Figure 8.
Test of the accuracy of the retrieved partial wave phase shift (a) η 0 and (b) η 2 for Ar 3p channel. The streaking trace and corresponding AC pattern are generated using TDSE with the known input parameters. The two phases retrieved from the AC pattern are compared to the input phases. The center energy of the XUV pulse is at 60 eV. Only chirped XUV pulse is used. phase shifts η 0 and η 2 can be accurately extracted even the XUV has a large atto-chirp. With the retrieved W 0 , W 2 , η 0 and η 2 we can construct the angular dependent TDM phase by equations (17) and (18 figure 10 we compare the input (a) and output (b) two-dimension angular dependent photoionization time delays (difference) Δτ Ar-Ne (E, θ k ). As expected, the retrieved Δτ Ar-Ne (E, θ k ) is in good agreement with the input one. Thus we can conclude that the chirped XUV photoionization time delays (difference) can be retrieved accurately if the central photon energy is about 60 eV and higher.
We should note here that the Ar 3p photoionziation data used here in figures 7(a) and (b) are very different from the relevant experimental data in [61], the difference here may come from the invalidity of  the single-electron approximation, which can not be removed in our model (or in any methods based on single-electron approximation). As a result, if the multi-electron effect is evident, a more accurate input data can not improve the validity of our method. This is indeed the limitation of most photon-ionization time delay retrieval methods.

Error effect
All the input AC patterns or streaking trace in above sections are calculated from SFA or solving TDSE. When applying our method to real experimental data, there are inherent noises from the measurement. To test the robustness of our method with respect to noise, we artificially generate the 2D spectrograms S(E, t d ) contaminated by random errors at each grid points. Here we also give two examples: TL XUV case and chirped XUV case. The two XUV pulses and IR laser fields are the same as that in section 4.2, spectrograms and corresponding AC patterns are calculated by SFA model. Starting from the Ne and Ar spectrograms, we add random noise and treat the new AC patterns from these noise-added spectrograms as the input of our retrieval. Here the noise at each data point has a mean-zero normal distribution with a standard deviation of five percent or ten percent of the original value. Figure 11 demonstrates that our method is stable for random errors up to ten percent, and therefore it can actually be applied to real experimental data.

Conclusion
In this article, we apply PROBP AC method which was previous proposed to reconstruct broadband attosecond pulses to retrieve photoionization time delays (difference). In this method, we successfully remove the so-called WPA which would introduce a prominent error to the retrieved photoionization time delays (difference) if the attosecond pulses have a large atto-chirp. We simulate the input streaking traces and corresponding AC patterns for photoelectrons from Ne (2s and 2p channels) and Ar (3p channel) by applying SFA model and solving 3D TDSE. The photoionization group delay (difference) retrieved from PROBP AC has been compared with the result from PROBP and FROG-CRAB which have a limitation in chirped XUV cases. Our results indicate that PROBP AC is a more accurate and stable method than PROBP and FROG-CRAB. The angular dependent photoionization time delays can also be retrieved by our method. Since the attosecond pulses used to generate the streaking trace have atto-chirps that cannot be fully compensated by materials over a broad spectral range, the previous reported quantity of photoionization time delays (difference) in atoms need to be re-examed carefully.