Attosecond delays of high harmonic emissions from isotopes of molecular hydrogen measured by Gouy phase XUV interferometer

Mumta Hena Mustary Centre for Quantum Dynamics, Griffith University, Brisbane, Queensland 4111, Australia. Liang Xu Key Laboratory for Laser Plasmas (Ministry of Education) and School of Physics and Astronomy, Collaborative Innovation Center for IFSA (CICIFSA), Shanghai Jiao Tong University, Shanghai 200240, China. Shanghai Key Lab of Modern Optical System, University of Shanghai for Science and Technology, 200093 Shanghai, People’s Republic of China. Wanyang Wu Key Laboratory for Laser Plasmas (Ministry of Education) and School of Physics and Astronomy, Collaborative Innovation Center for IFSA (CICIFSA), Shanghai Jiao Tong University, Shanghai 200240, China. Nida Haram Centre for Quantum Dynamics, Griffith University, Brisbane, Queensland 4111, Australia. Dane E. Laban Centre for Quantum Dynamics, Griffith University, Brisbane, Queensland 4111, Australia. Han Xu Centre for Quantum Dynamics, Griffith University, Brisbane, Queensland 4111, Australia. Feng He Key Laboratory for Laser Plasmas (Ministry of Education) and School of Physics and Astronomy, Collaborative Innovation Center for IFSA (CICIFSA), Shanghai Jiao Tong University, Shanghai 200240, China. CAS Center for Excellence in Ultra-intense Laser Science, Shanghai, 201800, China. Igor V. Litvinyuk Centre for Quantum Dynamics, Griffith University, Brisbane, Queensland 4111, Australia. R. T. Sang Centre for Quantum Dynamics, Griffith University, Brisbane, Queensland 4111, Australia.

High harmonic spectroscopy[1, 2, 3] can access structural and dynamical information on molecular systems encoded in amplitude and phase of high harmonic generation (HHG) signals[4, 5, 6, 7, 8, 9]. However, measurement of the harmonic phase is a daunting task. Here we present a precise measurement of HHG phase difference between two isotopes of molecular hydrogen using the advanced extreme-ultraviolet (XUV) Gouy phase interferometer[10]. The measured phase difference is about 200 mrad, corresponding to $\sim$ 3 attoseconds (1 as = 10^-18 s) time delay which is nearly independent of harmonic order. The measurements agree very well with numerical calculations of a four-dimensional time-dependent Schödinger equation. Numerical simulations also reveal the effects of molecular orientation and intra-molecular two-centre interference on the measured phase difference. This technique opens a new avenue for measuring the phase of harmonic emission for different atoms and molecules. Together with isomeric or isotopic comparisons it also enables the observation of subtle effects of molecular structures and nuclear motion on electron dynamics in strong laser fields.

The HHG process is a sensitive probe of molecular dynamics and structures [11]. The hydrogen molecule, being the lightest, exhibits the fastest nuclear motion. Also, being the simplest neutral molecule, it allows accurate ab initio quantum mechanical simulations without resorting to severe approximations. Those simulations can be used to understand and validate the experimental results. Moreover, availability of different isotopes adds isotopic comparison as a valuable benchmarking and validation tool. A number of studies with such comparisons have already been performed, including effects of nuclear dynamics on relative HHG yields and tunnel ionization rates in H₂ and D₂. The strong field ionization of H₂ launches a free electron with H ${}_{2}^{+}$ in the $1s\sigma_{g}$ state. The electron wavepacket driven by the laser field accelerates away from the parent ion. At the same time, the nuclear wavepacket evolves on the potential energy surface of H ${}_{2}^{+}$ (see Fig. 1(a)). After a certain time delay, the electron wavepacket returns to the parent ion upon the reversal of the driving laser field. Its consequent recombination induces an oscillating dipole, which imprints its amplitude and phase on the emitted HHG photons.

The investigation of the nuclear dynamics of hydrogen isotopes was envisioned theoretically [12, 13, 14] and implemented experimentally[6, 7, 15] by comparing their high-harmonic yields. However, only two studies explored the isotope effects on the high-harmonic phase[16, 17]. The first one estimated the relative phase by measuring harmonic yields in a mixed gas cell with H₂ and D₂[17]. The second study relied on the determination of group delay by using the reconstruction of attosecond beating by interference of two-photon transitions (RABITT) [16]. However, ensuring that both isotopes are at the same number density in the HHG interaction region is very challenging. Additionally, a sign ambiguity[17] and large experimental uncertainties of these measurements made it difficult to determine the phase difference and absolute delays accurately[16, 17]. Neither of those studies actually reproduced the measured phase difference theoretically.

Refer to caption — Figure 1: (a) Three step model of HHG with the potential energy surfaces of H₂ molecule. Tunneling ionization launches an electron wavepacket together with a nuclear vibrational wavepacket. The ionized electron is driven by the laser electric field and it returns to the parent ion at a certain recollision time. Simultaneously, the nuclear wavepacket evolves on the potential energy surface of the ground state of H ${}_{2}^{+}$ . Upon electron recombination a photon is emitted with an amplitude and phase dependent on the correlated electron-nuclear dynamics. (b) Schematics of the Gouy phase interferometric technique. Two gas jets are placed sequentially with a separation $\Delta z$ in a single laser focus. Two coherent high harmonic pulses are generated simultaneously from the top (H₂) and bottom (D₂) gas jets. A phase difference between the two pulses includes contributions from the Gouy phase shift of the driving laser field at the two jets positions and the intrinsic phase difference between the two species. By measuring HHG yields for different gas configurations and jet separations these two contributions can be disentangled.

