Effect of the finite speed of light in ionization of extended molecular systems

We study propagation effects due to the finite speed of light in ionization of extended molecular systems. We present a general quantitative theory of these effects and show under which conditions such effects should appear. The finite speed of light propagation effects are encoded in the non-dipole terms of the time-dependent Shrödinger equation and display themselves in the photoelectron momentum distribution projected on the molecular axis. Our numerical modeling for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}_{2}^{+}$$\end{document}H2+ molecular ion and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {Ne}_2$$\end{document}Ne2 dimer shows that the finite light propagation time from one atomic center to another can be accurately determined in a table top laser experiment which is much more readily accessible than the ground breaking synchrotron measurement by Grundmann et al. (Science 370:339, 2020).

Every quantum system evolves on its characteristic time scale which varies widely between molecules (femtoseconds-10 −15 s) 1 , atoms (attoseconds-10 −18 s) 2 and nuclei (zeptoseconds-10 −21 s) 3 . A crossover between these time scales is very rare in nature. It was therefore quite unexpected to discover an ionization process in the hydrogen molecule that evolved on a zeptosecond time scale 4 . An explanation of this phenomenon appeared to be quite simple. While it takes tens of attoseconds for an electron to trespass the H 2 molecule, the incoming light wave sweeps from one molecular end to another orders of magnitude faster. This results in one of the constituent hydrogen atoms getting ionized a fraction of the attosecond sooner than its counterpart. Such a tiny ionization delay manifests itself quite noticeably in the two-slit electron interference that the H 2 molecule readily displays 5 . To discover a zeptosecond delay in molecular photoionization in 4 the authors needed to deploy an extremely bright synchrotron source of highly energetic photons 6 . We demonstrate that even a table top laser experiment is capable of detecting a similar effect making it much more readily affordable.
In this work we present a general quantitative theory of the delay due to the finite speed of light propagation and we show under which conditions such effects should manifest themselves. In our numerical demonstrations, we consider the H + 2 molecular ion and the Ne 2 dimer. The H + 2 molecular ion has been scrutinized since the early days of quantum mechanics 7 and was recently used as a model for the study of the interference effects in photon-momentum transfer for the process of molecular ionization 8 . Photoionization studies of Ne 2 is a novelty 9 . We subject both targets to an attosecond laser pulse that can be readily produced in high-order harmonics generation sources 10,11 . The photoelectron flux encodes the timing information about the ionization process. This flux is reconstructed by solving the laser-driven time-dependent Schrödinger equation (TDSE). The numerical results obtained for H + 2 using the TDSE can be interpreted in a transparent qualitative way by considering a very simple heuristic tight-binding model (TBM). This gives us a tool for understanding the time delay caused by the finite speed of light propagation. We apply this tool to H + 2 and Ne 2 and find the time delay of a fraction of attosecond that depends sensitively on the orientation of the molecular axis relative to the propagation and polarization directions.
Atomic units (a.u.) are used throughout the paper with the electron charge e, mass m and the Planck constant all set to unity m = e = = 1 . The speed of light in these system of units is c ≈ 137.036 a.u. www.nature.com/scientificreports/ where Ĥ mol is the field-free one-electron Hamiltonian and Ĥ int (t) describes the field-target interaction. We consider the hydrogen molecule ion H + 2 in the ground state as a target. The ion interacts with a short (total duration of four optical cycles) pulse with the base frequency ω = 4.04 a.u. (photon energy of 110 eV) with the peak field strength E 0 = 0.1 a.u. (intensity of 3.51 × 10 14 W/cm 2 ). The pulse described by the vector potential A(t − x/c) propagates in the positive x-direction and is polarized in the z-direction.

