
\documentclass[fleqn,25pt]{wlscirep}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{multirow}

\newcommand{\Rmnum}[1]{\expandafter\@slowromancap\romannumeral #1@}
\newcommand\subtitle[1]{{\small #1}}

\renewcommand\thesection{\Roman{section}}
\usepackage[version=3]{mhchem}




\usepackage{caption}
\captionsetup[figure]{name={Supplementary Figure }}


\title{Supplementary Information for "Directly imaging excited states-resolved coherent nuclear motions of water with picometer-femtosecond precision" }


\author[1,+]{Zhenzhen Wang}
\author[2,+]{Xiaoqing Hu}
\author[1]{Shengpeng Zhou}
\author[1]{Xitao Yu}
\author[1]{Yizhang Yang}
\author[1]{Banchi Zhao}
\author[2]{Zheng Shu}
\author[2]{Zhenpeng Wang}
\author[1]{Xiaokai Li}
\author[2,3*]{Yong Wu}
\author[2]{Jianguo Wang}
\author[1*]{Chuncheng Wang}
\author[1*]{Dajun Ding}


\affil[1]{Institute of Atomic and Molecular Physics and Jilin Provincial Key Laboratory of Applied Atomic and Molecular Spectroscopy, Jilin University, Changchun 130012, China}
\affil[2]{Key Laboratory of Computational Physics, Institute of Applied Physics and Computational Matematics, Beijing 100088, China}
\affil[3]{HEDPS, Center of Applied Physics and Technology, Peking University, 100871 Beijing, China.}

\affil[*]{Correspondence and requests for materials should be addressed to Chuncheng Wang(email: ccwang@jlu.edu.cn) or to Dajun Ding (email: dajund@jlu.edu.cn) or to Yong Wu (email: wuyong@jlu.edu.cn)}

\affil[+]{these authors contributed equally to this work}


\begin{document}

	
	\flushbottom
	\maketitle	
	\thispagestyle{empty}


\maketitle
\tableofcontents
\newpage




\section{Kinetic energy correlation and Dalitz diagram }

\begin{figure}[h]
	\centering
	\includegraphics[width=0.9\linewidth]{figS1.ai}
	\caption{(Color online)  Kinetic energy correlation of two D$\rm ^+$ and Dalitz diagram. (a) and (b) present energy correlation maps of the two D$\rm ^+$ for linear and circular polarisation, respectively, and the distributions with strong polarisation-dependence are labelled \uppercase\expandafter{\romannumeral1} and \uppercase\expandafter{\romannumeral2},  which originate from TERCE.  (c) simulated and (d) experimental Dalitz diagrams for events from regions \uppercase\expandafter{\romannumeral1} and \uppercase\expandafter{\romannumeral2}.}
	\label{fig1}
\end{figure}

The energy correlations between two D$^+$ are presented in Supplementary Fig. \ref{fig1} (a) and (b) for linearly and circularly polarised light,  respectively.  The absence of regions \uppercase\expandafter{\romannumeral1} and \uppercase\expandafter{\romannumeral2} in the case of circularly polarised light suggests that the three-body electron recollision-assisted Coulomb explosion (TERCE) is the dominant underlying mechanism for their formation.
Region \uppercase\expandafter{\romannumeral1} is concentrated along the diagonal of the energy correlation map,  which means that the two D$^+$ have the same energy, and their geometry before CE would thus be expected to be symmetrical.  However,  region \uppercase\expandafter{\romannumeral2} presents a distribution along the inverse diagonal,  which suggests that the two D$^+$ are strongly correlated but with different energies;  thus,  we can infer that the geometry before CE is asymmetrical.
Region \uppercase\expandafter{\romannumeral1} can be roughly separated from the strong field enhanced ionisation by selecting events with KER values higher than 25 eV\cite{guillemin_selecting_2015}.  In region \uppercase\expandafter{\romannumeral2},  the events exhibit larger energy differences between the two D$^+$, which are also distinguishable from other events.  We present the Dalitz diagram of the events from regions \uppercase\expandafter{\romannumeral 1} and \uppercase\expandafter{\romannumeral2} in Supplementary Fig. \ref{fig1} (d) to show the distinct differences between them.  The X- and Y-axes of the Dalitz diagram are defined as
%\begin{gather*}

\begin{center}
	\begin{equation}
	X  =  (\epsilon_{D^+_{1}}- \epsilon_{D^+_{2}})/(\sqrt{3}{\epsilon_K}),    ~~~~~  Y  =   \epsilon_{O^+}  / {\epsilon_K}  -  1/3
	\end{equation}
\end{center}
%%

where $\rm {\epsilon_K}$ denotes the total kinetic energy release\cite{neumann_fragmentation_2010}.  As indicated by the solid elliptical curves in the simulated Dalitz diagram (Supplementary Fig. \ref{fig1}(c)), the central distribution (region \uppercase\expandafter{\romannumeral1}) centred at zero (X-axis) represents the symmetrical concerted CE,  whereas the arm-like distributions located at 0.2 and -0.2 (region \uppercase\expandafter{\romannumeral2})indicate that those events originate from concerted CE processes with asymmetrical geometry.


\section{Two-dimensional fitting procedure }

A two-dimensional (2D) multiple Gaussian peak fitting program was developed to fit the experimental data.  The fitting function can be expressed as
\begin{center}
	\begin{equation}
	f(x,y)   =    a_{0}+$$\sum_{i=1}^n$$ a_{i} exp\left(- \left(\frac{(x-c_{x,i})^2}{2\sigma_{x,i}^2} + \frac{(y-c_{y,i})^2}{2\sigma_{y,i}^2}\right)\right)
	\end{equation}
\end{center}





Here, n is the assumed number of components, c and $\rm \sigma$ denote the peak position and variance of the 2D Gaussian distributions, respectively.  According to the principle of nonlinear least squares (NLS),  several sets of parameters,  which correspond to the local minima for the given initial values of the parameters,  can satisfy the NLS.  To determine the global minimum,  we implemented the following improvements. \\
1.	The parameters were assigned suitable initial values. \\
2.	A set of parameters was obtained according to NLS.  Its fitting error was recorded as the initial error of the global minimum,  and these parameters were recorded as the best parameters. \\
3. Some of the peak positions in the set of best parameters were reset to random values within the allowed range.  These values were designated to be the initial values of the parameters for the next cycle of NLS.  Then,  the next NLS computation was performed. \\
4. A new set of parameters was obtained from the new NLS cycle,  and a new fitting error was also obtained.  If the new fitting error was smaller than the error in step 2,  the present fitting error and the corresponding parameters were assumed to be the new value of the error and the new best parameters for the global minimum. Then,  the program proceeded to Step 3 for the next iterative cycle. The loop is terminated if the error from the new loop is less than the preset global minimum error or if the number of cycles exceeds the predefined number.\\


\begin{figure}[h]
	\centering
	\includegraphics[width=0.9\linewidth]{figS2.ai}
	\caption{2D global fit of experimental results with six components, and the five components in region \uppercase\expandafter{\romannumeral1} marked with circles.  (b) Summary of the residual errors of the 2D fit with an increasing number of components used in the fitting procedure.}
	\label{fig2}
\end{figure}

Supplementary Fig. \ref{fig2} (a) presents the outputs of our 2D fit for Fig. 1(c)  in the main text.  In the present case, six components were ultimately shown to fit the experimental results very well, and this number was further verified by estimating the residual errors, as shown in Supplementary Fig. \ref{fig2} (b).  The residual error $\rm E_{res.}$ is defined as:
\begin{center}
	\begin{equation}
	E_{res.} =   $$\frac{1}{n}\sum_{i=1}^n (Y_{i}-\overset{\wedge}Y_{i})^2$$
	\end{equation}
\end{center}


where $\rm Y_{i}$ and $\rm  \overset{\wedge}Y_{i}$ represent the true and predicted values, respectively.
We found that fitting with five components significantly increased the residual errors, whereas increasing the number of components to more than six did not further decrease the residual errors; thus, the global optimum can be reached with six components, which is an appropriate choice in the present case. Five dominant components are discussed in the main text because they can be roughly distinguished among the raw experimental data. In contrast, the sixth component (marked by the grey dashed circle) is weak and cannot be easily resolved in the raw data (Fig. 1(a) in the main text); thus, we refrained from discussing this component in the present work to avoid over-interpretation.




\section{Simulation details }


\begin{figure}[h]
	\centering
	\includegraphics[width=0.9\linewidth]{figS3.ai}
	\caption{(Color online) Calculated strong-field ionisation probability of cation. This probability was calculated by MO-ADK for different geometries, where $\theta_{DOD}$ ranges from 100$^\circ$ to 180$^\circ$ and $\rm R_{OD}$ is set as 100 pm. (a) Calculated time-dependent ionisation probability of $\rm D_2O^{+}$ within the laser pulse duration; (b) magnification of the ionisation probability in the central region of the laser pulse ; (c) total ionisation probability as a function of $\rm \theta_{DOD}$, normalized to 180$^\circ$ }
	\label{fig3}
\end{figure}

In the present simulation, the electric field E(t) is equal to  $\rm E_0sin^2(\frac{\pi t}{\tau})sin( \omega t+\varphi)$, where E$\rm _{0}$ is the peak amplitude of the electric field, $\rm \omega$ is the carrier frequency, $\rm \tau$ is the pulse duration, and $\rm \varphi$ is the carrier envelope phase, which was set as $\rm \frac{\pi}{2}$ in this study. To simplify the calculation, the instants of single ionisation and double ionisation induced by strong-field tunnelling, t$\rm _1$ and t$\rm _2$, are only selected from t = n T, where T is $\frac{8}{3}$ fs and n = 12, 12.5, 13,$\cdots$, 21, $\cdots$, 28.5, 29. Furthermore, the time interval $\rm t{_2 }$ $-$ $\rm t{_1}$ between single and double ionisation was counted from 0 to 10 T. The calculated strong-field ionisation probability with MO-ADK is shown in Supplementary Fig. \ref{fig3}(a), and the probability from the most intense region of the laser pulse are shown in Supplementary Fig. \ref{fig3} (b). As we can see, because of the lower ionisation potential in the large $\rm \theta_{DOD}$ region, the ionisation probability at 180$^\circ$ is five times larger than that at 100$^\circ$, which can be seen in the total ionisation probability (Supplementary Fig. \ref{fig3}(c)). The simulations suggest that geometries with larger bond angles have the most significant influence on the projected wave-packets of D$\rm _2$O$\rm ^{2+}$.

The simulated time-dependent vibrational wave-packets in the X and A states at different time intervals are shown in Supplementary Figs. \ref{fig4} and \ref{fig5}. The initial molecular wave-packet is set as the ground vibrational state of the neutral molecule and is then projected to the eigen vibrational states in the A and X states. The vibrational wave-packet before secondary ionisation is readily obtained by multiplying by the evolution factor $\rm exp(-iE_n(t_2 - t_1))$. The projected time-dependent wave-packets in the dication are shown in the second row, where the wave-packets evolutions during electron recollision were considered. The time-dependent wave-packets in the quartet and doublet states of the trication are also shown in the third and fourth rows, respectively. The calculated average distributions of $\rm R_{OD}$ and $\rm \theta_{DOD}$ are shown in Supplementary Fig. \ref{fig6}. The amplitude of the oscillations in the A state is larger than that in the X state. The oscillation period is approximately 20 fs for both the bond length and bond angle for the A state. This period is approximately 14 fs for the bond length and 32 fs for the bond angle for the X state. Because of the fast large-amplitude oscillation of the wave-packets in the A state, the ionisation probability of double ionisation significantly affects the final wave-packet distributions before the CE. The total ionisation intensity before CE is summed for all different time intervals between single and double ionisation; the results are presented in Supplementary Fig. \ref{fig7}. For the X state, a broad peak appears at approximately 2 optical cycles, whereas for the A state, an intense and well-resolved peak appears at 3 optical cycles. Thus, for the A state, the transient structure captured by TERCE mainly reflects the wave-packet distribution at approximately 8 fs after the initial tunnelling ionisation from a neutral molecule.

The evolutions of the most probable structures along the potential energy surfaces (PESs) of cation and dication can be obtained from the simulations and the results are shown as Supplementary Fig. \ref{fig8} for different time intervals between single and double ionisation. According to Supplementary Fig. \ref{fig7}, the most probable time intervals for X and A states are 2 and 3 optical cycles, respectively, thus the most probable value of $\rm R_{OD}$ and $\rm \theta_{DOD}$ at these time instants can be acquired, and their differences between neutral to cation and cation to dication are marked by the arrows in the Supplementary Fig. \ref{fig8}.  For the X state, the $\rm R_{OD}$ stretches 7 pm relative to the neutral equilibrium structure, and it further extends 5 pm during the electron recollision process along the PESs of dication, these numbers are inserted in Supplementary Fig. \ref{fig8}(a). The $\rm \theta_{DOD}$ stays similar as the neutral state. In contrast, for the A state, and $\rm \theta_{DOD}$ increases more than 40$^\circ$ and the $\rm R_{OD}$ stretches 10 pm relative to the neutral equilibrium structure, and the $\rm R_{OD}$ further extends 7 pm during the electron recollision process along the PESs of dication, as shown in Supplementary Fig. \ref{fig8}(c) and (d).





\newpage
\begin{figure}[h]
	\centering
	\includegraphics[width=0.9\linewidth]{figS4.ai}
	\caption{(Color online) Simulated time-dependent vibrational wave-packets of the X state of cation (Ca.) and the projected wave-packets distribution for dication (Dica.) and trication (Trica.) during TERCE.}
	\label{fig4}
\end{figure}

\newpage
\begin{figure}[!h]
	\centering
	\includegraphics[width=0.9\linewidth]{figS5.ai}
	\caption{(Color online) Simulated time-dependent vibrational wave-packets of the A state of cation (Ca.) and the projected wave-packets distribution for dication (Dica.) and trication (Trica.) during TERCE. }
	\label{fig5}
\end{figure}
%


\newpage
\begin{figure}[h]
	\centering
	\includegraphics[width=0.9\linewidth]{figS6.ai}
	\caption{(Color online) Calculated average $\rm R_{OD}$ and $\rm \theta_{DOD}$ based on the time-dependent vibrational wave-packet for the X and A states of D$\rm {_2}O^{+}$, respectively.}
	\label{fig6}
\end{figure}



\newpage

\begin{figure}[h]
	\centering
	\includegraphics[width=0.9\linewidth]{figS7.ai}
	\caption{(Color online)  Integrated total ionisation probability before CE, presented for different time intervals between single and double ionisation. A broad peak centred at 2 optical cycles is clearly visible for the X state((a) and (c)) ; in contrast, the A state is detected as an intense and well-resolved peak located at 3 optical cycles ((b) and (d)).   }
	\label{fig7}
\end{figure}



\newpage
\begin{figure}[h]
	\centering
	\includegraphics[width=0.9\linewidth]{figS8.ai}
	\caption{(Color online) Most probable structural evolutions for cation and dication. For different time intervals between single ionisation and double ionisation, the distributions of $\rm R_{OD}$ and $\rm \theta_{DOD}$ for X and A  states are presented, respectively. The changes of $\rm R_{OD}$ and $\rm \theta_{DOD}$ relative to the neutral equilibrium structure are marked as dashed arrows, and the extensions of $\rm R_{OD}$ during the electron recollision are marked as solid arrows in (a) and (c). 		
		}
	\label{fig8}
\end{figure}



\newpage
\section{Experimental setup }
\begin{figure}[h]
	\centering
	\includegraphics[width=\linewidth]{figS9.ai}
	\caption{(Color online)  Schematic diagram of the experimental setup that was used for TERCE. The D$\rm _2$O molecules are introduced into the chamber via supersonic expansion and irradiated by the linearly and circularly polarized laser pulses. The three ions produced via the Coulomb explosion are detected in coincidence by the microchannel plate and delay-line detector and their three dimensional momentums can be obtained. }
	\label{ fig9}
\end{figure}



%\bibliography{supplement_RE_D2O}
\begin{thebibliography}{1}
\urlstyle{rm}
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\expandafter\ifx\csname doiprefix\endcsname\relax\def\doiprefix{DOI: }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem{guillemin_selecting_2015}
\bibinfo{author}{Guillemin, R.} \emph{et~al.}
\newblock \bibinfo{journal}{\bibinfo{title}{Selecting core-hole localization or
  delocalization in {CS}{\textsubscript{2}} by photofragmentation dynamics}}.
\newblock {\emph{\JournalTitle{Nature Communications}}}
  \textbf{\bibinfo{volume}{6}}, \bibinfo{pages}{6166} (\bibinfo{year}{2015}).

\bibitem{neumann_fragmentation_2010}
\bibinfo{author}{Neumann, N.} \emph{et~al.}
\newblock \bibinfo{journal}{\bibinfo{title}{Fragmentation {dynamics} of
  {CO}{\textsubscript{2}}{\textsuperscript{3+}} {investigated} by {multiple}
  {electron} {capture} in {collisions} with {slow} {highly} {charged} {ions}}}.
\newblock {\emph{\JournalTitle{Physical Review Letters}}}
  \textbf{\bibinfo{volume}{104}}, \bibinfo{pages}{103201}
  (\bibinfo{year}{2010}).

\end{thebibliography}


\end{document}