To retrieve the phase information, the high-harmonic interferometers either split the driving laser beam into two paths and focus at different locations within a single gas jet that produce two phase-locked HHG photons [18, 19], or use a single driving laser beam in a gas mixture [2, 17]. The first method is limited to a single gas species, while the second has low resolution and issues with determining the sign of the relative harmonic phase. Recently, the phase difference between two atomic species was measured by an all-optical attosecond interferometric technique[20], where the delay between two extreme-ultraviolet (XUV) pulses is controlled by a two-segment mirror which can provide the temporal resolution of $\sim$ 6 attoseconds. Such an interferometer must maintain sub-cycle stability of its path difference (displacement between the two mirror segments) for XUV frequencies for the duration of the experiment while this path difference is being scanned. This is an exceedingly difficult task which severely limits practical utility of the device.

Building an interferometer in the XUV region is quite challenging for two reasons: firstly, it is challenging to control the delay of the XUV pulses precisely between the two arms with sub-cycle precision; secondly, the highly reflective XUV optics is yet to be developed. Our passively stabilized Gouy phase interferometer [21, 10], on the other hand, is an all-optical direct XUV interferometric technique. It does not require calibration of gas pressures to ensure the same number densities. Additionally, it does not require any XUV optics. The technique provides an elegant way to generate two coherent high harmonic pulses without splitting the driving laser and XUV beams. The two mutually coherent XUV pulses are generated by exploiting the inherent properties (Gouy phase) of a single Gaussian focused laser beam. Its unprecedented resolution of $\sim$ 300 $\mu$ rad ( $\sim$ 100 zeptoseconds) is a result of any instability in the distance between the two arms (the gas jets in this case) of the interferometer is determined relative to the Rayleigh length, $z_{R}$ of the fundamental laser beam as opposed to the XUV wavelength in a conventional optical interferometer [10].

Here, we apply the technique to investigate the effect of nuclear dynamics on the electron motion in molecular hydrogen by precise measurement of high harmonic phase difference (and corresponding HHG phase delays) produced in H₂ and D₂. Since the ionization potentials of H₂ ( $I_{p}$ = 15.43 eV) and D₂ ( $I_{p}$ = 15.46 eV) are almost identical [6], the difference in the phase accumulated by electron in the continuum is considered to be negligible. However, due to the nuclear mass difference, the evolution of nuclear wave packet while electron propagates in the continuum and then recombines to the ground state may differ substantially. The harmonic intensity from the heavy isotope was shown to be higher compared to the lighter molecule as harmonic emission is sensitive to the nuclear motion [6] though for D₂ the ionization probability is less than for H₂ [22]. The aim of this work is to measure a small phase difference of HHG signals and to gain an insight into the correlated electron–nuclear dynamics for the two isotopes of molecular hydrogen.

The high harmonics are generated by a linearly polarized, 800 nm, 9 fs pulses with a peak intensity of $\sim 5\times 10^{14}$ W/cm². Our advanced Gouy phase interferometer consists of two spatially separated gas jets along the propagation direction in a single laser focus. One of the gas jets (bottom) is fixed in position at the center of laser focus while the second jet (top) can move along the laser propagation direction before the focus of the laser beam, as shown in Fig. 1(b) (see Methods for details). The Gouy phase modulates the carrier envelope phase (CEP) of a focused laser pulse which results in a phase shift of XUV pulses generated from the two gas jets given by

\begin{split}\Delta\phi&=q\Delta\phi_{Gouy}\\ &=-q\tan^{-1}\Big{(}\dfrac{\Delta z}{z_{R}}\Big{)},\end{split}

(1)

where $q$ is the harmonic order and $\Delta\phi_{Gouy}$ is the Gouy phase difference at two gas jets positions (separated by $\Delta z$ ). The diameter of the gas jets is 200 $\mu$ m and gas pressure in both jets is kept at 100 Torr even though optimal phase matching happens at a slightly higher pressure. The short interaction region and low gas density help to reduce the macroscopic phase matching effects and to minimize the reabsorption of XUV photons by the second gas jet. The gases are switched in the jets by electronically actuated microvalves. That allows faster measurement and thus minimizes the errors due to laser fluctuation over time and also compensates for a small difference of intensities (and associated phases) at the two jets positions.

The interference fringes of the high harmonics from H23 to H37 generated with H₂ in both jets obtained by varying the separation between the jets are shown in Fig. 2(a). To optimize the resolution for all the observable harmonic orders, the isotopic phase difference measurements were performed at jet separations of 0.63 mm and 1.27 mm where the fluctuation of laser intensities is only 0.2% and 0.8% respectively.

To extract the relative phase difference between the harmonics generated from the two gases, for each of the two jet separations measurements are performed for two configurations. In the first configuration, D₂ and H₂ are delivered from the top and bottom jets respectively. The gases are swapped in the second configuration.

