Nondipole effects on the double-slit interference in molecular ionization by xuv pulses

The double-slit interference in single-photon ionization of the diatomic molecular ion $\mathrm{H}_2^+$ is theoretically studied beyond the dipole approximation. Via simulating and comparing the interactions of the prealigned $\mathrm{H}_2^+$ and the hydrogen atom with the xuv pulses propagating in different directions, we illustrate two kinds of effects that are encoded in the interference patterns of the photoelectrons from $\mathrm{H}_2^+$: (i) the photon-momentum transfer and (ii) the finite speed of light. While both effects could modify the maxima of the interference fringes, we show that the former one hardly affects the interference minima. Our results and analysis show that the interference minima rule out the influences of the photon-momentum transfer and, potentially, the multielectron effect, thus performing a better role in decoding the zeptosecond time delay for the pulse hitting one and the other atomic centers of the molecule.


I. INTRODUCTION
The development of few-cycle infrared laser pulses, attosecond pulses, and free-electron lasers in the past two decades has enabled scientists to detect and control the ultrafast electronic and nuclear dynamics on the attosecond time scale (see [1] and the references therein). While our understanding of the dynamic mechanisms of single atoms and molecules driven by strong fields is getting deeper and deeper, researchers have already been seeking new challenges in strong-field physics. On one hand, for instance, efforts have been made in studying the high harmonic generation in chiral molecules [2], liquids [3], and solids [4] as well as the photoemission process in condensed matter [5] and nano-particles [6], bringing us unprecedented insights into the strong-field phenomena in complex systems. On the other hand, scientists have been exploring the ever-faster responses of atoms and molecules below one attosecond [7].
Recently, the zeptosecond (10 −21 s) time delay between the ionization bursts from two centers of a diatomic molecule was measured [8]. In the experiment, the 800-eV xuv laser pulse was applied to trigger the singlephoton ionization of the hydrogen molecule. Due to the finite speed of light, it takes time for the pulse to hit one and the other nuclei of the molecule, leading to a phase difference between the ionizing wave packets from two centers and eventually resulting in the distorted interference patterns in the photoelectron momentum distributions (PMDs). In principle, the information of the time delay could be extracted from the tilted angle of the double-slit interference maximum (theoretical details will be discussed in the discussion section). Yet, a noticeable deviation was found between the measurement and the theoretical prediction for the time delay [8]. The deviation, as well as the cutting-edge experimental research itself, has immediately attracted broad interest and discussion [9][10][11]. * liukunlong@hust.edu.cn The nondipole effect in molecular strong-field ionization has been previously studied in theory [12][13][14], mainly focusing on the photoelectron momentum offset in the light propagation direction. Such momentum-transfer effect takes place in atomic ionization as well, and it will be referred to as the single-atom nondipole effect in the following discussion. In this paper, we revisit the double-slit interference in the single-photon ionization of the diatomic molecules [15] by numerically solving the three-dimensional time-dependent Schrödinger equations beyond the dipole approximation. We focus on the two-center-interference nondipole effect induced by the time delay for the pulse hitting one and the other nuclei. By comparing the PMDs for the hydrogen atom and the hydrogen molecular ion, we illustrate how the photon recoil affects the interference patterns when the pulse is travelling at different directions. In particular, our results and analysis show that the interference minima are associated with the time delay for the pulse travelling from one to the other nuclei, but they would be hardly affected by the single-atom nondipole effect. We demonstrate that applying the two neighboring minima of the zeroth-order interference maxima is more accurate for extracting the zeptosecond time delay between the ionization bursts from two atomic centers of the molecule.