Scientific
We apply the procedure previously used in 12 to study non-dipole effects in atomic photoionization. The leading order relativistic corrections to Ĥ int (t) come from the linear term of the expansion of A(t − x/c) in powers of c −1 . Such an expansion leads to the following Hamiltonian describing the target-field interaction: More detailed description of the derivation of the Hamiltonian (2) and details of the numerical procedure we used to solve the TDSE (1) can be found in the "Methods" section. The interaction Hamiltonian (2) can be related by a gauge transformation to the Kramers-Henneberger Hamiltonian used in 13 .
An additional source of relativistic corrections in Ĥ int is the interaction of the magnetic field of the pulse with the electron spin. These corrections cannot be obtained by a simple generalization of the minimal coupling Hamiltonian. To obtain them, one needs to consider systematically the transition of the Dirac equation to the non-relativistic limit 14 . The inclusion of the electron spin, however, is not necessary for the present study where we keep only the c −1 terms as was done, for instance, in 12,15 . The Breit-Pauli relativistic corrections to Ĥ mol , such as the effects of the relativistic kinematics, spin--orbit interaction and the so-called Darwin term, are all of the order of c −2 and can also be safely omitted 16 .
For a weak electromagnetic field that we employ, we can also obtain an analytical expression for the ionization amplitude by using the perturbation theory (PT) and treating the operator (2) as a perturbation: Here φ 0 and φ − p are the initial and final molecular states with the corresponding energies ε 0 and E p . We split the ionization amplitude (3) into the nonrelativistic part a (0) p and the first order relativistic correction a (1) p . By introducing the Fourier transforms A(t) = (2π) −1 a(�)e −it� d� and A 2 (t) = (2π) −1 b(�)e −it� d� we obtain for these amplitudes: The second term on the right hand side of Eq. (5) can be neglected because the ratio of the second term and the first term in this expression is approximately b(�)/(a(�)p) , where p is the typical value of the momentum of the ionized electron. For the pulse parameters that we consider, this value can be estimated as E 0 /(ωp) ≈ 0.01.
It can be seen from Eq. (5) that the relativistic correction a (1) p vanishes in certain cases. In particular, it is zero for an axially symmetric system with the field-free Hamiltonian that is invariant under rotations about the z-axis. Indeed, for such systems the scattering state φ − p (r) is a function of the arguments z, p z , ρ, p ρ and ρ · p ρ , where z, p z and ρ , p ρ are the components of the r and p vectors parallel and perpendicular, respectively, to the z-axis of the coordinate system that we employ. For an axially symmetric system, the ground state wave function φ 0 in Eq. (5) (which we assume to be non-degenerate) is invariant under rotations about the z-axis. We obtain then from Eq. (5): a (1) (p x , p y , p z ) = −a (1) (−p x , p y , p z ) . Therefore the amplitude a (1) p is zero in the plane p x = 0 which encompasses the most important region of the momenta near the maximum of the momentum distribution.
The effect of the propagation correction (5) can be most easily illustrated within the tight-binding model (TBM) in which the ground state of H + 2 and Ne 2 is represented by a Heitler-London wave function: In the TBM, the overlap of the two terms is small and φ(r) is represented by a spherically symmetric atomiclike state. Under these conditions, and by employing the Born approximation for the scattering state φ − p , one can show (the formal derivation is given in the section "Methods") that Eq. (4) and Eq. (5) lead to the following expression for the ionization amplitude: where δ = −κ · R/2 and κ = c −1 e x is the photon momentum, and a d (p) is the dipole amplitude which we would obtain if electrons were emitted from the state described by the single-center initial wave-function φ(r) centered at the origin. It is this phase factor δ in Eq. (7) that is responsible for the modified interference pattern observed and decoded in 4 . www.nature.com/scientificreports/ Appearance of this additional phase in Eq. (7) is due to an extra propagation time of the light wave from one end of the molecule to the other. This interpretation becomes yet more transparent if we use coordinate representation for the part of the wave-function describing the ionized wave packet: Here φ − p (r) are the molecular scattering states. To evaluate this integral in the limit t → ∞ , we rely on the saddle-point method (SPM) that is commonly used in description of ionization 17,18 or scattering 19 processes. By writing a d (p) = |a d (p)|e iη(p) in Eq. (7) and employing the SPM, we obtain from Eq. (8) (the details of the derivation are given in the section "Methods"): The value p = p 0 in Eq. (9) is the momentum values satisfying energy conservation p 2 0 /2 = ε 0 + ω , τ = ∂η/∂E is the usual Wigner photoemission time-delay 20 and τ 1 = ∂δ/∂E = −R x /(2c) with δ = −κ · R/2 from Eq. (7). The τ 1 term represents an additional time-delay that it takes for the light pulse to cover the distance R x = R · e x . The function G(u) in Eq. (9) is sharply peaked near the origin. Therefore the two terms in Eq. (9) describe two wave packets emitted from the two atomic centers r = −R/2 and r = R/2 at the times τ + τ 1 and τ − τ 1 , respectively.
For transparency of derivation, we omitted in Eq. (9) the Coulomb terms which would only add slowly varying (logarithmic) corrections in the arguments of G(u) 17,18 . These additional logarithmic terms are the same for the two wave packets and they would therefore cancel in the time delay difference between these wave packets.
We note that the dipole amplitude a d (p) in the Eq. (7) is essentially a product of two factors, the factor p z responsible for the angular dependence of the amplitude and the Gaussian factor exp −b(p − p 0 ) 2 representing the energy conservation p 2 0 /2 = ε 0 + ω . Such a Gaussian representation of the ionized wave packets emitted in the single-center problems is often used in studying temporal dynamics of atomic ionization 17,18 . The parameter b determines the width of the wave packet and depends on the pulse parameters (more details are given in the section "Methods"). By employing the Gaussian ansatz for a d (p) we finally obtain from Eq. (7): We use Eq. (10) to evaluate the photoelectron emission pattern and to compare it with the fully ab initio TDSE calculations for various orientations and different inter-nuclear distances of H + 2 . This comparison is presented in Fig. 1. The top row of panels illustrates the geometry of the ionization process. It is assumed that the molecular axis is confined to the xz-plane where the propagation and polarization vectors of the pulse belong, making an angle θ with the propagation direction. In (a-b), θ = 0 while in (c) θ = π/4 . The photoelectron momentum distribution is projected on the xz plane and computed as P(p x , p z ) = |a p | 2 dp y with the amplitude a p obtained by projecting the TDSE solution on the set of the scattering states of H + 2 . The number and location of the bright spots in Fig. 1 reflect a simple two-center interference pattern governed by the cosine term in Eq. (10).