The HHG spectra are recorded for two different jet separations and for two gas configurations to account for systematic uncertainties. For each configuration, three HHG spectra are recorded: for both gas jets activated and for individual jet activated. The measured spectra are presented in Fig. 3. The relative phase difference between H₂ and D₂ is calculated as

\Delta\phi_{H_{2}-D_{2}}={\rm sin}^{-1}\Bigg{[}\dfrac{I_{H_{2}}+I_{D_{2}}+2\sqrt{I_{H_{2}}I_{D_{2}}}}{4\sqrt{I_{H_{2}}I_{D_{2}}}}\Big{(}\dfrac{\Delta I_{N}}{sin(q\Delta\phi_{Gouy})}\Big{)}\Bigg{]},

(2)

where, $I_{H_{2}}$ and $I_{D_{2}}$ are the harmonic intensities from individual gas jets of H₂ and D₂ respectively, and $\Delta I_{N}$ is the difference in the normalized intensity obtained from the two gas configurations. The results are presented in Fig. 4(a). The left axis corresponds to the phase difference averaged over both separations and both gas configurations. It should be noted, that our method measures both the value and sign of the phase difference and the two different configurations of gases induce opposite signs. Their absolute values were averaged and the shown phases correspond to emission from H₂ being delayed relative to D₂ as our measurements indicate (see Methods for details). The corresponding phase delays of harmonics from H₂ relative to D₂ are depicted on the right axis of Fig. 4(a) and are determined by

\Delta t_{H_{2}-D_{2}}=\dfrac{\Delta\phi_{H_{2}-D_{2}}}{\omega_{q}},

(3)

where, $\omega_{q}$ is the angular frequency of the $q^{th}$ harmonic order. The harmonics from H₂ are found to be $\sim$ 3 attoseconds delayed in phase with respect to D₂ for all the observable harmonic orders.

The experimental measurements are supported by numerical solutions of the four-dimensional time-dependent Schrödinger equation (4D-TDSE) within the single active electron (SAE) approximation (see Methods for details). In this model, the electronic and nuclear motions are both confined to a plane containing the molecular vibration and rotation. The simulated results (red crosses) are presented in Fig. 4(a) and they agree well with the experimental data (blue stars). Additional simulations showed that $\Delta\phi_{H_{2}-D_{2}}$ is relatively insensitive to the laser intensity, so the laser focal-volume averaging was not included.

According to the well-known three-step model, the total HHG phase includes contributions from the ionization, propagation, and recombination processes. Thus, for molecules aligned along different angles $\theta_{r}$ towards the laser polarization direction, the dynamics of ionization and recombination should be different, resulting in the $\theta_{r}$ -dependent $\Delta\phi_{H_{2}-D_{2}}$ . Numerically, by excluding molecular rotations and restricting the nuclear motion to a single vibrational degree of freedom we reduced our 4D-TDSE model to a three-dimesional one. We used that 3D-TDSE model to obtain $\theta_{r}$ -dependent $\Delta\phi_{H_{2}-D_{2}}$ , presented in Fig. 4(b). Only the molecules aligned at small angles (0-30 degrees) from the laser polarisation direction exhibit significant dependence of $\Delta\phi_{H_{2}-D_{2}}$ on harmonic order (energy of recolliding electrons) which could be attributed to two-center interference in the recombination step. Due to the different nuclear masses, the molecular bond elongation between the ionization and recombination steps is different for H₂ and D₂. In the recombination step, the complex two-center destructive interference occurs at a specific internuclear distance $R$ when the recombining electron momentum $k$ meets the condition $kR$ cos( $\theta_{r})=\pi$ , and the HHG phase undergoes a sudden change at the corresponding frequency $\omega=k^{2}/2+I_{p}$ . It is the difference in two-center interference between H₂ and D₂ that is responsible for the angle and energy dependence of the curves in Fig. 4(b). A more detailed analysis based on the Lewenstein model and the two-center interference can be found in the supplementary information.

In summary, we have demonstrated a novel all-optical Gouy phase interferometric technique for measuring HHG phase difference between two atomic or molecular species. We used this technique to measure the HHG phase difference between H₂ and D₂ to be about 0.2 radian, corresponding to 3 attosecond phase delay for all the harmonics. We also simulated this phase difference by numerically solving the time-dependent Schrödinger equation. The theoretical results obtained at the highest level of simulations agree quite well with the experiment. The simulations reveal that $\Delta\phi_{H_{2}-D_{2}}$ depends on the molecular alignment, bond stretching, and the two-center interference. The phase difference of harmonics from molecular isotopes can be used as a sensitive probe of ultrafast correlated electron-nuclear dynamics in molecules.

Authors contribution: D. L, I. L. and R. S. conceived, planned and lead the project. M. M and D. L. set up the interferometer and the experiments. M. M. carried out the measurement, analyzed and interpreted experimental data, improved the analysis and interpretation in discussion with D. L., I. L and R. S.. M. M wrote up experimental parts of this manuscript. L. X., F. H and W. U conducted numerical simulations. L. X and F. H wrote up the theoretical sections. N. H and H. X. contributed to important discussion on results and theoretical analysis. All authors contributed in discussions on the results and prepared the manuscript for submission.