II. NUMERICAL METHODS
We solve numerically the nondipole three-dimensional (3D) time-dependent Schrödinger equations (TDSEs) in Cartesian coordinate system for the H + 2 prealigned along the x axis as well as for the hydrogen atom. The xuv laser pulse is linearly polarized in the y direction and propagates in the direction of e k = cos βe x + sin βe z , where β is the angle between the propagation direction and the x axis, as shown in the sub-diagram in Fig. 1. The TDSE is given by (in atomic units) [16]  where V 0 (r) is the Coulomb potential of the chosen system. The vector potential of the xuv laser pulse is given by for 0 ≤ η ≤ ωN T and otherwise A(η) = 0, with η = ω[t − (x cos β + z sin β)/c], where c is the speed of light in vacuum. Here, ω, T = 2π/ω, E, and N indicate the laser frequency, the optical cycle, the electric field amplitude, and the number of the optical cycles of the full pulse, respectively. In the present simulations, the laser parameters are chosen as ω = 29.40 a.u. (corresponding to the photon energy of 800 eV), E = 3 a.u., and N = 50. The nondipole TDSE is numerically solved using the split-operator spectral method [17] with modifications, in order to deal with the space-dependent vector potential. The details of the procedure are following. The analytic solution of Eq. (1) is given by a Dyson's time ordering operator P as (in the Coulomb gauge) with the momentum operatorp = −i∇ = (p x ,p y ,p z ). The approximate evaluation of the exponential with time ordering in Eq. (3) can be written as [20] Ψ(r, t + δt) ≈ exp −i δt 2 with t = t + δt/2. The first exponential operation on the wave function can be calculated directly and one can obtain a temporary wave function. Then, the y dimension of the temporary wave function is Fourier transformed to the momentum space as Ψ(x, y, z) → Ψ(x, p y , z), so that the second operation can be calculated by direct multiplication. In the next step, the other two dimensions of the wave function are transformed to the momentum space as Ψ(x, p y , z) → Ψ(p x , p y , p z ), for the calculation of the third operation. Finally, the rest of the operations can be analogously evaluated via transforming the temporary wave function back to the position space step by step. By repeating the processes above, one can obtain the evolution of the wave function numerically.
The initial stationary wave functions are obtained by the imaginary-time propagation method. The real-time propagation includes two parts: the interaction and the free propagation afterwards. For the interaction part, the evolution starts when the xuv pulse enters the box of the 3D grid and ends when the tail of the pulse lefts the box. We have chosen a 3D grid large enough to contain the majority of the ionizing wave packets until the interaction process ends. This is guaranteed and has been verified Illustration of the coordinates in our simulation. The molecule is prealigned along the x axis. The xuv laser pulse is linearly polarized in the y direction and propagates in the direction given by the angle β. The angular distributions of the photoelectrons are calculated as a function of θ (0 ≤ θ ≤ π) and φ (0 ≤ φ < 2π) indicated in the figure. in our calculations. After the interaction, we continue the evolution of the wave function without external fields and apply the absorbing potential to split the outgoing wave packets. In this case, one can solve the TDSE and obtain the photoelectron momentum distributions in the conventional way [18,19]. Regarding the simulation parameters, there are 2496 × 2496 × 520 grid points of the box. The spacing steps are ∆x = ∆y = 1/6 a.u. for R = 2 and 10 a.u., ∆x = ∆y = 0.14 a.u. for R = 1.4 a.u., and ∆z = 0.25 a.u. The spacing steps are chosen differently for the given internuclear distances R, so that the nuclei are placed in the middle of two grid points in the x axis to avoid the singularity of the Coulomb potential [19]. The time step for the evolution of the wave function is chosen as δt = 0.005 a.u.