Discussion
Comparison of the TDSE calculations (the middle row of panels in Fig. 1) and results based on Eq. (10) (the bottom row of panels) shows that the TBM reproduces the spectra adequately for the geometries that we consider. The TBM certainly has some limitations if we demand an accurate global representation of the spectra. For instance, for the case of the non-collinear orientation of the polarization vector and the molecular axis shown in Fig. 1, the ab initio TDSE and the TBM spectra exhibit some differences. One could improve, in principle, the visual agreement of the TDSE and the TBM results in this case by adjusting the value of the fitting parameter b in Eq. (10). The value of this parameter affects the tilt of the TBM distribution with respect to the x-axis, so we can rotate the TBM spectrum by varying this parameter. Achieving a good global agreement between the TDSE and TBM is, however, not needed for the purposes of the present work. For the procedure of the extraction of the relativistic delays that we describe below to work, we need Eq. (10) to accurately describe the spectra locally, in the vicinity of a maximum. We will see below that this lesser goal can be achieved. We will, therefore, analyze and interpret our TDSE results using the transparent TBM that is equally applicable to both H + 2 and Ne 2 . We will focus on the photoelectron momentum distribution projected on the (xz) plane and integrated over the momentum component p ⊥ which is perpendicular to the molecular axis. Such a momentum distribution is a function of the momentum component p || which is parallel to the molecular axis. By employing the Gaussian ansatz (10) we obtain: We use the analytical expression (11) with adjustable parameters B, C and δ to fit P(p ) obtained from the numerical TDSE calculations. The accuracy of the fitting procedure is illustrated in Fig. 2 where we display P(p ) for various geometries and internuclear distances. In the cases of θ = 0 P(p ) is simply the projection of the 2D (11) P(p || ) = P(p) dp y dp ⊥ ≈ B exp −Cp 2 || cos 2 p || R 2 + δ . www.nature.com/scientificreports/ momentum distribution on the horizontal axis. We use the interval of the p � ∈ (−0.5, 0.5) for the fitting procedure. As is seen in Fig. 2, the analytical fit with Eq. (11) represents the central maximum of the TDSE calculations with the corresponding R value quite accurately. This allows us to extract the phase shift δ accumulated due to the finite speed of light propagation. In Fig. 3 we display thus extracted phase shifts δ for the molecular orientations and internuclear distances illustrated in Fig. 2. To enhance the relativistic effects, we artificially decrease the speed of light in the TDSE calculations from its physical value c 0 = 137.036 a.u. down to c = c 0 /10 . According to TBM, the phase shift should scale linearly as δ = αc 0 /c with the slope α = −R x κ/2 . The predicted linear scaling is confirmed very accurately by the numerical values shown in Fig. 3 with only a small error margin. The time delay values corresponding to the physical speed of light c 0 are shown in Table 1 in comparison with the estimate �t = 2R x /c 0 provided by the TBM. Agreement of the results can be deemed quite satisfactory given the relative simplicity of TBM. More importantly, the linear dependence of the TDSE results on the parameter c 0 /c clearly demonstrates existence of the finite speed of light effect in ionization of diatomic molecules.