Acknowledgements: The experimental research was supported by the Australian Research Council (LP140100813 and DP190101145) and Griffith University. M. M. was supported by a Griffith University International Postgraduate Research (GUIPRS) and a Griffith University Postgraduate Research (GUPRS) Scholarships. Theoretical work was supported by National Key R&D Program of China (2018YFA0404802, No. 2018YFA0306303), Innovation Program of Shanghai Municipal Education Commission (2017-01-07-00-02-E00034), National Natural Science Foundation of China (NSFC) (Grant No. 11574205, No. 91850203)

References

[1] O. Smirnova, Y. Mairesse, S. Patchkovskii, N. Dudovich, D. Villeneuve, P. Corkum, and M. Y. Ivanov. High harmonic interferometry of multi-electron dynamics in molecules. Nature, 460(7258):972, 2009.
[2] B. K. McFarland, J. P. Farrell, P. H. Bucksbaum, and M. Gühr. High-order harmonic phase in molecular nitrogen. Phys. Rev. A, 80:033412, 2009.
[3] P. Peng, C. Marceau, and D. M. Villeneuve. Attosecond imaging of molecules using high harmonic spectroscopy. Nat. Rev. Phys., 1:144, Feb 2019.
[4] T. Kanai, S. Minemoto, and H. Sakai. Quantum interference during high-order harmonic generation from aligned molecules. Nature, 435(7041):470, 2005.
[5] Wen Li, X. Zhou, R. Lock, S. Patchkovskii, A. Stolow, H. C. Kapteyn, and M. M. Murnane. Time-Resolved Dynamics in N2O4 Probed Using High Harmonic Generation. Science, 322(5905):1207–1211, 2008.
[6] S. Baker, J. S. Robinson, C. A. Haworth, H. Teng, R. A. Smith, C. C. Chirilă, M. Lein, J. W. G. Tisch, and J. P. Marangos. Probing proton dynamics in molecules on an attosecond time scale. Science, 312(5772):424–427, 2006.
[7] S. Baker, J. S. Robinson, M. Lein, C. C. Chirilă, R. Torres, H. C. Bandulet, D. Comtois, J. C. Kieffer, D. M. Villeneuve, J. W. G. Tisch, and J. P. Marangos. Dynamic two-center interference in high-order harmonic generation from molecules with attosecond nuclear motion. Phys. Rev. Lett., 101:053901, 2008.
[8] A. Zaïr, T. Siegel, S. Sukiasyan, F. Risoud, L. Brugnera, C. Hutchison, Z. Diveki, T. Auguste, J. W. G. Tisch, Pascal. Saliéres, M. Y. Ivanov, and J. P. Marangos. Molecular internal dynamics studied by quantum path interferences in high order harmonic generation. Chem. Phys., 414:184 – 191, 2013.
[9] J. Itatani, J. Levesque, D. Zeidler, H. Niikura, H. Pépin, J. C. Kieffer, Paul B. Corkum, and D. M. Villeneuve. Tomographic imaging of molecular orbitals. Nature, 432(7019):867, 2004.
[10] M. H. Mustary, D. E. Laban, J. B. O. Wood, A. J. Palmer, J. Holdsworth, I. V. Litvinyuk, and R. T. Sang. Advanced gouy phase high harmonics interferometer. Journal of Physics B: Atomic, Molecular and Optical Physics, 51(9):094006, apr 2018.
[11] Shubhadeep Biswas, Benjamin Förg, Lisa Ortmann, Johannes Schötz, Wolfgang Schweinberger, Tomáš Zimmermann, Liangwen Pi, Denitsa Baykusheva, Hafiz A Masood, Ioannis Liontos, et al. Probing molecular environment through photoemission delays. Nature Physics, pages 1–6, 2020.
[12] M. Lein. Attosecond probing of vibrational dynamics with high-harmonic generation. Phys. Rev. Lett., 94:053004, 2005.
[13] Xue-Bin Bian and André D. Bandrauk. Probing nuclear motion by frequency modulation of molecular high-order harmonic generation. Phys. Rev. Lett., 113:193901, Nov 2014.
[14] Lixin He, Qingbin Zhang, Pengfei Lan, Wei Cao, Xiaosong Zhu, Chunyang Zhai, Feng Wang, Wenjing Shi, Muzi Li, Xue-Bin Bian, et al. Monitoring ultrafast vibrational dynamics of isotopic molecules with frequency modulation of high-order harmonics. Nature communications, 9(1):1–7, 2018.
[15] P. Lan, M. Ruhmann, L. He, C. Zhai, F. Wang, X. Zhu, Q. Zhang, Y. Zhou, M. Li, M. Lein, and P. Lu. Attosecond probing of nuclear dynamics with trajectory-resolved high-harmonic spectroscopy. Phys. Rev. Lett., 119:033201, 2017.
[16] S. Haessler, W. Boutu, M. Stankiewicz, L. J. Frasinski, S. Weber, J. Caillat, R. TaÃ¯eb, A. Maquet, P. Breger, P. Monchicourt, B. Carré, and P. Salières. Attosecond chirp-encoded dynamics of light nuclei. Journal of Physics B: Atomic, Molecular and Optical Physics, 42(13):134002, jun 2009.
[17] T. Kanai, E. J. Takahashi, Y. Nabekawa, and K. Midorikawa. Observing the attosecond dynamics of nuclear wavepackets in molecules by using high harmonic generation in mixed gases. New Journal of Physics, 10(2):025036, feb 2008.
[18] M. Bellini, C. Lyngå, A. Tozzi, M. B. Gaarde, T. W. Hänsch, A. L’Huillier, and C.-G. Wahlström. Temporal coherence of ultrashort high-order harmonic pulses. Phys. Rev. Lett., 81:297–300, 1998.
[19] C. Lyngå, M. B. Gaarde, C. Delfin, M. Bellini, T. W. Hänsch, A. L’ Huillier, and C.-G. Wahlstrom. Temporal coherence of high-order harmonics. Physical Review A, 60:4823–4830, 1999.
[20] Doron Azoury, Omer Kneller, Shaked Rozen, Barry D Bruner, Alex Clergerie, Yann Mairesse, Baptiste Fabre, Bernard Pons, Nirit Dudovich, and Michael Krüger. Electronic wavefunctions probed by all-optical attosecond interferometry. Nature Photonics, 13(1):54–59, 2019.
[21] D. E. Laban, A. J. Palmer, W. C. Wallace, N. S. Gaffney, R. P. M. J. W. Notermans, T. T. J. Clevis, M. G. Pullen, D. Jiang, H. M. Quiney, I. V. Litvinyuk, D. Kielpinski, and R. T. Sang. Extreme ultraviolet interferometer using high-order harmonic generation from successive sources. Phys. Rev. Lett., 109:263902, 2012.
[22] X. Wang, H. Xu, A. Atia-Tul-Noor, B. T. Hu, D. Kielpinski, R. T. Sang, and I. V. Litvinyuk. Isotope effect in tunneling ionization of neutral hydrogen molecules. Phys. Rev. Lett., 117:083003, Aug 2016.
[23] F. He, C. Ruiz, and A. Becker. Absorbing boundaries in numerical solutions of the time-dependent schrödinger equation on a grid using exterior complex scaling. Phys. Rev. A, 75:053407, May 2007.
[24] Kosloff R. and Tal-Ezer H. A direct relaxation method for calculating eigenfunctions and eigenvalues of the schrodinger equation on a grid. Chemical Physics Letters, 127:223, June 1986.
[25] Crank J. and Nicolson P. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Math. Proc. Cambridge Philos. Soc., 43:50, January 1947.
[26] M. Lein, N. Hay, R. Velotta, J. P. Marangos, and P. L. Knight. Role of the intramolecular phase in high-harmonic generation. Phys. Rev. Lett., 88:183903, Apr 2002.