III. RESULTS AND DISCUSSION
We have calculated the 3D photoelectron momentum distributions Y 3D for the hydrogen atom and H + 2 interacting with the xuv pulses. We show in Fig. 2 of the photoelectrons (see Fig. 1 for the definition of θ and φ). We also show in the figures the corresponding 2D distributions, Y xy (p x , p y ) on the bottom, Y yz (p y , p z ) on the left-hand site, and Y xz (p x , p z ) on the right-hand site, which are integrated over the third dimensions of Y 3D (p x , p y , p z ), respectively. As shown in Fig. 2(a), the peak of the distribution Y xz drifted towards the light propagation direction, due to the momentum transfer from the photon to the photoelectron. For the cases of H + 2 shown in Figs. 2(b) and 2(c), the double-slit interference patterns appear. The interference fringes shown by Y xz are perpendicular to the molecular axis for both cases, but the distributions of the fringes seem different.
To quantitatively compare the distributions of the interference patterns, we calculated the photoelectron angular distributions (PADs) given by for the atom and H + 2 (R = 1.4 a.u.) in the cases of β = 0, 45 • , and 90 • , respectively. The results after normalization to the maximum of Y + θ are shown in Figs. 3(a)-3(c). We can see that PADs for the interference patterns are enveloped by those PADs for the hydrogen atom. It indicates that, when the xuv pulse hits each center of the molecule, the ionizing wave packet burst is approximately equivalent to that from a single atom and two wave packets coherently interfere with each other in the continuum. Thus, the distorted PAD from each atomic center will be encoded in the interference patterns, modifying the locations and yields of the fringe maxima. This can be seen in both Fig. 2 and Fig. 3.
When the pulse travels at β = 90 • (i.e., perpendicular to the molecular axis), the photon-momentum transfer leads to more photoelectron yield for p z > 0 than that for p z < 0. The global maxima of Y + θ and Y − θ are no longer equal, although they both appear at θ = 90 • for the atom and the molecule, as shown in Fig. 3(c). In this particular case, the pulse interacts with two atomic centers at the same time. The situation becomes complicated when the pulse no longer travels perpendicular to the molecular axis. As shown in Figs. 3(a) and 3(b), the PADs for the atom are tilted towards the light direction and they peak at the angle (referred to as θ A ) far from θ = 90 • . For the molecule, the interference maximum is also tilted but peaks at the angle (referred to as α 0 ) slightly smaller than 90 • . At a first glance, the tilted interference maximum is affected by the photon-momentum transfer during the interaction. This is true, but it was shown that the finite speed of light plays a more important role in modifying the interference patterns [8]. In details, beyond the dipole approximation, there is a time delay for the pulse interacting with one and the other atomic centers of the molecule due to the finite speed of light, resulting in a phase difference between the ionizing wave packets generated from two centers. Such phase difference eventually modifies the double-slit interference pattern of the photoelectrons.
It was proposed that by finding the interference maxima, one can obtain the interacting time delay that is as small as hundreds of zeptoseconds [8]. However, as illustrated in Fig. 3 and also shown by previous works [9], the photon-momentum transfer has impact on the interference maxima. Moreover, in multielectron systems, the multielectron dynamics would potentially affect the interference patterns, as discussed in [10]. These factors might limit the accuracy of extracting the time delay from the interference maximum. In the following, we will show that adopting the interference minima instead of the maxima will rule out the side effects and provides us with a better way to decode the interacting time delay.
We start with the wave function of the photoelectron ionizing from the atom, which is written as where A(p) describes the amplitude of the yield distribution and E k = p 2 /2 indicates the final kinetic energy of the photoelectron. In Eq. (8), we have made two assumptions: (A1) the electronic wave packet originates from the origin and (A2) the Coulomb effect on the spreading wave packet is neglectable. By further assuming that the ionizing wave packet generated from each center of the molecule is equivalent to that from a single atom (referred to as the assumption A3), we can write the wave function for the photoelectron emitted from the left and right atomic centers, respectively, as with r 1,2 = r ± R/2 and R = (R, 0, 0). In Eq. (9), t d indicates the interacting time delay. The term ωt d is the extra phase difference between the wave packets generated from the left and right atomic centers due to the time delay. This phase difference can be understood as following. We assume that the xuv pulse interacts with the atom on the left side earlier. Then, by absorbing the photon energy, the total energy of the wave packet about to be ionized from the left core is lifted by ω, while that on the right side remains unchanged. When a duration of t d passes, an initial phase difference of ωt d between two sides arises. After the interaction, the ionizing wave packets from both sides experience the same process: carrying the same energy, escaping from the parent core, and ending up with the same kinetic energy, which will accumulate the same phase for both wave packets. Thus, the overall phase difference (regarding the energy part) between Ψ L and Ψ R is ωt d .
Since we have E k = ω − I p in single-photon ionization, with I p being the ionization potential, we define E 0 = E k + I p ≡ ω as the initial kinetic energy of the photoelectron. Then, the wave function of the photoelectron from the diatomic molecules is written as On the righthand site of Eq. (11), the first line equals to Ψ A (p), i.e. the photoelectron wave function from a single atom. The second one simply introduces an extra global phase to the wave function. The last one is the interference term which modulates the PADs. In principle, the interacting time delay can be decoded from the interference pattern as t d appears in the interference term. However, it could be tricky for measurements. Let us look into the conditions for the extrema of the interference. The interference term leads to its maxima when the following condition is satisfied: For the given photoelectron momentum p = 2(ω − I p ) with I p being the ionization potential of the system, the interference maxima can be found at the angles θ n satisfying In particular, for the zeroth-order interference maximum (i.e., n = 0), we define the tilt parameter as by assuming t d = (R cos β)/c. Note that the conditions shown above only guarantee the maxima of the interference term in Eq. (11). The global maximum of the PAD relies on the envelope A(p) as well. In general cases, A(p) peaks beyond θ 0 , unless the pulse travels perpendicular to the molecular axis, as already shown in Fig. 3. Therefore, the angle α 0 extracted from the maximum of the observed interference fringe is expected to differ from θ 0 . To verify our expectation, we compare the values of T Theo and cos α 0 (extracted from the simulated PADs) at three internuclear distances in the case of β = 0. Note that for the calculation of T Theo , we have assumed E 0 = ω, i.e., the initial kinetic energy of the photoelectron equals to the photon energy. The results are shown in Fig. 4. It is noticeable that cos α 0 tends to deviate from T Theo more significantly as the internuclear distance becomes smaller. This can be intuitively understood from Fig. 5, where we show the PADs obtained from numerical solutions of the TDSEs for R = 1.4 and 10 a.u. and β = 0. For larger R, the frequency of the modulation of the interference fringes becomes higher. In this case, the interference fringes are so narrow that the tilted envelope hardly modifies the interference maxima. For R = 1.4 a.u., however, the broad distribution of the interference fringes would be  affected more significantly, leading to the deviation of cos α 0 from T Theo .
To rule out the effect of the envelope A(p), we turn to the interference minima. According to Eq. (11), the condition for the interference minima is given by The two neighboring minima of the zeroth-order maxima are located at the angles θ ± given by Then, we can obtain the relation between θ ± and θ 0 as following: Theoretically, the interference term in Eq. (11) equals to zero at θ ± , so the effect of the envelope A(p) on the interference minima should be vanishing. For demonstration, we extract the angles α ± for the neighboring minima of the zeroth-order maxima from the PADs and show the tilt parameter given by (cos α + + cos α − )/2 in Fig. 4. One can see that the theoretical prediction is generally in agreement with the observation given by (cos α + +cos α − )/2. In particular, the deviation for small internuclear distances is reduced remarkably. In Figs. 4, we notice that T Theo is slightly higher than the tilt parameter given by (cos α + + cos α − )/2. One of the possible reasons is that we assumed a constant initial kinetic energy equal to the photon energy in Eqs. (9) and (17). In fact, the photoelectron experiences a decelerating process even in the duration of t d , and the deceleration would become more pronounced under the stronger Coulomb attraction. As a result, the theoretical T Theo has been overestimated by adopting E 0 = ω in Eq. (14), especially for smaller internuclear distances. Besides, the assumptions A1-A3 mentioned in the previous discussion are also the potential reasons leading to the deviation between the theoretical prediction and the numerical calculations. Nevertheless, the assumptions A1-A3 are fairly justified as the deviation is less than 1.6%.
Furthermore, the theoretical model given by Eq. (14) indicates that the tilt parameter is linearly proportional to cos β. In Fig. 6, we show the dependence of the tilt parameter on the direction of the light propagation. We can see that the tilt parameters obtained from the interference minima of the PADs agree with the theoretical results very well for all three internuclear distances.
So far, the present work based on the single-activeelectron system does not reproduce the results measured in the experiment [8], even if we look at the tilt parameters given by the interference maxima (see Fig. 6). On the bright side, however, we have shown that the interference minima of the observable PADs perform a better role in extracting the information regarding the interacting time delay, as they mathematically rule out the influence from the photon-momentum transfer. The remaining question is whether the multielectron effect could be ruled out as well at the interference minima. To answer this question, let us revisit the theoretical model given by Eq. (11). On one hand, the photonmomentum transfer relies on the Coulomb interaction of the photoelectron with the ion, as shown by previous studies [12]. If we further consider an additional active bound electron remaining in the ion, then the electronelectron Coulomb interplay would most likely modify the momentum transfer process. This may lead to the further modification of the overall angular distribution of the photoelectron given by A(p).
This would potentially affect the tilt parameters corresponding to the interference maxima. However, as shown in our previous discussion, it can be safely ruled out if we adopt the minima to calculate the tilt parameters. On the other hand, the active electron remaining in the ion would change the deceleration of the photoelectron and thus might affect the interference term in Eq. (11). To estimate the impact, we may compare it to the case where an additional nucleus is present or not. For instance, the Coulomb attraction on the electron escaping from the first core at the beginning of the ionization for R = 1.4 a.u. is almost twice stronger than that for R = 10 a.u. Yet, in both cases the deviation of the tilt parameters between the theoretical model and the numerical simulations (for the minima) is rather small, as shown in Fig. 4. Therefore, by qualitatively analyzing the electron-electron Coulomb interaction, we expect that the multielectron effect would have hardly impact on the locations of the interference minima.