Conclusions
In conclusion, our work was motivated by the synchrotron based experiment 4 which discovered a zeptosecond time delay in photoionization of the H 2 molecule. We aimed to demonstrate that a similar delay, that is caused by the finite speed of light propagation from one constituent atom to another, can be detected in a much more accessible table top laser settings. In doing so, we developed a general theory of the finite speed of light propagation effects in ionization of extended systems. As a simple case study, we considered the H + 2 molecular ion interacting with a short laser pulse. This target affords a very accurate numerical treatment within the timedependent Schrödinger equation. At the same time, an heuristic tight-binding model employing the Heitler-London molecular ground state produces very similar results. TBM can be easily adopted to the Ne 2 by a simple increase of the inter-atomic distance to R ≃ 6 a.u. Notably, the corresponding interference pattern displayed in Fig. 1b) is very similar to that obtained in the experiment 9 and the earlier SPM modeling 21 .
Our simulations demonstrate that the speed of light delay in ionization of diatomic molecular targets can be decoded from the photoelectron momentum distribution projected on the molecular axis. This method has a clear advantage over the earlier synchrotron measurement based on decoding the 2D interference pattern 4 . Indeed, the bright interference spots have a finite angular width. To detect their shift due to the finite speed of light requires a sufficiently large photon momentum that should not be vanishingly small in comparison with the photoelectron momentum. This in turn requires very high photon energy (800 eV in 4 ). Fitting of a   www.nature.com/scientificreports/ one-dimensional momentum distribution produces significantly reduced error bars. Thus smaller values of the photon energy of the order of 100 eV can be used. Driving pulses with the parameters employed in our simulations (110 eV photon energy, intensity of 3.51 × 10 14 W/cm 2 ) can be obtained by using relativistic high harmonic generation. It was shown 22 that even the Schwinger intensity ( 10 29 W/cm 2 ) can be reached by reflecting the currently available PW laser beam with intensity of 10 22 W/cm 2 on a plasma mirror. Intensities of the order of 10 21 W/cm 2 can be reached by using a table-top 100 TW system 23 . In addition, by projecting the momentum distribution on the molecular axis, we increase the count rate and improve statistics of the measurement. Hence the photon flux can be significantly reduced. This reduction of both the photon flux and energy makes the proposed method much more readily accessible in desk-top conventional laser settings. This we hope will stimulate further speed of light delay determinations in diatomic molecules and other extended systems. Incidentally, the finite c delay measured in H 2 deviates from accurate theoretical modeling 24 . The proposed novel technique may shed an additional light on the nature of this disagreement.