Methods

Experiment: The experiment is performed with Quantronix Ti:Sapphire laser system consisting of Ti-light modelocked oscillator and Odin-II chirp pulsed amplifier delivering 3 mJ, 37 fs pulses centred at 800 nm wavelength at 1 kHz repetition rate. The laser pulse is coupled into Neon gas filled hollow core fiber that spectrally broadens the pulse by self-phase modulation. A pair of dielectric chirp compensation mirrors are then used to introduce negative GDD to generate 9 fs pulses. The laser beam is then focused into the HHG chamber by a 750 mm focal length spherical mirror, which leads to an estimated intensity of $\sim$ 5.0 $\times$ 10¹⁴ W/cm² in the HHG interaction region.

The Gouy phase interferometer consists of two spatially separated gas jets; fixed bottom jet pointing up towards the laser focus and moveable top jet pointing down before the focus that can translate along all three axes. Each jet is connected with two electromagnetically actuated pulsed microvalves, so that the gases in the jet can be switched according to the measurement requirement. H₂ and D₂ gases are supplied from separate cylinders connected to the pulsed microvalves.

Synchronization of microvalve opening delay with the laser pulse reduces the background pressure in the HHG generation and detection chambers, which facilitates better signal to noise ratio (for details, see [10]). The switching of gases in the jets allows faster measurement and thus minimizing the effect due to laser fluctuation over time and also compensates for a small difference of intensities (and associated phases) at the two different positions along the laser focus. The inner diameter of the jet is only 200 $\mu m$ , which is very thin compared to the Rayleigh length of the driving laser beam. The schematic of full experimental setup is shown in the supplementary Fig 1. More details about the instrumentation can be found in [10].

Derivation of phase shift: The electric field of the HHG signal generated from H₂ gas in top jets in figure 1 for the $q^{th}$ harmonic order is

\vec{E}_{H_{2}}=E_{H_{2}}e^{i(q\omega t+q\phi_{Gouy1}+\phi_{H_{2}})},

(1)

here $E_{H_{2}}$ is maximum field amplitude, $q$ is the harmonic order, $\omega$ is the angular frequency of the driving laser, $\phi_{Gouy1}$ and $\phi_{H_{2}}$ are the Gouy phase at the top jet position and intrinsic phase of H₂ molecule, respectively. Similarly, the electric field of the HHG signal generated from bottom jet with D₂ gas is

\vec{E}_{D_{2}}=E_{D_{2}}e^{i(q\omega t+q\phi_{Gouy2}+\phi_{D_{2}})},

(2)

where $E_{D_{2}}$ is maximum field amplitude, $\phi_{Gouy2}$ and $\phi_{D_{2}}$ are the Gouy phase at bottom jet position and intrinsic phase of D₂ molecules, respectively. If both jets are ON, two HHG pulses will be generated from top ( $H_{2}$ ) and bottom ( $D_{2}$ ) jet with single laser pulse and the total field will be the coherent sum of these two electric fields of equation 1 and 2