IV. CONCLUSION AND OUTLOOK
In summary, we have shown how two kinds of nondipole effects modify the double-slit interference in single-photon ionization of H + 2 when the xuv pulse travels in different directions. The single-atom nondipole effect leads to a tilted envelope of the overall interference pattern, whereas the interacting delay between the pulse and two nuclear centers further modifies the locations of the interference maxima and minima. Our results show that the locations of the interference minima are hardly changed by the single-atom nondipole effect. Further analysis indicates that the multielectron dynamics is likely to have insignificant effect on the interference minima as well. By adopting the two neighboring minima of the zeroth-order interference maxima to calculate the tilt parameters, we found good agreement between theoretical model and numerical simulations regarding the time delay between the pulse interacting with two centers of the molecule.
Yet, the present work, as well as the previous theoretical study having included the multielectron dynamics [10], fails to reproduce the experimental results in [8]. In the framework of Schrödinger equation, the laser field is treated classically and we can investigate how the given laser field induces the electronic dynamics. However, how would the electronic motion impact the imposing photon during the interaction? From another perspective, there is electron cloud between two nuclei of a molecule, so is the phase velocity of the external field still the speed of c inside the molecule? The answers to these questions are unknown yet. Further study is in progress as these questions might be the key to the gap between the measurement and the simple theoretical model regarding the zeptosecond time delay.

ACKNOWLEDGMENTS
Liu thanks Barth and Renziehausen for interesting discussion. This work is supported by the National Natural Science Foundation of China (12174133, 92150106, 12021004).
The computing work in this paper is supported by the public computing service platform provided by Network and Computing Center of HUST.