Methods
Leading order relativistic corrections to the interaction Hamiltonian and details of the numerical procedure. The leading order relativistic corrections to the non-relativistic Hamiltonian describing the target-field interaction arise from the vector potential being a traveling wave A(η) with η = t − x/c (we assume that the pulse propagates in the positive x-direction and is polarized in z-direction). The leading relativistic correction to non-relativistic interaction Hamiltonian can be obtained by substituting this expression for the vector potential into the standard minimal coupling Hamiltonian Ĥ min = 1 2 (p + e z A(x, t)) 216 , and keeping the terms linear in c −1 . One arrives then at the expression for the interaction Hamiltonian 12 : where A(t) no longer depends on the coordinates, and E(t) = − ∂A(t) ∂t is the electric field of the pulse. The last term on the r.h.s. of the Eq. (12) is a function of time only, and can therefore be removed by a unitary transformation of the wave-function.
To solve the TDSE with the interaction Hamiltonian (12) we use the procedure that we employed previously for the solution of the non-relativistic TDSE [25][26][27] . Solution of the TDSE is represented as an expansion: where Y lm (θ) are the spherical harmonics. The radial variable is treated by discretizing the TDSE on a grid with the step-size δr = 0.1 a.u. in a box of the size R max . The values of the parameters R max and l max in Eq. (13) were chosen (after the necessary convergence checks) as: l max = 10 , R max = 400 a.u.
For the vector potential in Eq. (12) we used the form: with ω = 4.04 a.u. (photon energy of 110 eV), peak field strength E 0 = 0.1 a.u. (intensity of 3.51 × 10 14 W/cm 2 ), and total duration T 1 = 4T , where T = 2π/ω -optical cycle corresponding to the base frequency ω. The potential energy in the field-free H + 2 Hamiltonian: was expanded in spherical harmonics. Using this expansion, Eq. (13) of the wave-function, and expressing all other operators in Eqs. (12) and (15) as irreducible tensor operators 16 , one obtains a system of coupled equations for the functions f lm (r, t) in Eq. (13). This system was solved using the Lanczos algorithm 28 (we used propagation stepsize t = 0.05 a.u., and propagation algorithm of the order N = 30 ). The initial ground state of the H + 2 was prepared using the relaxation procedure 29 , propagating variational estimate for the initial wave-function in the imaginary time. Differential ionization probabilities P(p) are calculated by letting the wave-function of the system evolve for some time after the end of the pulse so that the wave-packet describing ionized electron state leaves the region where r < R 0 with some R 0 (we use the value R 0 = 80 a.u. in the calculations). In the region r > R 0 the true scattering states of H + 2 can be well approximated by the single-center Coulomb scattering states, and to calculate the differential ionization probabilities we can project the TDSE solution on the ingoing Coulomb scattering states, calculating all the the radial integrals starting from r = R 0 . This procedure proposed in 30 helps to avoid possible non-orthogonality of molecular states belonging to bound and continuous spectra which is difficult to avoid in numerical calculations 31 . We checked that results we obtain are stable against further increase of the cutoff parameter R 0 .
Details of the tight binding model. A general expression for the amplitude we obtained was: www.nature.com/scientificreports/ with and where � = E p − ε 0 , and a(�) , and b(�) are Fourier transforms of A(t) and A 2 (t) . These expressions can be simplified if we assume first that the momentum is large enough so we can replace the scattering state φ − p with the plane wave, and second, if we use a tight binding description for the bound state of H + 2 , i.e., we assume initial state wave-function can be represented as φ 0 (r) = (φ(r − R/2) + φ(r + R/2))/ √ 2 , where support of φ(r) is some (sufficiently small) ball centered at the origin. We assume that φ(r) is spherically symmetric, i.e. describes an s-state. It is easy to see that under these approximations: and In Eq. (19): is the dipole amplitude which we would obtain if electrons were emitted from the single-center initial wavefunction φ(r) centered at the origin.
From Eqs. (16), (19), and (20) we obtain: Origin of the Gaussian factor in Eq. (10). Derivation of the Eq. (10) was based on the assumption that the dipole amplitude a d (p) in the Eq. (7) can be represented as a a product of two factors, the factor p z responsible for the angular dependence of the amplitude and the Gaussian factor exp −b(p − p 0 ) 2 , peaked around the expressing energy conserving momentum p 2 0 /2 = ε 0 + ω . Possibility of such a representation can be seen if we examine a perturbative expression for a d (p) , which looks quite similar to the Eq. (4) above, with the difference that we should use a single-center wave-function φ(r) to calculate the matrix element: where � = E p − ε 0 , and a(�) is the Fourier transform of the pulse vector potential. Assuming that the scattering state in Eq. (23) can be approximated by a plane wave, and assuming for the moment Gaussian shape A(t) = A 0 e −βt 2 sin ωt for the pulse vector potential we obtain from Eq. (23): where φ(p) is Fourier transform of the single-center wave-function φ(r) . We do not write explicitly unimportant constant factors in Eq. (24), and we omitted in Eq. (24) exponentially small term proportional to . For long enough pulse (small value for the β-parameter) exponential factor in Eq. (23) is sharply peaked about the energy conserving value of E p , and we can neglect the much slower momentum dependence of the Fourier transform of φ(r) (since φ(r) is assumed to represent an s-state its Fourier image φ (p) depends only on |p| ). For the same reason we can rewrite the exponential factor in Eq. (24) as: with p 0 such that p 2 0 /2 = ε 0 + ω , is the energy conserving value of momentum, and g(p) is a function which varies much slower than the exponential factor in the region where the total expression is non-negligible, so that we can neglect its momentum dependence. That gives us the expression for the amplitude a d (p) we used above. In the present calculation we used the pulse shape (14) different from the Gaussian shape, but one can see that near the photo-electron distribution maximum (p close to p 0 ) we can still use the Gaussian shape for the amplitude (this corresponds again to assuming that we neglect p-dependence of all slowly varying functions near p = p 0 ). www.nature.com/scientificreports/ Derivation of Eq. (9). To derive Eq. (9) we substitute expression (10) for the amplitude into Eq. (8), and assume that at large distances, which interest us, the scattering states can be approximated by the plane waves. As we mentioned above, this assumption is not true if the Coulomb interaction is present. We make it only to simplify the derivation. The Coulomb interaction would add slowly varying logarithmic corrections which could easily be included into consideration. These Coulomb corrections, however, are not important for us as long as we are interested in the time delay difference between the wave packets. Representing the dipole amplitude as a d (p) = |a d (p)|e iη , we obtain then the following integral for the wave-function describing ionized wave-packet (ignoring the unimportant constant factors): where and we use notation I 1 , I 2 for the integrals with the exponential functions e iS 1 (p,r,t) and e iS 1 (p,r,t) , correspondingly. In Eq. (26) η is the phase of the dipole amplitude a d (p) and δ = −κ · R/2 -is the phase due to the non-dipole effects from Eq. (7). For large t the integrals in Eq. (25) can be computed using the saddle point method. Computing derivatives with respect to p and taking into account that both δ and η (for the case of the single-center wave-function φ(r) describing an s-state) are functions of |p| , we can write the saddle point equation for the integral I 1 as follows: where p sp is the saddle point momentum, τ = ∂η/∂E is the Wigner photoemission time-delay 20 and τ 1 = ∂δ/∂E = −R x /(2c) is the time delay due to the relativistic non-dipole effects. Then, for t → ∞ we obtain for the integral I 1 : where p sp is the root of the saddle-point equation (27). As we mentioned above, the dipole amplitude a d (p) is essentially a product of two factors, the factor p z giving the angular dependence of the amplitude and the factor depending on the momentum magnitude |p| , which is peaked around the momentum value p = p 0 satisfying the energy conservation p 2 0 /2 = ε 0 + ω . We can write therefore, using the expression (27) for p sp : where we absorbed all the constant factors in the definition of the function G. Basing on the properties of the a d (p) we discussed above, we may conclude that the function G(u) in Eq. (29) is sharply peaked around the value u = 0 . It describes, therefore, a wave-packet propagating from the point r 0 = −R/2 and delayed by the amount τ + τ 1 . For the integral I 2 in Eq. (25) we obtain quite analogously: with similar interpretation. Adding contributions (29) and (30) we obtain Eq. (9).