\vec{E}_{H_{2}D_{2}}=\vec{E}_{H_{2}}+\vec{E}_{D_{2}}.

(3)

In terms of intensity, the irradiance of the sum of these electric fields is

\begin{split}I_{H_{2}D_{2}}&=I_{H_{2}}+I_{D_{2}}+2\sqrt{I_{H_{2}}I_{D_{2}}}cos(q\Delta\phi_{Gouy}+\Delta\phi_{H_{2}-D_{2}}),\end{split}

(4)

where $\Delta\phi_{Gouy}=\phi_{Gouy1}-\phi_{Gouy2}$ and $\Delta\phi_{H_{2}-D_{2}}=\phi_{H_{2}}-\phi_{D_{2}}$ . Data is then normalized to the in-phase sum of the HHG signals from the H₂ jet and D₂ jet

N_{H_{2}D_{2}}=\dfrac{I_{H_{2}}+I_{D_{2}}+2\sqrt{I_{H_{2}}I_{D_{2}}}cos(q\Delta\phi_{Gouy}+\Delta\phi_{H_{2}-D_{2}})}{I_{H_{2}}+I_{D_{2}}+2\sqrt{I_{H_{2}}I_{D_{2}}}}.

(5)

When the ordering of gases between top and bottom jets is swapped, i.e. D₂ is in the top gas jet and H₂ in the bottom gas jet, then the normalized intensity will be

N_{D_{2}H_{2}}=\dfrac{I_{H_{2}}+I_{D_{2}}+2\sqrt{I_{H_{2}}I_{D_{2}}}cos(q\Delta\phi_{Gouy}+\Delta\phi_{D_{2}-H_{2}})}{I_{H_{2}}+I_{D_{2}}+2\sqrt{I_{H_{2}}I_{D_{2}}}}.

(6)

The sign of the molecular phases, $\Delta\phi_{H_{2}-D_{2}}$ and $\Delta\phi_{D_{2}-H_{2}}$ will be opposite in this case,i.e. $\Delta\phi_{H_{2}-D_{2}}=-\Delta\phi_{D_{2}-H_{2}}$ . Finally, by subtracting eq(5) from eq(6), the relative high harmonic phase difference between H₂ and D₂ is expressed as

\Delta\phi_{H_{2}-D_{2}}=sin^{-1}\Bigg{[}\dfrac{I_{H_{2}}+I_{D_{2}}+2\sqrt{I_{H_{2}}I_{D_{2}}}}{4\sqrt{I_{H_{2}}I_{D_{2}}}}\Big{(}\dfrac{\Delta I_{N}}{sin(q\Delta\phi_{Gouy})}\Big{)}\Bigg{]},

(7)

where $I_{H_{2}}$ and $I_{D_{2}}$ are HHG yield with H₂ and D₂ gas respectively. $\Delta I_{N}=N_{D_{2}H_{2}}-N_{H_{2}D_{2}}$ and $\Delta\phi_{Gouy}$ is Gouy phase difference of the driving laser.

Data analysis: The acquired data is processed and analyzed by using MATLAB. For each spectra shown in Fig. 3(a) in the main text, the data points around each harmonic peak are plotted. The peak intensity is determined by fitting these data points by a Gaussian function using the least square fitting and the fitting function is

f(x)=a\exp^{\Bigg{(}\dfrac{-(x-b)^{2}}{2c^{2}}\Bigg{)}}

(8)

where $a,b$ and $c$ are the fitting parameters. The peak amplitude of the harmonic signal is given by the fitting parameter $a$ . The normalized intensities, as described by equations 5 and 6, are calculated from these peaks.

The fitting procedure estimates $1\sigma$ $(68\%)$ confidence bound to quantify the uncertainties, which are used to calculate the standard error. Since $N_{H_{2}D_{2}}$ is a function of $I_{H_{2}},I_{D_{2}},I_{H_{2}D_{2}}$ and corresponding error is the combined uncertainty of these $\delta N_{H_{2}D_{2}}(\delta I_{H_{2}},\delta I_{D_{2}},\delta I_{H_{2}D_{2}})$ , therefore the error bar of $\Delta N_{{D_{2}-H_{2}}}$ is obtained by the propagation of errors and is defined as $\delta(\Delta N_{{D_{2}}-{H_{2}}})=\sqrt{(\delta N_{H_{2}D_{2}}^{2}+\delta N_{D_{2}H_{2}}^{2})}$ .

The effective $\Delta\phi_{Gouy}$ for each harmonic order is obtained from the normalized intensity vs separation data of Fig. 2. These data are fitted with a decaying cosine function as

f(x)=W\exp^{-Yx}cos(Zx)+C

(9)

where $W$ , $Y$ , $Z$ and $C$ are the fitting parameters and the $Zx$ gives the $q\Delta\phi_{Gouy}$ . Finally, the total error bar of phase shift between H₂ and D₂ is the propagation error, defined as $\delta\phi_{{D_{2}}-{H_{2}}}=\sqrt{\delta_{\Delta N_{{D_{2}}-{H_{2}}}}^{2}+\delta_{sin(q\Delta\phi_{Gouy})}^{2}}$ .

TDSE simulations: We numerically simulate the four dimensional time-dependent Schrödinger equation (4D-TDSE) with the single-active-electron (SAE) approximation for H₂ and D₂ (Hartree atomic units are used $e=m=\hbar=1$ unless stated otherwise),

i\frac{\partial}{\partial t}\Psi(\textbf{{r}},\textbf{{R}};t)=\left[T+V(\textbf{{r}},\textbf{{R}})\right]\Psi(\textbf{{r}},\textbf{{R}};t),

(10)

where r = (x, y) and R = ( $\rho$ , $\theta$ ) are the electronic coordinate and the internuclear displacement, respectively. The total kinetic energy operator is given as in the velocity gauge

T=\frac{[\textbf{{p}}_{\textbf{{r}}}+\textbf{{A}}(t)]^{2}}{2}+\frac{\textbf{{p}}_{\textbf{{R}}}^{2}}{2\mu}.

(11)

Here, $\mu$ is the reduced nuclear mass, and the laser vector potential is $\textbf{A}(t)=-\int_{0}^{t}\textbf{E}(t^{\prime})dt^{\prime}$ with $\textbf{E}(t^{\prime})$ being the electric field vector. $\textbf{{p}}_{\textbf{{r}}}$ and $\textbf{{p}}_{\textbf{{R}}}$ are the electronic canonical momentum operator and the relative nuclear momentum operator, respectively.

This model reproduces molecular potential energy curves [12] and thus it is expected to capture the main features of nuclear dynamics of H₂ and D₂ in laser fields. The interaction between the nuclei and the laser field vanishes because of the homonuclear diatomic molecule. The effective potential energy is expressed as

V(\textbf{{r}},\textbf{{R}})=V_{BO}^{+}(R)-\sum_{j=1,2}\frac{Z(R,|\textbf{{r}}-\textbf{{R}}_{j}|)}{\sqrt{|\textbf{{r}}-\textbf{{R}}_{j}|^{2}+0.5}},

(12)

which describes the interaction between the outer active electron and the atomic core. $V_{BO}^{+}(R)$ is the lowest isotropic Born-Oppenheimer approximation potential energy surface 1 $s\sigma_{g}$ of H ${}_{2}^{+}$ and D ${}_{2}^{+}$ . R = $|\textbf{{R}}|$ is the internuclear distance and $\textbf{{R}}_{j}=\textbf{{R}}/2$ is the $j$ -th nuclear displacement in the center of molecular mass coordinate. $Z(R,|\textbf{{r}}-\textbf{{R}}_{j}|)=\{1+\textrm{exp}[-(\textbf{{r}}-\textbf{{R}}_{j})^{2}/\sigma^{2}(R)]\}$ is introduced as an effective nuclear charge, which includes the average screening of the nuclei by the second electron with the screening parameter $\sigma(R)$ . By choosing proper $\sigma(R)$ , the resulting potential energy surface based on Eq. 12 can match the real electronic ground energy of H₂ and D₂. The other details about this potential can been found in [12].

In our simualtions, to suppress the unphysical reflection from boundaries and eliminate the long trajectory contributions to high harmonics, a cos^1/8 masking function has been adopted [23]. The initial vibration-rotation ground state is obtained using the imaginary time propagation [24], and the Crank-Nicolson method is used to propagate the molecular wave packet in real time [25]. The harmonic spectrum is calculated via the Fourier transformed dipole acceleration

F_{x}(\omega)=|\frac{1}{(2\pi)^{3/2}}\int D_{x}(t)e^{-i\omega t}dt|^{2}

(13)

with the dipole acceleration

D_{x}(t)=\langle\Psi(\textbf{{r}},\textbf{{R}};t)\left|-\frac{\partial V(\textbf{{r}},\textbf{{R}})}{\partial x}\right|\Psi(\textbf{{r}},\textbf{{R}};t)\rangle.

(14)

Here we extract the phase of the complex integral in Eq. 13 as the phase of the emitted harmonic radiation [26].

In addition to the 4D calculations, we also performed the reduced 3D calculations by omitting the dynamical terms related to $\theta$ in the 4D model. This can be done because the molecular rotation can be neglected during the harmonic generation. Such a 3D model may directly reveal the contribution from the molecule aligned at a certain $\theta_{r}$ based on the $\theta_{r}$ -parameterized molecular wave packet. In all 3D calculations, all parameters and procedures are the same as those implemented in 4D calculations, which can be found in the supplementary information.

Supplementary information

The calculation parameters:

Here, the used laser field polarized along the x-axis is written as

E(t)=E_{0}\textrm{sin}^{2}(\pi t/\tau)\textrm{cos}(\omega t+\phi).

(1)

$E_{0}$ , $\tau$ , $\omega$ , and $\phi$ denote the amplitude, the total duration, the angular frequency, and the carrier envelope phase (CEP) of the laser pulse, respectively. In calculations, the laser pulse has the wavelength 800 nm ( $\omega=0.057$ a.u.) and lasts ten optical cycles ( $\tau=10T$ ) with zero CEP ( $\phi=0$ ), and its peak intensity is $5\times 10^{14}$ W/cm² ( $E_{0}=0.12$ a.u.). The spatial steps in simulations are $\Delta x=\Delta y=0.3$ a.u., $\Delta\rho=0.04$ a.u., and $\Delta\theta=2\pi/300$ and the time step is $\Delta t=0.05$ a.u.. The 4D simulation box has the corresponding grids points: 250 $\times$ 100 $\times$ 100 $\times$ 300. In 3D calculations, the space grids points are 250 $\times$ 100 $\times$ 100. We obtain the HHG spectra for H₂ and D₂, as shown in Fig. 1. Benefited from the better vibrational autocorrelation character, D₂ can generate more intense harmonics than H₂ and their intensity difference increases with the harmonic order [12]. The convergence of the calculations is tested by using finer spatial and time grids whereby almost identical results are obtained.

[Uncaptioned image] — List of suppfigures 1 The calculated HHG spectra for H₂ (black line) and D₂ (red line) from 4D-TDSE simulation.The laser parameters are the wavelength $\lambda$ = 800 nm, the intensity I = $5\times 10^{14}$ W/cm², the total duration $\tau$ = 10T, and the sin² envelope.

Further analyses based on the Lewenstein model : Though the TDSE simulation may quantitatively reproduce the experimental measurements, the Lewenstein model is able to analysis the problem transparently. The time-dependent dipole moment is expressed as

\begin{split}D(t)=2i\int_{0}^{t}dt^{\prime}E(t^{\prime})C(t-t^{\prime})\int\frac{d^{3}p}{(2\pi)^{3}}\bar{d}^{*}[p+A(t)]\\ e^{-i\int_{t^{\prime}}^{t}\{[p+A(t^{\prime\prime})]^{2}/2+I_{p}\}dt^{\prime\prime}}\bar{d}[p+A(t^{\prime})]+c.c.,\end{split}

(2)

where $p$ and $I_{p}$ are the electronic canonical momentum and molecular ionization energy, respectively. $t$ and $t^{\prime}$ denote the electronic ionization and recombination moments. As described by Eq. 2, the nuclear contribution to HHG is related to the vibrational autocorrelation function $C(t)=\int\chi(R,t=0)\chi(R,t)dR$ and the other terms represent the electronic contribution. $\chi(R,t)$ is approximately governed by the following 1D-TDSE

i\frac{\partial}{\partial t}\chi(R,t)=\left[\frac{p_{R}^{2}}{2\mu}+V(R)\right]\chi(R,t),

(3)

where $V(R)$ is the potential curve of $1s\sigma_{g}$ state for H ${}_{2}^{+}$ or D ${}_{2}^{+}$ and R is the internuclear distance.

Based on Eq. 2, the phase difference must be contributed by the electron dynamics and the nuclear relaxation after the molecular single ionization, as well as the electron-nuclei coupling. We first discuss the the simplest case, that is, $\theta_{r}=90^{\rm o}$ , where the nuclear relaxation and the electronic motion are mutually perpendicular. The phase difference $\Delta\phi_{H_{2}-D_{2}}$ is mainly contributed by the nuclear part $\Delta\phi_{n}$ as well as electronic ionization and recombination since the electron dynamical parts in the continuum are almost the same for H₂ and D₂. To pick up the phase differences induced by the nuclear relaxation, we solely propagate the nuclear wave packets of H ${}_{2}^{+}$ and D ${}_{2}^{+}$ for $\tau_{q}$ by using Eq. 3. Here the time interval $\tau_{q}$ between ionization and recombination is estimated from the simple-man model. The initial state is the ground nuclear state of H₂ or D₂ in simulations. The spatial and time steps are $\Delta R=0.02$ a.u. and $\Delta t=0.01$ a.u., respectively. The numerical method is the same as those used in 3D or 4D simulations in Methods. The obtained $\Delta\phi_{n}$ shown in Fig. 2(a) is very small about -0.01 radian. The gap between $\Delta\phi_{n}$ and 3D-TDSE results at 90 degrees is up to about 0.2 radian. We inferred that this gap mainly roots in ionization and recombination processes, which is missed in the extended Lewenstein model.

Next we turn to other angle cases. Just like Young’s double-slit interference in wave optics, when the returning electron momentum $k$ meets the condition $kR$ cos $(\theta_{r})=\pi$ , the high harmonic shows minimal radiation and its phase undergoes an abrupt change at the corresponding frequency $\omega=k^{2}/2+I_{p}$ . For example, at $\theta_{r}=0^{\rm o}$ case, the destructive interference mainly changes the phase differences and radiation intensities around the 35th harmonic order, as presented in Fig. 2(a) and (b) respectively, which is estimated as follows. For the 35th order harmonic, the time interval between ionization and rescattering is about 46 a.u. according to the simple-man model under our laser parameters. The nuclear wave packet distributions of H₂ and D₂ propagated for $\tau_{35}$ by Eq. 3 have the most conspicuous difference at the nuclear distance $R_{c}\approx 1.8$ a.u.. The corresponding momentum $k$ and nuclear distance $R_{c}$ at the rescattering moment right meet the destructive interference condition. Thus, the phase difference at 35th order has the maximal modulation depth with about 0.4 radian compared with $\theta_{r}=90^{\rm o}$ case. However, the nuclear wave packets of H₂ and D₂ distribute within a wide $R$ range, so the double-slit interference has also slight effect on the neighbouring order harmonics. Similarly, for $\theta_{r}=30^{\rm o}$ case, the destructive interference mainly changes the phase differences around 45th harmonic order, as shown in Fig. 2(a), which is beyond our measurement range.