Supporting Information

Quantum nuclear dynamics on a distributed set of ion-trap quantum computing systems

Anurag Dwivedi
Department of Chemistry, Indiana University, Bloomington, Indiana 47405, USA
Indiana University Quantum Science and Engineering Center, Bloomington, Indiana 47405, USA
A. J. Rasmusson
Department of Physics, Indiana University, Bloomington, Indiana 47405, USA
Indiana University Quantum Science and Engineering Center, Bloomington, Indiana 47405, USA
Philip Richerme
Department of Physics, Indiana University, Bloomington, Indiana 47405, USA
Indiana University Quantum Science and Engineering Center, Bloomington, Indiana 47405, USA
Srinivasan S. Iyengar
Indiana University Quantum Science and Engineering Center, Bloomington, Indiana 47405, USA
Department of Chemistry, Indiana University, Bloomington, Indiana 47405, USA

This Supporting Information document is organized as follows. In <ref> we present our theory for multi-dimensional quantum nuclear wavepacket dynamics that naturally affords itself to implementation on a family of distributed quantum systems. This approach is based on tensor networks, and a specialized version of this general theory is implemented on the IonQ's 11-qubit trapped-ion quantum computer, Harmony, and accessed through the cloud via Amazon Braket [Wright et al., 2019, Amazon Web Services, 2020]. A key aspect of this methodology is the development of a set of entangled low-dimensional quantum systems that are implemented here through a quantum-classical interface. The potential energy surfaces used for this purpose are described in <ref>, and experimental details are provided in <ref>. To compute vibrational properties we use an approach that was described in Ref. Debadrita-Mapping-1D-3Qubits and summarized in <ref>. The range of initial conditions considered for wavepacket dynamics is discussed in <ref>. The time-evolution of classical and quantum simulations are compared in <ref>, and the spectral features derived from these are compared in <ref>.

§ MULTI-DIMENSIONAL WAVEPACKET PROPAGATION ON A DISTRIBUTED STREAM OF QUANTUM COMPUTING SYSTEMS

psi = [rectangle, very thick, minimum width=2cm, minimum height = 0.8cm, inner sep=0.1cm,draw=purple, fill=purple!20]

phi = [circle, very thick, minimum size=0.8cm, inner sep=0.1mm, draw=orange, fill=orange!20]

phinew = [circle, very thick, minimum size=0.7cm, inner sep=0.1mm, draw=bblue, fill=bblue!20]

U = [rectangle, very thick, minimum size=0.8cm, inner sep=0.1cm,draw=green, fill=green!20]

arw/.style=-Stealth, thick
[>=latex, node distance=1cm,
every edge quote/.style=font=71,scale=0.5]
(psi1) [psi, font=91]$\psi(x_1,x_2)$;

(psix01) [below=0mm of psi1, xshift = -7.5mm, minimum size=0mm, inner sep=0mm];
(psix02) [below=0mm of psi1, xshift = +7.5mm, minimum size=0mm, inner sep=0mm];

(psix00) [below=2mm of psi1, minimum size=0mm, inner sep=0mm];
(ppx01) [below=5mm of psix00, minimum size=0mm, inner sep=0mm];

(psix1) [below=2mm of psi1, xshift = -7.5mm, minimum size=0mm, inner sep=0mm];
(psix2) [below=2mm of psi1, xshift = +7.5mm, minimum size=0mm, inner sep=0mm];

(pphi1) [phi,below=8mm of psi1,xshift = -9.5mm, font=71]$\phi^{[1]}_\alpha(x_1)$;
(pphi2) [phi, below=8mm of psi1,xshift = +9.5mm, font=71]$\phi^{[2]}_\alpha(x_2)$;

(ppx1) [below=2mm of pphi1, minimum size=0mm, inner sep=0mm];
(ppx2) [below=2mm of pphi2, minimum size=0mm, inner sep=0mm];

(u1) [U, below=2mm of ppx1, font=81]$U^{[1]}_\beta$;
(u2) [U, below=2mm of ppx2, font=81]$U^{[2]}_\beta$;

(ux1) [below=3mm of u1, minimum size=0mm, inner sep=0mm];
(ux2) [below=3mm of u2, minimum size=0mm, inner sep=0mm];

(eql0sgn) [below=1mm of ux1, xshift=9mm];
(eqlsgn) [below=5mm of eql0sgn];

(phit1) [phinew, below=12mm of u1, font=71]$\varphi^{[1]}_{\alpha,\beta}(x_1)$;
(phit2) [phinew, below=12mm of u2, font=71]$\varphi^{[2]}_{\alpha,\beta}(x_2)$;

(phitx1) [below=3mm of phit1, minimum size=0mm, inner sep=0mm];
(phitx2) [below=3mm of phit2, minimum size=0mm, inner sep=0mm];

[-To, very thick, black, to path=-| ()] (psix00) – (ppx01);
[-To, very thick, black, to path=-| ()] (eql0sgn) edge (eqlsgn);

(psix01) edge[very thick, purple] node[left,black, yshift=-1mm, font=71]  (psix1)
(psix02) edge[very thick, purple]
node[left,black, yshift=-1mm,font=101]  (psix2)

(pphi1) edge[very thick, orange] node[above,black,font=71] $\alpha$ (pphi2)

(pphi1) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx1)
(pphi2) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx2)

(u1) edge[very thick, green] node[above,black,font=71] $\beta$  node[below,black,font=71]  (u2)

(u1) edge[very thick, green] (ppx1)
(u1) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (ux1)
(u2) edge[very thick, green] (ppx2)
(u2) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (ux2)

(phit1) edge[very thick, bblue] node[above,black,font=71] ${\alpha,\beta}$ node[below,black]

(phit1) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx1)
(phit2) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx2)

psi = [rectangle, very thick, minimum width=2cm, minimum height = 0.8cm, inner sep=0.1cm,draw=purple, fill=purple!20]

phi = [circle, very thick, minimum size=0.8cm, inner sep=0.2mm, draw=orange, fill=orange!20]

phinew = [circle, very thick, minimum size=0.8cm, inner sep=0.2mm, draw=bblue, fill=bblue!20]

U = [rectangle, very thick, minimum size=0.8cm, inner sep=0.1cm,draw=green, fill=green!20]

arw/.style=-Stealth, thick

[>=latex,node distance=3cm,
every edge quote/.style=font=71,

(pphi1) [phi, xshift = -8mm, font=71]$\phi^{[1]}_\alpha(x_1)$;
(pphi2) [phi, right=5mm of pphi1, font=71]$\phi^{[2]}_\alpha(x_2)$;

(ppx01) [below=5mm of psix00, minimum size=0mm, inner sep=0mm];

(ppx1) [below=2mm of pphi1, minimum size=0mm, inner sep=0mm];
(ppx2) [below=2mm of pphi2, minimum size=0mm, inner sep=0mm];

(V1) [U, below=2mm of ppx1, font=71]$\mathcal{V}^{[1]}_\beta$;
(V2) [U, below=2mm of ppx2, font=71]$\mathcal{V}^{[2]}_\beta$;

(K1) [U, below=2mm of V1, font=71]$\mathcal{K}^{[1]}$;
(K2) [U, below=2mm of V2, font=71]$\mathcal{K}^{[2]}$;

(u1) [U, below=2mm of K1, font=71]$\mathcal{V}^{[1]}_\gamma$;
(u2) [U, below=2mm of K2, font=71]$\mathcal{V}^{[2]}_\gamma$;

(ux1) [below=3mm of u1, minimum size=0mm, inner sep=0mm];
(ux2) [below=3mm of u2, minimum size=0mm, inner sep=0mm];

(eql0sgn) [below=10mm of ux1, right=7.5mm of ux1];
(eqlsgn) [below=5mm of eql0sgn];

(phit1) [phinew, below=12mm of u1, font=71]$\varphi^{[1]}_{\mu}(x_1)$;
(phit2) [phinew, below=12mm of u2, font=71]$\varphi^{[2]}_{\mu}(x_2)$;

(phitx1) [below=3mm of phit1, minimum size=0mm, inner sep=0mm];
(phitx2) [below=3mm of phit2, minimum size=0mm, inner sep=0mm];

[-To, very thick, black, to path=-| ()] (eql0sgn) edge (eqlsgn);

(pphi1) edge[very thick, orange] node[above,black,font=81] $\alpha$ (pphi2)

(pphi1) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx1)
(pphi2) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx2)

(V1) edge[very thick, green] node[above,black,font=81] $\beta$ (V2)

(u1) edge[very thick, green] node[below,black,font=81] $\gamma$ (u2)

(V1) edge[very thick, green] (ppx1)
(K1) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (u1)
(V2) edge[very thick, green] (ppx2)
(K2) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (u2)

(V1) edge[very thick, green] (K1)
(V2) edge[very thick, green] (K2)

(u1) edge[very thick, green] (ux1)
(u2) edge[very thick, green] (ux2)

(phit1) edge[very thick , bblue] node[above,black,font=81] ${\mu}$ node[below,black]

(phit1) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx1)
(phit2) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx2)

Figure <ref> depicts  the time evolution of the tensor network representation of a wavepacket. Figure <ref> enumerates this action when the time-evolution is written in Trotterized form. Here, the initial vector is given by entanglement variables $\alpha$ while the Trotterized propagator has two entangled variables, $\beta$ and $\gamma$, which combine to determine the index $\mu \equiv \left( \alpha, \beta, \gamma\right)$ for the propagated system.

psi = [rectangle, very thick, minimum width=2cm, minimum height = 0.8cm, inner sep=0.1cm,draw=purple, fill=purple!20]

phi = [circle, very thick, minimum size=0.8cm, inner sep=0.1mm, draw=orange, fill=orange!20]

phinew = [circle, very thick, minimum size=0.7cm, inner sep=0.1mm, draw=bblue, fill=bblue!20]

U = [rectangle, very thick, minimum size=0.8cm, inner sep=0.1cm,draw=green, fill=green!20]

arw/.style=-Stealth, thick
[>=latex, node distance=1cm,
every edge quote/.style=font=71,scale=0.5]
(psi1) [psi, font=91]$\psi(x_1,x_2)$;

(psix01) [below=0mm of psi1, xshift = -7.5mm, minimum size=0mm, inner sep=0mm];
(psix02) [below=0mm of psi1, xshift = +7.5mm, minimum size=0mm, inner sep=0mm];

(psix00) [below=2mm of psi1, minimum size=0mm, inner sep=0mm];
(ppx01) [below=5mm of psix00, minimum size=0mm, inner sep=0mm];

(psix1) [below=2mm of psi1, xshift = -7.5mm, minimum size=0mm, inner sep=0mm];
(psix2) [below=2mm of psi1, xshift = +7.5mm, minimum size=0mm, inner sep=0mm];

(pphi1) [phi,below=8mm of psi1,xshift = -9.5mm, font=71]$\phi^{[1]}_\alpha(x_1)$;
(pphi2) [phi, below=8mm of psi1,xshift = +9.5mm, font=71]$\phi^{[2]}_\alpha(x_2)$;

(ppx1) [below=2mm of pphi1, minimum size=0mm, inner sep=0mm];
(ppx2) [below=2mm of pphi2, minimum size=0mm, inner sep=0mm];

(u1) [U, below=2mm of ppx1, font=81]$U^{[1]}_\beta$;
(u2) [U, below=2mm of ppx2, font=81]$U^{[2]}_\beta$;

(ux1) [below=3mm of u1, minimum size=0mm, inner sep=0mm];
(ux2) [below=3mm of u2, minimum size=0mm, inner sep=0mm];

(eql0sgn) [below=1mm of ux1, xshift=9mm];
(eqlsgn) [below=5mm of eql0sgn];

(phit1) [phinew, below=12mm of u1, font=71]$\varphi^{[1]}_{\alpha,\beta}(x_1)$;
(phit2) [phinew, below=12mm of u2, font=71]$\varphi^{[2]}_{\alpha,\beta}(x_2)$;

(phitx1) [below=3mm of phit1, minimum size=0mm, inner sep=0mm];
(phitx2) [below=3mm of phit2, minimum size=0mm, inner sep=0mm];

[-To, very thick, black, to path=-| ()] (psix00) – (ppx01);
[-To, very thick, black, to path=-| ()] (eql0sgn) edge (eqlsgn);

(psix01) edge[very thick, purple] node[left,black, yshift=-1mm, font=71]  (psix1)
(psix02) edge[very thick, purple]
node[left,black, yshift=-1mm,font=101]  (psix2)

(pphi1) edge[very thick, orange] node[above,black,font=71] $\alpha$ (pphi2)

(pphi1) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx1)
(pphi2) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx2)

(u1) edge[very thick, green] node[above,black,font=71] $\beta$  node[below,black,font=71]  (u2)

(u1) edge[very thick, green] (ppx1)
(u1) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (ux1)
(u2) edge[very thick, green] (ppx2)
(u2) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (ux2)

(phit1) edge[very thick, bblue] node[above,black,font=71] ${\alpha,\beta}$ node[below,black]

(phit1) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx1)
(phit2) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx2)

The figure depicts  the time evolution of the tensor network representation of a wavepacket.

psi = [rectangle, very thick, minimum width=2cm, minimum height = 0.6cm, inner sep=0.0cm,draw=purple, fill=purple!20]

phi = [circle, very thick, minimum size=0.6cm, inner sep=0.0mm, draw=orange, fill=orange!20]

phinew = [circle, very thick, minimum size=0.6cm, inner sep=0.0mm, draw=bblue, fill=bblue!20]

phiempty = [circle, very thick, minimum size=0.cm, inner sep=0.0mm, draw=bblue, fill=bblue!20]

U = [rectangle, very thick, minimum size=0.6cm, inner sep=0.1cm,draw=green, fill=green!20]

arw/.style=-Stealth, thick

[>=latex, node distance=0.7cm,
every edge quote/.style=font=71]

(pphiempty) [phiempty, font=71];

(pphi1) [phi, left=5mm of pphiempty, font=71]$\phi^{[1]}_\alpha(x_1)$;
(pphi2) [phi, right=10mm of pphi1, font=71]$\phi^{[2]}_\alpha(x_2)$;

(ppx01) [below=5mm of psix00, minimum size=0mm, inner sep=0mm];

(ppx1) [below=5mm of pphi1, minimum size=0mm, inner sep=0mm];
(ppx2) [below=5mm of pphi2, minimum size=0mm, inner sep=0mm];

(V1) [U, below=5mm of ppx1, font=71]$\mathcal{U}^{[1]}_\beta$;
(V2) [U, below=5mm of ppx2, font=71]$\mathcal{U}^{[2]}_\beta$;

(ux1) [below=10mm of V1, minimum size=0mm, inner sep=0mm];
(ux2) [below=10mm of V2, minimum size=0mm, inner sep=0mm];

(pphi1) edge[very thick, orange] node[above,black,font=101] $\alpha$ (pphi2)

(pphi1) edge[very thick, orange] node[left,black,yshift=-1mm,font=71]  (ppx1)
(pphi2) edge[very thick, orange] node[left,black,yshift=-1mm,font=71]  (ppx2)

(V1) edge[very thick, green] node[above,black,font=101] $\beta$ (V2)

(V1) edge[very thick, green] (ppx1)
(V1) edge[very thick, green] node[left,black,yshift=-0.5mm,font=71] (ux1)
(V2) edge[very thick, green] (ppx2)
(V2) edge[very thick, green] node[left,black,yshift=-0.5mm,font=71] (ux2)

[black] (1.2,-2.3) circle (0pt)node[below]Ion-trap-2;
[black] (-1.2,-2.3) circle (0pt)node[below]Ion-trap-1;
[black] (0.0,-3) circle (0pt)node[below]In general, $\alpha*\beta$ times;
[black] (0.0,-3.5) circle (0pt)node[below]number of dimensions;
[black] (0.0,-4) circle (0pt)node[below]distributed ion traps,;
[black] (0.0,-4.5) circle (0pt)node[below]working in parallel.;

[black] (0,0) circle (3pt)node[left]$\ket{\phi^\alpha_2}$;
[black] (0,2) circle (3pt)node[left]$\ket{\phi^\alpha_1}$;
[black] (2.5,0) circle (3pt)node[below]Ion-trap-2;
[black] (2.5,2) circle (3pt)node[above]Ion-trap-1;

Derived from Figure <ref> and depicts the simultaneous study of quantum dynamics for individual one-dimensional wavefunctions on distinct ion-trap quantum computers.

We present the theory here to construct quantum nuclear wavepacket dynamics calculations on a stream of quantum processors. At the outset, the initial quantum nuclear state is expressed as the matrix product state (MPS)[Orús, 2014], where the $N$-dimensional quantum nuclear state is written as a family of one-dimensional functions,

\begin{align}
\psi(\vb{\bar{x}}) = \sum_{\bar{\boldsymbol\alpha}}^{\bar{\boldsymbol\eta}} \tensor*{\phi}{^{[1]}_{}^{}_{\alpha_1}}(x_1) \qty[\prod_{j=2}^{N\text{-} 1} \tensor*{\phi}{^{[j]}_{}^{}_{\alpha_{j\text{-} 1}}_{\alpha_j}}(x_j)] \tensor*{\phi}{^{[N]}_{}_{\alpha_{N\text{-} 1}}}(x_N),
\label{Eq:MPS_WF}
\end{align}

Here, each term $\phi^{[j]}$ denotes a one-dimensional function, with the summation in <ref> extending over the so-called entanglement variables, bond dimensions, or Schmidt rank[Orús, 2014] $\bar{\boldsymbol\alpha} \equiv \qty{ \alpha_1, \alpha_2, \ldots, \alpha_{N-1} }$. The respective limits of this summation are denoted by $\bar{\boldsymbol\eta} \equiv \qty{\eta_1,\eta_2,\ldots,\eta_{N-1}}$.

This next paragraph in its current form argues that there is no need for quantum computing and hence I am removing it.
For an $N$-dimensional system, with each dimension represented within a $D$-dimensional linear vector space, the required size of the wavepacket to accurately define its quantum state scales as $O(D^N)$. However, in the MPS form of the TN representation described in <ref>, the complexity of storing the quantum state scales as $O(ND\eta^2_{max})$, where $\eta_{max} = \max\qty{\eta_1,\eta_2,\ldots,\eta_{N-1}}$.
Therefore, for systems where entanglement scales as area of the Hilbert space as opposed to volume[Orús, 2014], that is when $\eta_{max}$ not scale exponentially with $N$, the TN formalism may mitigate the computational complexity, scaling linearly with dimension $N$, in contrast to the exponential scaling of $O(D^N)$.

In the current publication, we have two quantum nuclear degrees of freedom, that is,  $N=2$. Thus, the above expression is a Schmidt decomposition[Pathak, 2013]:

\begin{align}
\psi(x_1,x_2) = \sum_{{\alpha}_1}^{{\eta}_1} \tensor*{\phi}{^{[1]}_{}^{}_{\alpha_1}}(x_1)  \tensor*{\phi}{^{[2]}_{}_{\alpha_{1}}}(x_2),
\label{Eq:Schmidt}
\end{align}

For the two dimensional case treated here, the bond dimension is $\alpha_1$, and in the diagramatic represenatation of Eq. (<ref>) in  <ref>(a),
$\alpha_1$ is replaced by the symbol $\alpha$.
Thus, in the two dimensional case, the system is defined by a set of one-dimensional functions: $\left\{ \tensor*{\phi}{^{[1]}_{}^{}_{\alpha}}(x_1) \right\};\left\{ \tensor*{\phi}{^{[2]}_{}_{\alpha}}(x_2) \right\}$ as illustrated using orange circles in <ref>(a).

psi = [rectangle, very thick, minimum width=2cm, minimum height = 0.8cm, inner sep=0.1cm,draw=purple, fill=purple!20]

phi = [circle, very thick, minimum size=0.8cm, inner sep=0.1mm, draw=orange, fill=orange!20]

phinew = [circle, very thick, minimum size=0.7cm, inner sep=0.1mm, draw=bblue, fill=bblue!20]

U = [rectangle, very thick, minimum size=0.8cm, inner sep=0.1cm,draw=green, fill=green!20]

arw/.style=-Stealth, thick
[>=latex, node distance=1cm,
every edge quote/.style=font=71,scale=0.5]
(psi1) [psi, font=91]$\psi(x_1,x_2)$;

(psix01) [below=0mm of psi1, xshift = -7.5mm, minimum size=0mm, inner sep=0mm];
(psix02) [below=0mm of psi1, xshift = +7.5mm, minimum size=0mm, inner sep=0mm];

(psix00) [below=2mm of psi1, minimum size=0mm, inner sep=0mm];
(ppx01) [below=5mm of psix00, minimum size=0mm, inner sep=0mm];

(psix1) [below=2mm of psi1, xshift = -7.5mm, minimum size=0mm, inner sep=0mm];
(psix2) [below=2mm of psi1, xshift = +7.5mm, minimum size=0mm, inner sep=0mm];

(pphi1) [phi,below=8mm of psi1,xshift = -9.5mm, font=71]$\phi^{[1]}_\alpha(x_1)$;
(pphi2) [phi, below=8mm of psi1,xshift = +9.5mm, font=71]$\phi^{[2]}_\alpha(x_2)$;

(ppx1) [below=2mm of pphi1, minimum size=0mm, inner sep=0mm];
(ppx2) [below=2mm of pphi2, minimum size=0mm, inner sep=0mm];

(u1) [U, below=2mm of ppx1, font=81]$U^{[1]}_\beta$;
(u2) [U, below=2mm of ppx2, font=81]$U^{[2]}_\beta$;

(ux1) [below=3mm of u1, minimum size=0mm, inner sep=0mm];
(ux2) [below=3mm of u2, minimum size=0mm, inner sep=0mm];

(eql0sgn) [below=1mm of ux1, xshift=9mm];
(eqlsgn) [below=5mm of eql0sgn];

(phit1) [phinew, below=12mm of u1, font=71]$\varphi^{[1]}_{\alpha,\beta}(x_1)$;
(phit2) [phinew, below=12mm of u2, font=71]$\varphi^{[2]}_{\alpha,\beta}(x_2)$;

(phitx1) [below=3mm of phit1, minimum size=0mm, inner sep=0mm];
(phitx2) [below=3mm of phit2, minimum size=0mm, inner sep=0mm];

[-To, very thick, black, to path=-| ()] (psix00) – (ppx01);
[-To, very thick, black, to path=-| ()] (eql0sgn) edge (eqlsgn);

(psix01) edge[very thick, purple] node[left,black, yshift=-1mm, font=71]  (psix1)
(psix02) edge[very thick, purple]
node[left,black, yshift=-1mm,font=101]  (psix2)

(pphi1) edge[very thick, orange] node[above,black,font=71] $\alpha$ (pphi2)

(pphi1) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx1)
(pphi2) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx2)

(u1) edge[very thick, green] node[above,black,font=71] $\beta$  node[below,black,font=71]  (u2)

(u1) edge[very thick, green] (ppx1)
(u1) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (ux1)
(u2) edge[very thick, green] (ppx2)
(u2) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (ux2)

(phit1) edge[very thick, bblue] node[above,black,font=71] ${\alpha,\beta}$ node[below,black]

(phit1) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx1)
(phit2) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx2)

psi = [rectangle, very thick, minimum width=2cm, minimum height = 0.8cm, inner sep=0.1cm,draw=purple, fill=purple!20]

phi = [circle, very thick, minimum size=0.8cm, inner sep=0.2mm, draw=orange, fill=orange!20]

phinew = [circle, very thick, minimum size=0.8cm, inner sep=0.2mm, draw=bblue, fill=bblue!20]

U = [rectangle, very thick, minimum size=0.8cm, inner sep=0.1cm,draw=green, fill=green!20]

arw/.style=-Stealth, thick

[>=latex,node distance=3cm,
every edge quote/.style=font=71,

(pphi1) [phi, xshift = -8mm, font=71]$\phi^{[1]}_\alpha(x_1)$;
(pphi2) [phi, right=5mm of pphi1, font=71]$\phi^{[2]}_\alpha(x_2)$;

(ppx01) [below=5mm of psix00, minimum size=0mm, inner sep=0mm];

(ppx1) [below=2mm of pphi1, minimum size=0mm, inner sep=0mm];
(ppx2) [below=2mm of pphi2, minimum size=0mm, inner sep=0mm];

(V1) [U, below=2mm of ppx1, font=71]$\mathcal{V}^{[1]}_\beta$;
(V2) [U, below=2mm of ppx2, font=71]$\mathcal{V}^{[2]}_\beta$;

(K1) [U, below=2mm of V1, font=71]$\mathcal{K}^{[1]}$;
(K2) [U, below=2mm of V2, font=71]$\mathcal{K}^{[2]}$;

(u1) [U, below=2mm of K1, font=71]$\mathcal{V}^{[1]}_\gamma$;
(u2) [U, below=2mm of K2, font=71]$\mathcal{V}^{[2]}_\gamma$;

(ux1) [below=3mm of u1, minimum size=0mm, inner sep=0mm];
(ux2) [below=3mm of u2, minimum size=0mm, inner sep=0mm];

(eql0sgn) [below=10mm of ux1, right=7.5mm of ux1];
(eqlsgn) [below=5mm of eql0sgn];

(phit1) [phinew, below=12mm of u1, font=71]$\varphi^{[1]}_{\mu}(x_1)$;
(phit2) [phinew, below=12mm of u2, font=71]$\varphi^{[2]}_{\mu}(x_2)$;

(phitx1) [below=3mm of phit1, minimum size=0mm, inner sep=0mm];
(phitx2) [below=3mm of phit2, minimum size=0mm, inner sep=0mm];

[-To, very thick, black, to path=-| ()] (eql0sgn) edge (eqlsgn);

(pphi1) edge[very thick, orange] node[above,black,font=81] $\alpha$ (pphi2)

(pphi1) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx1)
(pphi2) edge[very thick, orange] node[left,black,yshift=-1mm,font=81]  (ppx2)

(V1) edge[very thick, green] node[above,black,font=81] $\beta$ (V2)

(u1) edge[very thick, green] node[below,black,font=81] $\gamma$ (u2)

(V1) edge[very thick, green] (ppx1)
(K1) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (u1)
(V2) edge[very thick, green] (ppx2)
(K2) edge[very thick, green] node[left,black,yshift=-0.5mm,font=81] (u2)

(V1) edge[very thick, green] (K1)
(V2) edge[very thick, green] (K2)

(u1) edge[very thick, green] (ux1)
(u2) edge[very thick, green] (ux2)

(phit1) edge[very thick , bblue] node[above,black,font=81] ${\mu}$ node[below,black]

(phit1) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx1)
(phit2) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=81] (phitx2)

psi = [rectangle, very thick, minimum width=2cm, minimum height = 0.6cm, inner sep=0.0cm,draw=purple, fill=purple!20]

phi = [circle, very thick, minimum size=0.6cm, inner sep=0.0mm, draw=orange, fill=orange!20]

phinew = [circle, very thick, minimum size=0.6cm, inner sep=0.0mm, draw=bblue, fill=bblue!20]

phiempty = [circle, very thick, minimum size=0.cm, inner sep=0.0mm, draw=bblue, fill=bblue!20]

U = [rectangle, very thick, minimum size=0.6cm, inner sep=0.1cm,draw=green, fill=green!20]

arw/.style=-Stealth, thick

[>=latex, node distance=0.7cm,
every edge quote/.style=font=71]

(pphiempty) [phiempty, font=71];

(pphi1) [phi, left=5mm of pphiempty, font=71]$\phi^{[1]}_\alpha(x_1)$;
(pphi2) [phi, right=10mm of pphi1, font=71]$\phi^{[2]}_\alpha(x_2)$;

(ppx01) [below=5mm of psix00, minimum size=0mm, inner sep=0mm];

(ppx1) [below=5mm of pphi1, minimum size=0mm, inner sep=0mm];
(ppx2) [below=5mm of pphi2, minimum size=0mm, inner sep=0mm];

(V1) [U, below=5mm of ppx1, font=71]$\mathcal{U}^{[1]}_\beta$;
(V2) [U, below=5mm of ppx2, font=71]$\mathcal{U}^{[2]}_\beta$;

(ux1) [below=10mm of V1, minimum size=0mm, inner sep=0mm];
(ux2) [below=10mm of V2, minimum size=0mm, inner sep=0mm];

(pphi1) edge[very thick, orange] node[above,black,font=101] $\alpha$ (pphi2)

(pphi1) edge[very thick, orange] node[left,black,yshift=-1mm,font=71]  (ppx1)
(pphi2) edge[very thick, orange] node[left,black,yshift=-1mm,font=71]  (ppx2)

(V1) edge[very thick, green] node[above,black,font=101] $\beta$ (V2)

(V1) edge[very thick, green] (ppx1)
(V1) edge[very thick, green] node[left,black,yshift=-0.5mm,font=71] (ux1)
(V2) edge[very thick, green] (ppx2)
(V2) edge[very thick, green] node[left,black,yshift=-0.5mm,font=71] (ux2)

[black] (1.2,-2.3) circle (0pt)node[below]Ion-trap-2;
[black] (-1.2,-2.3) circle (0pt)node[below]Ion-trap-1;
[black] (0.0,-3) circle (0pt)node[below]In general, $\alpha*\beta$ times;
[black] (0.0,-3.5) circle (0pt)node[below]number of dimensions;
[black] (0.0,-4) circle (0pt)node[below]distributed ion traps,;
[black] (0.0,-4.5) circle (0pt)node[below]working in parallel.;

[black] (0,0) circle (3pt)node[left]$\ket{\phi^\alpha_2}$;
[black] (0,2) circle (3pt)node[left]$\ket{\phi^\alpha_1}$;
[black] (2.5,0) circle (3pt)node[below]Ion-trap-2;
[black] (2.5,2) circle (3pt)node[above]Ion-trap-1;

Figure <ref> depicts  the time evolution of the tensor network representation of a wavepacket. Figure <ref> enumerates the action of <ref> when the time-evolution is written in Trotterized form. Here, the initial vector is given by entanglement variables $\alpha$ while the Trotterized propagator has two entangled variables, $\beta$ and $\gamma$, which combine to determine the index $\mu \equiv \left( \alpha, \beta, \gamma\right)$ for the propagated system. Figure <ref> is derived from Figures <ref> and <ref> and concisely depicts the simultaneous study of quantum dynamics for individual one-dimensional wavefunctions on distinct ion-trap quantum computers.

It is noteworthy that writing the initial wavefunction into MPS stands as a pivotal step in this approach. Not only does it reduce the computational complexity in representing a quantum state, but it also furnishes us with a distributed quantum computing protocol. This paper essentially develops a parallel stream of one-dimensional propagators mapped onto independent ion-trap systems as depicted in <ref>(c).

§.§ Action of time-evolution operator on initial MPS wavefunction

We represent the Hamiltonian for the wavepacket dynamics in the coordinate representation, where the kinetic energy operator for each nuclear dimension is separately represented using “Distributed Approximating Functionals” (DAF)[Huang et al., 1994, Hoffman et al., 1991, Iyengar, 2006],

\begin{align}
K(x_i,x_i^{\prime}) =& K(\left\vert x_i-x_i^{\prime}\right\vert) \nonumber \\ =&
\frac{-\hbar^2}{4m\sigma^3\sqrt{2\pi}}
\exp \left\{ -\frac{ {\left( x_i - x_i^\prime \right)}^2}
{2 {\sigma}^2} \right\} %\nonumber \\ &
\sum_{n=0}^{M_{DAF}/2} {\left( \frac{-1}{4} \right)}^n \frac{1}{n!}
H_{2n+2} \left( \frac{ x_i - x_i^\prime }{ \sqrt{2} \sigma} \right).
\label{DAFderivative}
\end{align}

The DAF provides an analytic banded Toeplitz representation for the
the kinetic energy operator. The banded-Toeplitz representation of the DAF approximation for the kinetic energy operator,
where the property of its matrix elements, $K_{ij} \equiv K{(\vert} i-j {\vert )}$,
has a critical role in reducing the nuclear Hamiltonian to the form of the ion-trap Hamiltonian in Ref. Debadrita-Mapping-1D-3Qubits, and also here in obtaining circuit representations in <ref>. The quantities $\left\{ H_{2n+2}\left( \frac{ x_i - x_i^\prime }{ \sqrt{2} \sigma} \right) \right\}$ are the even order Hermite polynomials that only depend on the spread separating the grid basis vectors, $\ket{x_i}$ and $\ket{x_i^{\prime}}$, and $M_{DAF}$ and $\sigma$ are parameters that together determine the accuracy and efficiency of the resultant approximate kinetic energy operator. Appendix D in Ref. Debadrita-Mapping-1D-3Qubits provides a brief summary of the DAF approach for approximating functions, their derivatives and propagated forms. Thus,  the Hamiltonian describing the system under investigation is

\begin{align}
\hat{H}\left(\vb{\bar{x}},\vb{\bar{x}}^\prime\right) = \sum_i {K}\left({{x}_i},{x}_i^\prime\right) + \delta \left( \vb{\bar{x}}-\vb{\bar{x}}^\prime \right) \hat{V}\left(\vb{\bar{x}},\vb{\bar{x}}^\prime\right),
\end{align}

where the potential energy operator, $\hat{V}\left(\vb{\bar{x}},\vb{\bar{x}}^\prime\right)$ is diagonal in the coordinate representation and is obtained from electronic structure calculations. Furthermore, the potential energy operator is not naturally separable into dimensions like in the case of the kinetic energy operator. While this does present a challenge, we present a solution in <ref>. An additional challenge in quantum dynamics arises from the need to approximate the action of the time-evolution operator on a given state vector $\Psi$, that is, $e^{-i\hat{H}\Delta t/\hbar} \ket{\Psi}$. In this study, we employ the widely used Trotter-Suzuki factorization [Campbell, 1897, Baker, 1901, Trotter, 1959, Suzuki, 1976, Nelson, 1964, Iyengar and Jakowski, 2005]. This factorization method offers a straightforward approximation of the operator, where the accuracy depends on the size of the time-step $\Delta t$. Specifically, the Trotter symmetric split operator [Trotter, 1959, Nelson, 1964] expansion of $e^{-i\hat{H}\Delta t/\hbar}$, at second order in $\Delta t$, is

\begin{align}
e^{-i\hat{H}\Delta t/\hbar} &=
e^{-i\hat{V}\left(\vb{\bar{x}}\right)\Delta t/2\hbar} \left\{ \prod_i e^{-i\hat{K}\left({{x}_i},{x}_i^\prime\right)\Delta t/\hbar}\right\} e^{-i\hat{V}\left(\vb{\bar{x}^\prime}\right)\Delta t/2\hbar} + O(\Delta t^3) \nonumber \\ &=
e^{-i\hat{V}\left(\vb{\bar{x}}\right)\Delta t/2\hbar} \left\{ \prod_i\tensor*{\mathcal{K}}{^{[i]}}(x_i,x_i^\prime) \right\} e^{-i\hat{V}\left(\vb{\bar{x}^\prime}\right)\Delta t/2\hbar} + O(\Delta t^3)
%e^{-i(\hat{K} + \hat{V})\Delta t/\hbar}\nonumber\\&= e^{-i\hat{V}\Delta t/2\hbar}e^{-i\hat{K}\Delta t/\hbar}e^{-i\hat{V}\Delta t/2\hbar} + O(\Delta t^3),
\label{Eq:Trotter}
\end{align}

that is, the time-evolution operator is now expressed in terms of the potential propagator $e^{-i\hat{V}\Delta t/2\hbar}$ and the individual kinetic propagators $e^{-i\hat{K}\left({{x}_i},{x}_i^\prime\right)\Delta t/\hbar} \equiv \tensor*{\mathcal{K}}{^{[i]}}(x_i,x_i^\prime)$ using expression given in Eq. (<ref>).

Additionally, as noted above, the non-separable nature of the potential arising from electronic structure does present a problem and here, as in Eq. (<ref>), we represent $e^{-i\hat{V}\left(\vb{\bar{x}}\right)\Delta t/2\hbar}$ as a tensor network, which for two-dimensions is given by

\begin{align}
e^{-i\hat{V}\left(\vb{\bar{x}}\right)\Delta t/2\hbar} = \sum_{\beta} \tensor*{\mathcal{V}}{^{[1]}_{\beta}^{}}(x_1) \tensor*{\mathcal{V}}{^{[2]}_{\beta}}(x_2).
\label{Eq:UV-TN-1}
\end{align}

Here, $\left\{ \tensor*{\mathcal{V}}{^{[1]}_{}^{}_{\beta}} (x_1)\right\};\left\{ \tensor*{\mathcal{V}}{^{[2]}_{}^{}_{\beta}}(x_2) \right\}$ are the one-dimensional functions obtained from the Schmidt decomposition[Pathak, 2013] of the potential propagator $e^{-i\hat{V}\left(\vb{\bar{x}}\right)\Delta t/2\hbar}$ and the sum runs over the entanglement variable $\beta$.

Finally, the time-evolution of the MPS nuclear wavefunction described in Equation <ref> is achieved through Equations <ref> and <ref> as

\begin{align}
\nonumber
\psi(\bar{\vb{x}};t+\Delta t) &=  \sum_{\alpha,\beta,\gamma}\int \dd{x'_1}\dd{x'_2}\\
&\quad\times \qty[\tensor*{\mathcal{V}}{^{[1]}_{\gamma}}(x_1) \tensor*{\mathcal{K}}{^{[1]}}(x'_1,x_1) \tensor*{\mathcal{V}}{^{[1]}_{\beta}}(x'_1)]\phi^{[1]}_{\alpha}(x'_1)\nonumber\\
&\quad\times\qty[\tensor*{\mathcal{V}}{^{[2]}_{\gamma}}(x_2) \tensor*{\mathcal{K}}{^{[2]}}(x'_2,x_2) \tensor*{\mathcal{V}}{^{[2]}_{\beta}}(x'_2)]\phi^{[2]}_{\alpha}(x'_{2})
\label{Eq:MPS_tEvo}
\end{align}

and this process is described by <ref>(b).

As a result, the action of the time-evolution operator on the initial MPS, is transformed into a set of parallel streams of one-dimensional effective quantum propagations and this is concisely represented in <ref>(c).

§.§ Computing effective 1D Hamiltonians from 1D potential propagators

psi = [rectangle, very thick, minimum width=2cm, minimum height = 0.8cm, inner sep=0.1cm,draw=purple, fill=purple!20]

phi = [circle, very thick, minimum size=0.8cm, inner sep=0.2mm, draw=orange, fill=orange!20]

phinew = [circle, very thick, minimum size=0.8cm, inner sep=0.2mm, draw=bblue, fill=bblue!20]

U = [rectangle, very thick, minimum size=0.8cm, inner sep=0.1cm,draw=green, fill=green!20]

gate = [rectangle, very thick, minimum width=1cm, minimum height = 1.8cm, inner sep=0.1cm,draw=black, fill=white!20]

arw/.style=-Stealth, thick

[>=latex,node distance=3cm,
every edge quote/.style=font=71,

[black] (-6.77,-0.5) circle (0pt)node[above]Ion-trap;
[black] (-6.77,-0.9) circle (0pt)node[above]Hamiltonian;
[black] (-6.77,-1.3) circle (0pt)node[above]Simulator-1;

[black] (7.09,-0.5) circle (0pt)node[above]Ion-trap;
[black] (7.09,-0.9) circle (0pt)node[above]Hamiltonian;
[black] (7.09,-1.3) circle (0pt)node[above]Simulator-2;

(pphi1) [phi, xshift = -8mm, font=81]$\phi^{[1]}_\alpha(x_1)$;
(pphi2) [phi, right=8mm of pphi1, font=81]$\phi^{[1]}_\alpha(x_2)$;

(ppx01) [below=5mm of psix00, minimum size=0mm, inner sep=0mm];

(ppx1) [below=2mm of pphi1, minimum size=0mm, inner sep=0mm];
(ppx2) [below=2mm of pphi2, minimum size=0mm, inner sep=0mm];

(V1) [U, below=2mm of ppx1, font=81]$\mathcal{V}^{[1]}_\beta$;
(V2) [U, below=2mm of ppx2, font=81]$\mathcal{V}^{[2]}_\beta$;

(K1) [U, below=2mm of V1, font=81]$\mathcal{K}^{[1]}$;
(K2) [U, below=2mm of V2, font=81]$\mathcal{K}^{[2]}$;

(u1) [U, below=2mm of K1, font=81]$\mathcal{V}^{[1]}_\gamma$;
(u2) [U, below=2mm of K2, font=81]$\mathcal{V}^{[2]}_\gamma$;

(pptp01) [left=2mm of V1, minimum size=0mm, inner sep=0mm];
(pptp02) [right=2mm of V2, minimum size=0mm, inner sep=0mm];

(pptp1) [left=5mm of pptp01, minimum size=0mm, inner sep=0mm];
(pptp2) [right=5mm of pptp02, minimum size=0mm, inner sep=0mm];

(kptk01) [left=2mm of K1, minimum size=0mm, inner sep=0mm];
(kptk02) [right=2mm of K2, minimum size=0mm, inner sep=0mm];

(kptk1) [left=5mm of kptk01, minimum size=0mm, inner sep=0mm];
(kptk2) [right=5mm of kptk02, minimum size=0mm, inner sep=0mm];

(uptp01) [left=2mm of u1, minimum size=0mm, inner sep=0mm];
(uptp02) [right=2mm of u2, minimum size=0mm, inner sep=0mm];

(uptp1) [left=5mm of uptp01, minimum size=0mm, inner sep=0mm];
(uptp2) [right=5mm of uptp02, minimum size=0mm, inner sep=0mm];

(LV1) [U, left=10mm of V1, font=81]$V^{[1]}_\beta=\log \mathcal{V}^{[1]}_\beta$;
(LV2) [U, right=10mm of V2, font=81]$V^{[2]}_\beta=\log \mathcal{V}^{[2]}_\beta$;

(LK1) [U, left=10mm of K1, font=81]$K^{[1]}=\log \mathcal{K}^{[1]}$;
(LK2) [U, right=10mm of K2, font=81]$K^{[2]}=\log \mathcal{K}^{[2]}$;

(Lu1) [U, left=10mm of u1, font=81]$V^{[1]}_\gamma=\log \mathcal{V}^{[1]}_\gamma$;
(Lu2) [U, right=10mm of u2, font=81]$V^{[2]}_\gamma=\log \mathcal{V}^{[2]}_\gamma$;

(ptH01) [left=2mm of LK1, minimum size=0mm, inner sep=0mm];
(ptH02) [right=2mm of LK2, minimum size=0mm, inner sep=0mm];

(ptH1) [left=5mm of ptH01, minimum size=0mm, inner sep=0mm];
(ptH2) [right=5mm of ptH02, minimum size=0mm, inner sep=0mm];

(H1) [U, left=8.2mm of LK1, font=81, align=center]$\begin{aligned}H^{[1]}_{\gamma;\beta} = &\; K^{[1]} + \\ &\frac{1}{2}({V}^{[1]}_{\gamma}+{V}^{[1]}_{\beta})\end{aligned}$;
(H2) [U, right=8.2mm of LK2, font=81, align=center]$\begin{aligned}H^{[2]}_{\gamma;\beta} = &\; K^{[2]} + \\ &\frac{1}{2}({V}^{[2]}_{\gamma}+{V}^{[2]}_{\beta})\end{aligned}$;

(H1x) [below=20mm of H1, minimum size=0mm, inner sep=0mm];
(H2x) [below=20mm of H2, minimum size=0mm, inner sep=0mm];

(phit1) [phinew, below=10mm of u1, font=81]$\varphi^{[1]}_{\mu}(x_1)$;
(phit2) [phinew, below=10mm of u2, font=81]$\varphi^{[2]}_{\mu}(x_2)$;

(phitx1) [below=3mm of phit1, minimum size=0mm, inner sep=0mm];
(phitx2) [below=3mm of phit2, minimum size=0mm, inner sep=0mm];

(phitl) [left=5mm of phit1, minimum size=0mm, inner sep=0mm];
(phitr) [right=5mm of phit2, minimum size=0mm, inner sep=0mm];

[-To, very thick, black, to path=-> ()] (pptp01) edge (pptp1);
[-To, very thick, black, to path=-> ()] (pptp02) edge (pptp2);
[-To, very thick, black, to path=-> ()] (kptk01) edge (kptk1);
[-To, very thick, black, to path=-> ()] (kptk02) edge (kptk2);
[-To, very thick, black, to path=-> ()] (uptp01) edge (uptp1);
[-To, very thick, black, to path=-> ()] (uptp02) edge (uptp2);
[-To, very thick, black, to path=-> ()] (ptH01) edge (ptH1);
[-To, very thick, black, to path=-> ()] (ptH02) edge (ptH2);

(ux1) [below=3mm of u1, minimum size=0mm, inner sep=0mm];
(ux2) [below=3mm of u2, minimum size=0mm, inner sep=0mm];

(pphi1) edge[very thick, orange] node[above,black,font=101] $\alpha$ (pphi2)

(pphi1) edge[very thick, orange] node[left,black,yshift=-1mm,font=101]  (ppx1)
(pphi2) edge[very thick, orange] node[left,black,yshift=-1mm,font=101]  (ppx2)

(V1) edge[very thick, green] node[above,black,font=101] $\beta$ (V2)

(u1) edge[very thick, green] node[below,black,font=101] $\gamma$ (u2)

(V1) edge[very thick, green] (ppx1)
(u1) edge[very thick, green] node[left,black,yshift=-0.5mm,font=101] (ux1)
(V2) edge[very thick, green] (ppx2)
(u2) edge[very thick, green] node[left,black,yshift=-0.5mm,font=101] (ux2)

(V1) edge[very thick, green] (K1)
(V2) edge[very thick, green] (K2)

(phit1) edge[very thick, bblue] node[above,black,font=101] ${\mu}$ node[below,black]

(phit1) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=101] (phitx1)
(phit2) edge[very thick, bblue] node[left,black,yshift=-0.5mm,font=101] (phitx2)

(u1) edge[very thick, green] (K1)
(u2) edge[very thick, green] (K2)

[gate, left=38mm of phitl, yshift=4mm, font=71] (gateU) $U=e^{\frac{-\iota H^{[1]}_{\gamma;\beta}t}{\hbar}}$;

(q1) [above=5mm of gateU, minimum size=0mm, inner sep=1mm,font=71]$\ket{q_1}$;
(q0) [left=0.1mm of q1, minimum size=0mm, inner sep=1mm,font=71]$\ket{q_0}$;
(q2) [right=0.1mm of q1, minimum size=0mm, inner sep=1mm,font=71]$\ket{q_2}$;

(q1u) [below=5mm of q1, minimum size=0mm, inner sep=0mm];
(q0u) [left=6.5mm of q1u, minimum size=0mm, inner sep=0mm];
(q2u) [right=6.5mm of q1u, minimum size=0mm, inner sep=0mm];

(p1) [below=5mm of gateU, minimum size=0mm, inner sep=0mm];
(p0) [left=6.5mm of p1, minimum size=0mm, inner sep=0mm];
(p2) [right=6.5mm of p1, minimum size=0mm, inner sep=0mm];

(p1u) [above=5mm of p1, minimum size=0mm, inner sep=0mm];
(p0u) [left=6.5mm of p1u, minimum size=0mm, inner sep=0mm];
(p2u) [right=6.5mm of p1u, minimum size=0mm, inner sep=0mm];

(q0u) edge[very thick, black](q0)
(q1u) edge[very thick, black](q1)
(q2u) edge[very thick, black](q2)

(p0u) edge[very thick, black](p0)
(p1u) edge[very thick, black](p1)
(p2u) edge[very thick, black](p2)

[gate, right=38mm of phitr, yshift=4mm,font=71] (gateUr) $U=e^{\frac{-\iota H^{[2]}_{\gamma;\beta}t}{\hbar}}$;

(q1r) [above=5mm of gateUr, minimum size=0mm, inner sep=1mm,font=71]$\ket{q_1}$;
(q0r) [left=0.1mm of q1r, minimum size=0mm, inner sep=1mm,font=71]$\ket{q_0}$;
(q2r) [right=0.1mm of q1r, minimum size=0mm, inner sep=1mm,font=71]$\ket{q_2}$;

(q1ur) [below=5mm of q1r, minimum size=0mm, inner sep=0mm];
(q0ur) [left=6.5mm of q1ur, minimum size=0mm, inner sep=0mm];
(q2ur) [right=6.5mm of q1ur, minimum size=0mm, inner sep=0mm];

(p1r) [below=5mm of gateUr, minimum size=0mm, inner sep=0mm];
(p0r) [left=6.5mm of p1r, minimum size=0mm, inner sep=0mm];
(p2r) [right=6.5mm of p1r, minimum size=0mm, inner sep=0mm];

(p1ur) [above=5mm of p1r, minimum size=0mm, inner sep=0mm];
(p0ur) [left=6.5mm of p1ur, minimum size=0mm, inner sep=0mm];
(p2ur) [right=6.5mm of p1ur, minimum size=0mm, inner sep=0mm];

(q0ur) edge[very thick, black](q0r)
(q1ur) edge[very thick, black](q1r)
(q2ur) edge[very thick, black](q2r)

(p0ur) edge[very thick, black](p0r)
(p1ur) edge[very thick, black](p1r)
(p2ur) edge[very thick, black](p2r)

(fsl) [left=20mm of phitl, minimum size=0mm, inner sep=0mm, font=71];
(fsr) [right=20mm of phitr, minimum size=0mm, inner sep=0mm, font=71];

[-To, very thick, black, to path=-> ()] (fsl) edge (phitl);
[-To, very thick, black, to path=-> ()] (fsr) edge (phitr);

A detailed exposition of <ref> that precisely outlines how the distributed simulations are carried out. The individual one-dimensional potentials are derived from each one-dimensional propagator as shown, yielding a family of effective one-dimensional Hamiltonians for the system. These individual one-dimensional Hamiltonians are then be independently and concurrently simulated on separate ion-trap quantum computers.

We first note that the direct action of each Trotterized propagator block in <ref> and in <ref>(b) depicts the action of one-dimensional propagators, obtained from tensor-network decomposition of the potential-propagator,

on the MPS components of the initial wavepacket.
Towards simplifying this action and making direct analogues to ion-trap Ising Hamiltonains, as discussed in Ref. Debadrita-Mapping-1D-3Qubits, we write each of the potential-propagator terms, $\left\{\mathcal{V}^{[j]}_{\beta_{1}}(x_{j})\right\};\left\{ \mathcal{V}^{[j]}_{\gamma_{1}}(x_{j}) \right\}$, encoding information about the one dimensional potentials as,

\begin{align}
\nonumber\mathcal{V}^{[j]}_{\beta}(x_{j})
={\cal A}^{[j]}_{\beta}(x_{j}) \times \exp{-\imath {V}^{[j]}_{\beta}(x_{j}) \Delta t/2\hbar}\\
\mathcal{V}^{[j]}_{\gamma}(x_{\textit{\textit{\textit{\textit{}}}}j})
={\cal A}^{[j]}_{\gamma}(x_{j}) \times \exp{-\imath {V}^{[j]}_{\gamma}(x_{j}) \Delta t/2\hbar}
\label{Eq:V-A-S}
\end{align}

where, $\left\{{V}^{[j]}_{\beta}(x_{j})\right\};\left\{ {V}^{[j]}_{\gamma}(x_{j}) \right\}$
are effective one dimensional potentials obtained by taking the logarithmic of the potential propagators, that is,

\begin{align}
{\log \cal A}^{[j]}_{\beta}(x_{j}) &= Re \left[ \log \mathcal{V}^{[j]}_{\beta}(x_{j}) \right] \label{Eq:log_A}\\
{V}^{[j]}_{\beta}(x_{j}) &= -\frac{2\hbar}{\Delta t} Im \left[ \log \mathcal{V}^{[j]}_{\beta}(x_{j}) \right]
\label{Eq:log_V}
\end{align}

to obtain a family of reduced dimensional Hamiltonians

\begin{align}
H^{[j]}_{\gamma;\beta}(x_j, x'_j) =K^{[j]}\left({{x}_j},{x}_j^\prime\right)+\frac{1}{2}({V}^{[j]}_{\gamma}(x_{i})+{V}^{[j]}_{\beta}(x'_{j}))
\label{Eq:MPS_Ham}
\end{align}

as depicted within the pink and blue blocks in <ref>.

Using these in Eq. (<ref>), we obtain,

\begin{align}
\nonumber
\psi(\bar{\vb{x}};t+\Delta t) &=  \sum_{\alpha,\beta,\gamma}\int \dd{x'_1}\dd{x'_2}\\
&\quad\times \qty[\splitfrac{{\cal A}^{[1]}_{\gamma}(x_{1})\exp{-\imath {V}^{[1]}_{\gamma}(x_{1}) \Delta t/2\hbar}\exp{-iK^{[1]}(x_1, x'_1) \Delta t/\hbar}}%
{{\cal A}^{[1]}_{\beta}(x'_{1})\exp{-\imath {V}^{[1]}_{\beta}(x'_{1})\Delta t/2\hbar}}]\phi^{[1]}_{\alpha}(x'_1)\nonumber\\
&\quad\times\qty[\splitfrac{{\cal A}^{[2]}_{\gamma}(x_{2})\exp{-\imath {V}^{[2]}_{\gamma}(x_{2}) \Delta t/2\hbar}\exp{-iK^{[2]}(x_2, x'_2) \Delta t/\hbar}}%
{{\cal A}^{[2]}_{\beta}(x'_{2})\exp{-\imath {V}^{[2]}_{\beta}(x'_{2}) \Delta t/2\hbar}}]\phi^{[2]}_{\alpha}(x'_{2})
\label{Eq:MPS_tEvo-LogHam1}
\end{align}

Additionally, the quantities $\left\{{\cal A}^{[j]}{\beta}(x_j); {\cal A}^{[j]}{\gamma}(x_{j})\right\}$ in <ref> are responsible for redistributing probabilities among the various physical dimensions while ensuring the overall process remains unitary. These tensors, originating from the potential-based evolution operator (see <ref>), facilitate a crucial step of dissipating and regenerating probabilities among the propagated Schmidt vectors for each dimension. Notably, the potential energy surface (PES) on a $8\times8$ grid depicted in <ref>, satisfy the conditions ${\cal A}^{[j]}{\beta}(x_j) =1$ and ${\cal A}^{[j]}{\gamma}(x_{j}) =1$. Substituting these conditions in <ref> enable the unitary evolution of each effective one-dimensional system onto separate ion-traps quantum hardware, as illustrated in <ref>.

ψ(x̅;t+Δt)    =  ∑_α,β,γ∫x'_1x'_2

×[exp-V^[1]_γ(x_1) Δt/2ħexp-iK^[1](x_1, x'_1) Δt/ħ

×[exp-V^[2]_γ(x_2) Δt/2ħexp-iK^[2](x_2, x'_2) Δt/ħ
exp-V^[2]_β(x'_2) Δt/2ħ]ϕ^[2]_α(x'_2)

Hence, the aforementioned equations illustrate the Trotterized action (see <ref>) of the time-evolution operator acting independently on each dimension.

In this publication, ${\cal A}^{[j]}_{\beta}(x_{j}) = 1$, and hence <ref> does not contribute to the trajectories.
This feature enables the parallel computation of nuclear dynamics across distinct ion-trap quantum hardware. Consequently, the TN representation establishes effective one-dimensional subsystems, thereby defining corresponding `one-dimensional' Hamiltonians for each of these subsystems,

ψ(x̅;t+Δt)    = ∑_α,β,γ∫x'_1x'_2

×[exp-H_γ;β^[1](x_1, x'_1)Δt/ħ]ϕ^[1]_α(x'_1)

×[exp-H^[2]_γ;β(x_2, x'_2)Δt/ħ]ϕ^[2]_α(x'_2).

where, the effective `one-dimensional' Hamiltonians obey,

\begin{align}
\exp{-\imath H^{[j]}_{\gamma;\beta}(x_j, x'_j)\Delta t/\hbar} \equiv
\exp{{-i\left[K^{[j]}+
\frac{
\left({V}^{[j]}_{\gamma}(x_{j})+{V}^{[j]}_{\beta}(x'_{j})\right)
\Delta t/\hbar}}
\label{Eq:Ham1}
\end{align}

\begin{align}
\exp{-\imath H^{[j]}_{\gamma;\beta}(x_j, x'_j)\Delta t/\hbar} \equiv
\exp{{-i\left[K^{[j]}(x_j, x'_j)+
\frac{1}{2}\left({V}^{[j]}_{\gamma}(x_{j})+{V}^{[j]}_{\beta}(x'_{j})\right)\right]
\Delta t/\hbar}}
\label{Eq:Ham1}
\end{align}

The equations above produce a stream of one-dimensional unitary propagations that may be directly mapped to quantum simulators as depicted in <ref>.

The one dimensional Hamiltonians, $\left\{ H_{\gamma;\beta}^{[1]}\right\};\left\{H_{\gamma;\beta}^{[2]} \right\}$ represent independent and parallel, quantum computing steps, to be performed on independent quantum hardware systems (see <ref>) and each of these propagate the representative one-dimensional systems.

§.§ Block diagonalization of one-dimensional nuclear Hamiltonians, $\left\{ H_{\gamma;\beta}^{[1]}; H_{\gamma;\beta}^{[2]} \right\}$

As discussed in Ref. Debadrita-Mapping-1D-3Qubits, the choice of kinetic energy operator in Eq. (<ref>) and symmetric nature of potential energy operator, allows for the one-dimensional nuclear Hamiltonians, $\left\{ H_{\gamma;\beta}^{[1]}; H_{\gamma;\beta}^{[2]} \right\}$, to be represented in block-diagonal form.

This process aims to facilitate the independent propagation of each block across an additional parallel set of ion-trap systems. To confirm, while the algorithm presented above (see Figure <ref>) already yields as many multiple channels as entanglement dimensions times the number of physical dimensions, the block-diagonalization discussed here yields an additional factor of two in terms of number of distributed set of quantum processes.

The kinetic energy operator, when expressed in the Distributed Approximating Functional (DAF) representation [Kouri et al., 1995, Hoffman et al., 1991], takes the form of a banded Toeplitz structure where each sub-diagonal element has a constant value, as mentioned in  <ref>. Additionally, the Schmidt decomposed potentials, $\left\{ {V}^{[j]}_{\gamma}(x_{i})\right\}$

derived from <ref>

are symmetric as we note in Section <ref>.
Matrices with these two key properties (the Toeplitz form of the kinetic energy operator and the symmetric nature of the effective one-dimensional potentials) can be transformed into a block-diagonal form by applying a Givens unitary transformation[Saha et al., 2021]. Specifically, the $(i,l)$-the element of $H_{\gamma;\beta}^{[1]}$ and $H_{\gamma;\beta}^{[2]}$ may be transformed as per

\begin{align}
\tilde{H}{^{[j]}_{{\gamma;\beta},i,l}}= \frac{1}{2}(H^{[j]}_{i,l} + \alpha_{l}H^{[j]}_{i,n-l} + \alpha_{i}H^{[j]}_{n-i,l} + \alpha_{i}\alpha_{l}H^{[j]}_{n-i,n-l})
\end{align}

Here, $\alpha_i = \text{sgn}\left[i-\frac{n}{2}\right]$, and the grid point indices range from 0 to $n=2^{D}-1$, where $D$ represents the number of qubits used to simulate the respective one-dimensional Hamiltonian. The diagonal blocks of the unitary transformed Hamiltonian $\tilde{H}^{[j]}$ can then be expressed as:

\begin{align}
\tilde{H}^{[j]}_{i,l} = \left[K^{[j]}(x_i,x_l) + \alpha_{i}K(x_{i},x_{n-l})\right] + \frac{1}{2}\left[V(x_i)+V(x_{n-l})\right]\delta_{i,l}
\label{Eq:on-diagonal}
\end{align}

Meanwhile, the elements of the off-diagonal blocks are given by:

\begin{align}
\tilde{H}^{[j]}_{i,l} = \frac{1}{2}\left[V(x_i)+V(x_{n-l})\right]\delta_{i,l}
\label{Eq:off-diagonal}
\end{align}

In systems with symmetric potentials, as found in many hydrogen-bonded systems <cit.>, the off-diagonal blocks of the transformed Hamiltonian (<ref>) are uniformly zero. This formalism allows to simulate a $2^D\times2^D$ nuclear Hamiltonian on two-different ion-traps by recasting it into a Hamiltonian with two sub-blocks of size $2^{D-1}\times2^{D-1}$, with simulating each sub-block independently onto two ion-traps. Each sub-block of the transformed Hamiltonian operate on the transformed grid basis given by:
This formalism enables the simulation of a $2^D\times2^D$
nuclear Hamiltonian across two distinct ion-traps by restructuring it into a Hamiltonian featuring two sub-blocks, each of size $2^{D-1}\times2^{D-1}$. This restructuring allows for the independent simulation of each sub-block on separate ion-traps. The grid basis for each sub-block of the transformed Hamiltonian is defined as follows:

\begin{align}
\ket{\tilde{x}_{i}} &\equiv
\begin{cases}
\frac{1}{\sqrt{2}}\left(\ket{x_{i}} + \ket{x_{n-i}}\right) & , \quad \text{for } 0 \leq i < \frac{n+1}{2}\\
\frac{1}{\sqrt{2}}\left(\ket{x_{i}} - \ket{x_{n-i}}\right) & , \quad \text{for } \frac{n+1}{2} \leq i \leq n
\end{cases}
\label{Givens-transformed-basis}
\end{align}

The pictorial representation of the above transformed grid equation is shown in Figure 5 of Ref. Debadrita-Mapping-1D-3Qubits. This formulation enables a systematic transformation of the Hamiltonian, preserving its essential properties while facilitating parallel simulation of each diagonal blocks on quantum hardware.

Furthermore, the quantities: ${\cal A}^{[j]}_{\beta}(x_j); {\cal A}^{[j]}_{\gamma}(x_{j})$
in <ref> have a very critical role in the entire process. Equations. (<ref>)perform unitary evolution on the individual MPS Schmidt vectors for each dimension. Following this, the tensors in Eq. (<ref>)  redistribute probabilities among the physical various dimensions such that the overall process is still unitary. Since the the tensors in Eq. (<ref>) arise from the potential-based evolution operator, (see Eq. (<ref>)) this redistribution of probabilities amounts to a dissipation/regeneration step between the various propagated Schmidt vectors for each dimension. The effect of such dissipation/regeneration can be easily understood if we begin with a completely uncorrelated, initial product state.
Equation (<ref>) then reduces to

\begin{align}
&\psi(\bar{\vb{x}};t+\Delta t) =
\nonumber \\ & \sum_{\bar{\vb{\beta}},\bar{\vb{\gamma}}}\qty{{\cal A}^{[1]x_{1}}_{\gamma_{1}}{\cal A}^{[1]x_{1}}_{\beta_{1}} \exp{-\imath \tilde{H}_{\gamma_1;\beta_1}^{[1]x_{1}}\Delta t/\hbar}}\phi^{[1]x_1} \nonumber \\
& \times \prod_{j=2}^{N-1} \left\{{\cal A}^{[j]x_{j}}_{\gamma_{j-1};\gamma_{j}}{\cal A}^{[j]x_{j}}_{\beta_{j-1};\beta_{j}} \exp{-\imath \tilde{H}_{\gamma_{j-1};\gamma_{j};\beta_{j-1};\beta_{j}}^{[j]x_{j}}\Delta t/\hbar} \right\} \phi^{[j]x_j} \nonumber \\
& \times {\cal A}^{[N]x_{N}}_{\gamma_{N}}{\cal A}^{[N]x_{N}}_{\beta_{N}}\qty{\exp{-\imath \tilde{H}^{[N]x_{N}}_{\gamma_N;\beta_N}\Delta t/\hbar}}\phi^{[N]x_N}
\end{align}

where the $\bar{\vb{\alpha}}$ index has now been omitted. If we now insist that the tensors in Eq. (<ref>) are all identity, then the propagation reduces to

\begin{align}
&\psi(\bar{\vb{x}};t+\Delta t) =
\nonumber \\ & \sum_{\bar{\vb{\beta}},\bar{\vb{\gamma}}}\qty{\exp{-\imath \tilde{H}_{\gamma_1;\beta_1}^{[1]x_{1}}\Delta t/\hbar}}\phi^{[1]x_1} \nonumber \\
& \times \prod_{j=2}^{N-1} \left\{ \exp{-\imath \tilde{H}_{\gamma_{j-1};\gamma_{j};\beta_{j-1};\beta_{j}}^{[j]x_{j}}\Delta t/\hbar} \right\} \phi^{[j]x_j} \nonumber \\
& \times \qty{\exp{-\imath \tilde{H}^{[N]x_{N}}_{\gamma_N;\beta_N}\Delta t/\hbar}}\phi^{[N]x_N}
\label{prop-noA}
\end{align}

where each dimension is independently norm-conserving and hence there is no population redistribution as a function of time. Furthermore, Eq. (<ref>) then reduces to

\begin{align}
&\tensor*{\mathcal{V}}{*^{x'_1}_{x_1}^{\ldots}_{\ldots}^{x'_N}_{x_N}} \nonumber \\ & = \sum_{\vb{\bar{\boldsymbol\beta}}}
\exp{-\imath {\cal W}^{[N]x_{N}}_{\beta_{N-1}} \Delta t/2\hbar} \prod_{j=2}^{N-1} \exp{-\imath {\cal W}^{[j]x_{j}}_{\beta_{j\text{-} 1}\beta_{j}} \Delta t/2\hbar} \nonumber \\ & \times  \exp{-\imath {\cal W}^{[N]x_{N}}_{\beta_{N-1}} \Delta t/2\hbar} = \exp{-\imath {\cal W}^{\bar{x}}_{\bar{\beta}} \Delta t/2\hbar}
\label{pot-prop-MPS-2}
\end{align}

\begin{align}
{\cal W}^{\bar{x}}_{\bar{\beta}} = {\cal W}^{[N]x_{N}}_{\beta_{N-1}} + \sum_{j=2}^{N-1} {\cal W}^{[j]x_{j}}_{\beta_{j\text{-} 1}\beta_{j}} + {\cal W}^{[N]x_{N}}_{\beta_{N-1}}
\label{pot-prop-MPS-S}
\end{align}

which is an additive (and thus uncorrelated) potential. The resultant state is still entangled in the sense that a measurement at variable $x_1$, probabilistically constrains the sum over $\beta_1$ which affects the remaining states as is clear from

\begin{align}
&\bra{\chi^{[1]x_1}}\ket{\psi(\bar{\vb{x}};t+\Delta t)} =
\nonumber \\ & \sum_{\bar{\vb{\beta}},\bar{\vb{\gamma}}}\bra{\chi^{[1]x_1}}\ket{\qty{\exp{-\imath \tilde{H}_{\gamma_1;\beta_1}^{[1]x_{1}}\Delta t/\hbar}}\phi^{[1]x_1}} \nonumber \\
& \times \prod_{j=2}^{N-1} \left\{ \exp{-\imath \tilde{H}_{\gamma_{j-1};\gamma_{j};\beta_{j-1};\beta_{j}}^{[j]x_{j}}\Delta t/\hbar} \right\} \phi^{[j]x_j} \nonumber \\
& \times \qty{\exp{-\imath \tilde{H}^{[N]x_{N}}_{\gamma_N;\beta_N}\Delta t/\hbar}}\phi^{[N]x_N}
\label{prop-noA-proj}
\end{align}

Here, consistent with the definition of entanglement, we have projected the final state of $x_1$ onto some signature state, $\ket{\chi^{[1]x_1}}$, and then we ask of the final result, function of $\left[ x_2,\cdots,x_N \right]$, bears a signature of $\ket{\chi^{[1]x_1}}$. The distinction between Eq. (<ref>) and the appropriate $\ket{\chi^{[1]x_1}}$-projected form of Eq. (<ref>) is that in Eq. (<ref>) depending on the distribution of the projecting state in the $\phi^{[1]x_1}_{\alpha_1}$ state space, the result of such a projection is different. There exists no such flexibility in Eq. (<ref>). This essentially emphasizes the entanglement provided by the dissipation/regeneration terms, across $\left[ x_1,\cdots,x_N \right]$ in Eq. (<ref>), and that provided in its absence in Eq. (<ref>).

The presence of the dissipation/regeneration terms in Eq. (<ref>) alters the description, where the individual dimensional propagators in Eq. (<ref>) are not unitary anymore. This provides an entanglement that is fundamentally different as compared to Eq. (<ref>) through the ability to allow population flux between the dimensions.

Figure (a) shows the protonated 2,2'-bipyridine system with the two dimensions, $x_1$ and $x_2$, correspond to the vibrational degree of freedom of transferring proton along the internuclear axis joining the two nitrogen atoms, and the torsional degree of freedom due to planar rotations of the pyridyl rings about the C-C bond. Figure (b) shows the colormap plot of the potential energy surface for the molecule.

§ POTENTIAL ENERGY SURFACES (PES) FOR QUANTUM DYNAMICS

The protonated 2,2'-bipyridine system with the hydrogen bond dimension, $x_1$ and gating dimension, $x_2$, noted.

In this section, we provide a concise overview of the construction process for the Potential Energy Surface (PES) (see <ref>(b)) tailored to the protonated 2,2’-bipyridine system (shown in <ref>(a)), crucial for elucidating the TN propagation scheme. To capture the system's relevant degrees of freedom accurately, we employ cylindrical coordinates $(x_1,x_2)$, where $x_1$ signifies the relative position of the shared hydrogen nucleus along the NN axis, and $x_2$ denotes the dihedral angle between the two pyridyl rings surrounding the C-C bond.

The PES is constructed by identifying stationary points for each fixed dihedral angle $x_2\in[\ang{-35},\ang{35}]$, positioning the shared proton at the midpoint of the internuclear axis between the nitrogen nuclei. Notably, each of these torsional geometries yields a transition state with one imaginary frequency upon diagonalizing the Hessian matrix. This imaginary frequency corresponds to the vibrational mode associated with proton transfer along the internuclear axis connecting the nitrogen nuclei.

No. of grid points along the $x_1$-dimension      8

No. of grid points along the $x_2$-dimension     8

Grid size along the $x_1$-dimension     1.1Å

Range of angles along the $x_2$-dimension     $[-35^{\circ}$,$35^{\circ}]$

Level of theory     B3LYP/6-311++G(d,p)//B3LYP/6-311++G(d,p)

Characteristics of the grid over which 2D potential surface is created.

Subsequently, maintaining a fixed torsional geometry (that is, a constant value of $x_2$), we compute a series of one-dimensional effective potentials. These potentials, contingent upon the proton's position along the internuclear axis, are computed via a density functional theory (DFT) methodology. Specifically, we employ a hybrid generalized gradient approximation (GGA) functional B3LYP in conjunction with an atom-centered Gaussian basis set 6-311++G(d,p). To compute the Potential Energy Surface (PES), we employ a grid consisting of $2^3=8$ points along the internuclear axis connecting the two nitrogen atoms for each of the 8 specific dihedral angles within the range of [-35°, 35°]. We provide a summary of the primary features utilized in the calculation of the PES in <ref>.

The resultant rectangular $8\times 8$ grid along $(x_1,x_2)$ has potential values as shown in <ref>(b). This potential has inversion symmetry, $V(x_1,x_2) = V(-x_1,-x_2)$, and is a result of the symmetry in the molecule shown in <ref>(a). This symmetry of the potential is exploited here, using the methods discussed in Ref. Debadrita-Mapping-1D-3Qubits, to block-diagonalize the molecular Hamiltonian an obtain an efficient quantum computing algorithm. The process is discussed below.

§.§ Block diagonalization of one-dimensional nuclear Hamiltonians, $\left\{ H_{\gamma;\beta}^{[1]}\right\};\left\{ H_{\gamma;\beta}^{[2]} \right\}$ owing to symmetric nature of PES in <ref>(b)

As discussed in Ref. Debadrita-Mapping-1D-3Qubits, the choice of kinetic energy operator in Eq. (<ref>) and symmetric nature of potential energy operator highlighted above, allows for the one-dimensional nuclear Hamiltonians, $\left\{ H_{\gamma;\beta}^{[1]}\right\};\left\{ H_{\gamma;\beta}^{[2]} \right\}$, to be represented in block-diagonal form.
This process aims to facilitate the independent propagation of each block across an additional parallel set of ion-trap systems.
While the algorithm presented above (see Figure <ref>(c)) already yields as many multiple channels as entanglement dimensions times the number of physical dimensions, the block-diagonalization discussed here yields an additional factor of two in terms of number of distributed set of quantum processes.
The kinetic energy operator, when expressed in the Distributed Approximating Functional (DAF) representation [Kouri et al., 1995, Hoffman et al., 1991], takes the form of a banded Toeplitz structure where each sub-diagonal element has a constant value, as mentioned in  <ref>. Additionally, the Schmidt decomposed potentials, $\left\{ {V}^{[j]}_{\gamma}(x_{i})\right\}$ derived from <ref> are symmetric as we note in <ref>.
Matrices with these two key properties (the Toeplitz form of the kinetic energy operator and the symmetric nature of the one-dimensional effective potentials) can be transformed into a block-diagonal form by applying a Givens unitary transformation[Saha et al., 2021]. Specifically, the $(i,l)$-the element of $H_{\gamma;\beta}^{[1]}$ and $H_{\gamma;\beta}^{[2]}$ may be transformed as per

\begin{align}
\tilde{H}{^{[j]}_{{\gamma;\beta},i,l}}= \frac{1}{2}(H^{[j]}_{{\gamma;\beta},i,l} + \alpha_{l}H^{[j]}_{{\gamma;\beta},i,n_j-l} + \alpha_{i}H^{[j]}_{{\gamma;\beta},n_j-i,l} + \alpha_{i}\alpha_{l}H^{[j]}_{{\gamma;\beta},n_j-i,n_j-l})
\end{align}

Due to the symmetric nature of the potential energy surface, as discussed in <ref>, and consequently, the symmetric nature of the effective one-dimensional potential noted in <ref>, combined with the choice of kinetic energy operator mentioned in <ref>, the one-dimensional nuclear Hamiltonians $\left\{ H_{\gamma;\beta}^{[1]}; H_{\gamma;\beta}^{[2]} \right\}$ can be conveniently represented in a block-diagonal form[Saha et al., 2021]. This procedure is designed to enable the independent propagation of each block across an additional parallel set of ion-trap systems. While the algorithm presented above (see <ref>) already yields as many multiple channels as entanglement dimensions times the number of physical dimensions, the block-diagonalization discussed here yields an additional factor of two in terms of number of distributed set of quantum processes.

The kinetic energy operator, when formulated in the Distributed Approximating Functional (DAF) representation [Kouri et al., 1995, Hoffman et al., 1991], takes a banded Toeplitz structure with constant values along each sub-diagonal element, as mentioned in <ref>. Furthermore, the Schmidt decomposed potentials, $\left\{{V}^{[j]}_{\gamma}(x_{i})\right\}$ derived from <ref>, are symmetric as we note in <ref>. Matrices featuring these two crucial properties—the Toeplitz form of the kinetic energy operator and the symmetric nature of the effective one-dimensional potentials—can be transformed into a block-diagonal form through the application of a Givens unitary transformation [Saha et al., 2021]. Specifically, the $(i,l)$-the element of $H_{\gamma;\beta}^{[1]}$ and $H_{\gamma;\beta}^{[2]}$ may be transformed as follows:
Because of the symmetric nature of potential energy surface as noted in <ref>, and hence symmetric nature of effective one-dimensional potential as derived in <ref>, and with the choice of kinetic energy operator in <ref>, allows for the one-dimensional nuclear Hamiltonians, $\left\{ H_{\gamma;\beta}^{[1]}; H_{\gamma;\beta}^{[2]} \right\}$, to be represented in block-diagonal form. This process aims to facilitate the independent propagation of each block across an additional parallel set of ion-trap systems. To confirm, while the algorithm presented above (see Figure <ref>) already yields as many multiple channels as entanglement dimensions times the number of physical dimensions, the block-diagonalization discussed here yields an additional factor of two in terms of number of distributed set of quantum processes.
The kinetic energy operator, when expressed in the Distributed Approximating Functional (DAF) representation [Kouri et al., 1995, Hoffman et al., 1991], takes the form of a banded Toeplitz structure where each sub-diagonal element has a constant value, as mentioned in  <ref>. Additionally, the Schmidt decomposed potentials, $\left\{ {V}^{[j]}_{\gamma}(x_{i})\right\}$ derived from <ref> are symmetric as we note in <ref>. Matrices with these two key properties (the Toeplitz form of the kinetic energy operator and the symmetric nature of the effective one-dimensional potentials) can be transformed into a block-diagonal form by applying a Givens unitary transformation[Saha et al., 2021]. Specifically, the $(i,l)$-the element of $H_{\gamma;\beta}^{[1]}$ and $H_{\gamma;\beta}^{[2]}$ may be transformed as per

\begin{align}
\tilde{H}{^{[j]}_{{\gamma;\beta},i,l}}= \frac{1}{2}(H^{[j]}_{i,l} + \alpha_{l}H^{[j]}_{i,n-l} + \alpha_{i}H^{[j]}_{n-i,l} + \alpha_{i}\alpha_{l}H^{[j]}_{n-i,n-l})
\end{align}

Here, $\alpha_i = \text{sgn}\left[i-\frac{n_j}{2}\right]$, and the grid point indices range from 0 to $n_j=2^{D_j}-1$, where $D_j$ represents the number of qubits used to simulate the respective one-dimensional Hamiltonian. The diagonal blocks of the unitary transformed Hamiltonian $\tilde{H}_{\gamma;\beta}^{[j]}$ can then be expressed as:

\begin{align}
\tilde{H}^{[j]}_{\gamma;\beta,i,l} = \left[K^{[j]}_{i,l} + \alpha_{i}K^{[j]}_{i,n_j-l}\right] + \frac{1}{4}\left[(V^{[j]}_{\gamma,i}+V^{[j]}_{\beta,i})+(V^{[j]}_{\gamma,n_j-l}+V^{[j]}_{\beta,n_j-l})\right]\delta_{i,l}
\label{Eq:on-diagonal}
\end{align}

\begin{align}
\tilde{H}^{[j]}_{\gamma;\beta,i,l} = \left[K^{[j]}(x_i,x_l) + \alpha_{i}K^{[j]}(x_{i},x_{n-l})\right] + \frac{1}{2}\left[V(x_i)+V(x_{n-l})\right]\delta_{i,l}
\label{Eq:on-diagonal}
\end{align}

Meanwhile, the elements of the off-diagonal blocks are given by:

\begin{align}
\tilde{H}^{[j]}_{\gamma;\beta,i,l} = \frac{1}{4}\left[(V^{[j]}_{\gamma,i})-V^{[j]}_{\gamma,n_j-l})+(V^{[j]}_{\beta,i}-V^{[j]}_{\beta,n_j-l})\right]\delta_{i,l}
\label{Eq:off-diagonal}
\end{align}

Due to the inherent symmetry of the potential energy surface, as discussed in <ref>, and observed in other hydrogen-bonded systems[Sumner and Iyengar, 2007, Vendrell et al., 2007, Diken et al., 2005], the off-diagonal blocks of the transformed Hamiltonian (<ref>) are uniformly zero. In addition to the strategies outlined in <ref>, the symmetry of the potential here, enables further parallelization of simulation of these effective one-dimensional Hamiltonians.

This formalism enables the simulation of a $2^D\times2^D$
nuclear Hamiltonian across two distinct ion-traps by restructuring it into a Hamiltonian featuring two sub-blocks, each of size $2^{D-1}\times2^{D-1}$. This restructuring allows for the independent simulation of each sub-block on separate ion-traps. The grid basis for each sub-block of the transformed Hamiltonian is defined as follows:

\begin{align}
\ket{\tilde{x}_{i}} &\equiv
\begin{cases}
\frac{1}{\sqrt{2}}\left(\ket{x_{i}} + \ket{x_{n-i}}\right) & , \quad \text{for } 0 \leq i < \frac{n+1}{2}\\
\frac{1}{\sqrt{2}}\left(\ket{x_{i}} - \ket{x_{n-i}}\right) & , \quad \text{for } \frac{n+1}{2} \leq i \leq n
\end{cases}
\label{Givens-transformed-basis}
\end{align}

The pictorial representation of the above transformed grid equation is shown in Figure 5 of Ref. Debadrita-Mapping-1D-3Qubits.

Simulation type    4cInitial Wavepackets[Four different initial wavepackets were considered as shown on the four columns below. Notation: $\delta(x_1-x_1^1)$ implies the initial wavepacket has a non-zero value at grid point “1” (the superscript) for dimension “1” (the subscript)]

$x_1$-upper block[Using the block diagonalization scheme discussed in Section <ref>.]      $\delta(x_1-x_1^1)$     $\delta(x_1-x_1^2)$     $\delta(x_1-x_1^3)$     $\delta(x_1-x_1^4)$

$x_1$-lower block[2]      $\delta(x_1-x_1^5)$     $\delta(x_1-x_1^6)$     $\delta(x_1-x_1^7)$     $\delta(x_1-x_1^8)$

$x_2$-upper block[2]      $\delta(x_2-x_2^1)$     $\delta(x_2-x_2^2)$     $-$     $-$

$x_2$-lower block[2]     $\delta(x_2-x_2^6)$     $\delta(x_2-x_2^7)$     $\delta(x_2-x_2^8)$     $-$
$x_1$-full[Without the block diagonalization scheme discussed in Section <ref>. That is using the full Hamiltonian without block diagonalization.]     $\delta(x_1-x_1^1)$     $\delta(x_1-x_1^2)$     $-$     $-$

$x_2$-full[3]     $\delta(x_2-x_2^1)$     $-$     $-$     $-$

Initial wavepacket  details. The $x_1$ and $x_2$ simulations are conducted on a distributed set of two ion-traps

Hamiltonian     Size     4c|Quantum     4cClassical

(lr)3-4 (lr)5-6
$\Delta t$ (fs)     Total time (fs)     $\Delta \omega$ (THz)     Max. frequency (THz)     $\Delta t$ (fs)     Total time (fs)     $\Delta \omega$ (THz)     Max. frequency (THz)

$x_1$-upper/lower block     $4\times 4$     $2.50$     $400$     1.250     200.0     $0.25$     $400$     -     -

$x_2$-upper/lower block     $4\times 4$    $1.47$     $235$     2.125     340.4     $0.25$     $235$    -     -

$x_1$-full Hamiltonian     $8\times 8$    $100.00$     $16000$     0.031     5.0     $0.25$     $16000$    -     -

$x_2$-full Hamiltonian     $8\times 8$    $7.00$     $1120$     0.445     71.4     $0.25$     $1120$    -     -

Simulation details.

Simulation type     Size     4c|Quantum Simulation     4cClassical Simulation

(lr)3-6 (lr)7-10
$\Delta t$
(fs)       Total
time (fs)      $\Delta \omega$[Nyquist theorem]
(THz)      Max.
frequency[Max. frequency = $\frac{1}{2\Delta t}$]
(THz)      $\Delta t$
(fs)       Total
time (fs)      $\Delta \omega$[1]
(THz)      Max.

$x_1$-upper/lower block     $4\times 4$     2.50     400     1.250     200.0     0.25     400     1.250     2000

$x_2$-upper/lower block     $4\times 4$     1.47     235     2.125     340.4     0.25     235     2.125     2000

$x_1$-full Hamiltonian     $8\times 8$     100.00     16000     0.031     5.0     0.25     16000     0.031     2000

$x_2$-full Hamiltonian     $8\times 8$     7.00     1120     0.445     71.4     0.25     1120     0.445     2000

Simulation details.

Hamiltonian      $\Delta \omega$ (THz)     Max. frequency (THz)

$x_1$-upper/lower block     $1.250$     $200.0$

$x_2$-upper/lower block     $2.125$     $340.4$

$x_1$-full Hamiltonian    $0.031$     $5.0$

$x_2$-full Hamiltonian    $0.445$     $71.4$

Fourier transformation characteristics.

§ INITIAL CONDITIONS FOR WAVEPACKET DYNAMICS

Simulation type     4cInitial Wavepackets[Four different initial wavepackets were considered as shown on the four columns below. Notation: $\delta(x_1-x_1^1)$ implies the initial wavepacket has a non-zero value at grid point “1” (the superscript) for dimension “1” (the subscript)]

$x_1$-upper block[Using the block diagonalization scheme discussed in Section <ref>.]     $\delta(x_1-x_1^1)$     $\delta(x_1-x_1^2)$     $\delta(x_1-x_1^3)$     $\delta(x_1-x_1^4)$

$x_1$-lower block[2]     $\delta(x_1-x_1^5)$     $\delta(x_1-x_1^6)$     $\delta(x_1-x_1^7)$     $\delta(x_1-x_1^8)$

$x_2$-upper block[2]     $\delta(x_2-x_2^1)$     $\delta(x_2-x_2^2)$     $-$     $-$

$x_2$-lower block[2]     $\delta(x_2-x_2^6)$     $\delta(x_2-x_2^7)$     $\delta(x_2-x_2^8)$     $-$

$x_1$-full[Without the block diagonalization scheme discussed in Section <ref>. That is using the full Hamiltonian without block diagonalization.]     $\delta(x_1-x_1^1)$     $\delta(x_1-x_1^2)$     $-$     $-$

$x_2$-full[3]     $\delta(x_2-x_2^1)$     $-$     $-$     $-$

Initial wavepacket  details. The $x_1$ and $x_2$ simulations are conducted on a distributed set of two ion-traps

The dynamics of many initial wavepackets were studied and these initial wavepackets are shown in Table <ref>. In all cases, the initial wavepackets were chosen as Dirac delta functions localized at specific grid points as shown in Table <ref>. The reason for this is as follows. Since our goal was to obtain spectral properties from the dynamics, it is needed that all eigenstates be captured within a specific initial wavepacket state. Delta functions by their choice have this property. Secondly, delta functions typically have a high energy expectation value and hence in some sense these choices provide the hardest test cases for our simulation. As we see from the agreement in vibrational properties, the quantum simulations are highly robust.

Simulation type     Size     4c|Quantum Simulation     4cClassical Simulation

$\Delta t$
$(\unit{\femto\second})$       $T$
$(\unit{\femto\second})$       $\Delta \omega$[Nyquist theorem]
$(\unit{\tera\hertz})$       $\omega_{\text{max}}$[$\omega_{\text{max}} = \frac{1}{2\Delta t}$]
$(\unit{\tera\hertz})$       $\Delta t$
$(\unit{\femto\second})$       $T$
$(\unit{\femto\second})$       $\Delta \omega$[1]
$(\unit{\tera\hertz})$       $\omega_{\text{max}}$[2]

$x_1$-upper/lower block     $4\times 4$     2.50     400     1.250     200.0     0.25     400     1.250     2000

$x_2$-upper/lower block     $4\times 4$     1.47     235     2.125     340.4     0.25     235     2.125     2000

$x_1$-full Hamiltonian     $8\times 8$     100.00     16000     0.031     5.0     0.25     16000     0.031     2000

$x_2$-full Hamiltonian     $8\times 8$     7.00     1120     0.445     71.4     0.25     1120     0.445     2000

Simulation details.

Furthermore, given the block-diagonalization scheme discussed above, two separate sets of simulations were conducted as indicated by the top half and bottom half of Table <ref>. In essence, the symmetry of the potential allows us to further reduce complexity. Each of these initial wavepackets, depending on which Hamiltonian blocks they are associated with, undergo simulation on IonQ's ion-trap quantum hardware. Details regarding the total propagation time and simulation time-steps are provided in <ref>.

Additionally, the time evolution of the initial wavepacket on classical hardware is computed as follows:

\begin{align}
\nonumber
\ket{\psi_{\gamma;\beta}(x_{j};t)} &= e^{-iH^{[j]}_{\gamma;\beta}t/\hbar}\ket{\psi_{\gamma;\beta}(x_{j};0)} \\ &= \sum_i e^{-iE^{[j]}_{\gamma;\beta,i}t/\hbar}\ket{\chi^{[j]}_{\gamma;\beta,i}}\braket{\chi^{[j]}_{\gamma;\beta,i}}{\psi_{\gamma;\beta}(x_{j};0)},
\label{Eq:Classical}
\end{align}

Here, $H_{\gamma;\beta}^{[j]}$ represents the effective one-dimensional Hamiltonians (as derived in <ref>). $\ket{\chi^{[j]}_{\gamma;\beta,i}}$ represents the eigenfunctions that satisfy the eigenvalue equation $H^{[j]}_{\gamma;\beta}\ket{\chi_{\gamma;\beta,i}} = E^{[j]}_{\gamma;\beta,i}\ket{\chi^{[j]}_{\gamma;\beta,i}}$ with eigenvalue $E^{[j]}_{\gamma;\beta,i}$. The ket $\ket{\psi_{\gamma;\beta}(x_{j};0)}$ represents the initial wavepackets as detailed in <ref> for each simulation.

On the contrary for the quantum computation, multiple quantum circuits ${\cal U}(n\Delta t)$ are created and these act on the initial wavepacket to arrive at the wavepacket at a certain time step. In principle this last step can also be parallelized but is not considered in this publication.

§ VIBRATIONAL SPECTRAL BEHAVIOR FROM TIME-EVOLUTION OF ANY INITIAL STATE

The computation of the shared-proton wavepacket dynamics provides a means to accurately determine its vibrational frequencies and is given by the fourier transform of the density-density time auto-correlation function, ${\text{Tr}}[\rho(0)\rho(t)]$ as

\begin{align}
\int_{-\infty}^{+\infty} dt\, %\exp{{\imath \omega t}}
e^{{\imath \omega t}} \; {\text{Tr}}[\rho(0)\rho(t)] &= \int_{-\infty}^{+\infty} dt\, e^{{\imath \omega t}}
%\exp{{\imath \omega t}}
\; {\text{Tr}}\left[\ket{\chi(0)}\bra{\chi(0)}\sum_{i,j}c_{i}(0)c_{j}^{*}(0) %\exp(-\frac{i(E_{i}-E_{j})t}{\hbar})
e^{\imath (E_{i}-E_{j}) t/{\hbar}}
\ket{\phi_{i}}\bra{\phi_{j}}\right] \nonumber \\
&=\sum_{i,j}|c_{i}(0)|^{2}|c_{j}(0)|^{2} \; \delta \left(\omega-(E_{i}-E_{j})/\hbar \right),
\label{Eq:Density-timecorrelation-FT}
\end{align}

where time-dependent density matrix is $\rho(t) = \ket{\chi(t)}\bra{\chi(t)}$. Then, one may rewrite the right side of <ref> using convolution theorem[Press et al., 1992] as described in Ref. Sandia-1D-3Qubits, to obtain

\begin{align}
& \int dx  {
\left\{
\left\vert
\int_{-\infty}^{+\infty} dt  e^{{\imath \omega t}} \;
%\right. \right. \nonumber \\ & \left. \left.
\sum_{i,j}c_{i}(0)c_{j}^{*}(0) %\exp\left(-\frac{i(E_{i}-E_{j})t}{\hbar}\right)
e^{\imath (E_{i}-E_{j}) t/{\hbar}}
\right\vert
\right\}
}^2 +
\nonumber \\ &
\int_{x\neq x^\prime} dx \,dx^\prime \left\vert \int_{-\infty}^{+\infty} dt  e^{{\imath \omega t}}
%\right. \nonumber \\ & \phantom{\int_{x\neq x^\prime} dx dx^\prime \vert }\left.
\sum_{i,j}c_{i}(0)c_{j}^{*}(0) %\exp\left(-\frac{i(E_{i}-E_{j})t}{\hbar}\right)
e^{\imath (E_{i}-E_{j}) t/{\hbar}}
{\phi_{j}(x^\prime)} \right\vert^2.
\label{Density-timecorrelation-FT-2}
\end{align}

\begin{align}
\int_{-\infty}^{+\infty} dt\,
e^{{\imath \omega t}} \; {\text{Tr}}[\rho(0)\rho(t)]  &= \int dx  {
\left\{
\left\vert
\int_{-\infty}^{+\infty} dt  e^{{\imath \omega t}} \; \rho(x,x;t)
\right\vert
\right\}
}^2 \nonumber \\ &+
\int_{x\neq x^\prime} dx \,dx^\prime \left\vert \int_{-\infty}^{+\infty} dt  e^{{\imath \omega t}} \; \rho(x,x^\prime;t) \right\vert^2.
\label{Eq:Density-timecorrelation-FT-2}
\end{align}

While standard approaches may compute the results from Eq. (<ref>), a quantum computer, such as the ones utilized in Ref. Sandia-1D-3Qubits and in this publication, have the ability to obtain wavepacket density information (that is the diagonal elements of the wavepacket density, $\rho(x,x;t)$) and compute the individual terms inside the first term of the integral above as,

\begin{align}
{\cal I} (\omega; x) & = \int_{-\infty}^{+\infty} dt\,  e^{{\imath \omega t}} \; \rho(x,x; t) \\
& = \int_{-\infty}^{+\infty} dt\,  e^{{\imath \omega t}} \;
%\right. \right. \nonumber \\ & \left. \left.
\sum_{i,j}c_{i}(0)c_{j}^{*}(0) %\exp\left(-\frac{i(E_{i}-E_{j})t}{\hbar}\right)
e^{\imath (E_{i}-E_{j}) t/{\hbar}}
{\phi_{j}(x)} \nonumber \\
& = \sum_{i,j} \delta(\omega-(E_{i}-E_{j})) c_{i}(0)c_{j}^{*}(0) {\phi_{i}(x)}
\label{Eq:Density-timecorrelation-FT-3}
\end{align}

which makes the first term in <ref>,

\begin{align}
{\cal P} (\omega) = \int dx \; {\left\vert {\cal I} (\omega; x) \right\vert}^2
\label{Eq:integ-Density-timecorrelation-FT-3}
\end{align}

By measuring <ref> on an ion-trap quantum computer in Ref. Sandia-1D-3Qubits, we have provided a new approach to determine spectroscopic features in complex systems using quantum computing platforms.

However, as can be seen from <ref>, there are a variety of initial wavepackets used here, and hence each resultant trajectory would independently lead to the quantities in <ref>. Thus, in the results section, we have first individual computed the quantity in <ref> to obtain ${\cal I} (\omega; x) \rightarrow {\cal I}_{\chi_0} (\omega; x)$ and ${\cal P} (\omega) \rightarrow {\cal P}_{\chi_0} (\omega)$, where the respective quanities are now designated using the choice of initial wavepacket. We then compute the cumulative spectrum from all of these wavepackets and denote these as

\begin{align}
\sum_{\forall \chi_0} {{\cal P}_{\chi_0} (\omega)} \rightarrow {\cal P} (\omega)
\label{Eq:integ-Density-timecorrelation-FT-4}
\end{align}

to compute vibrational properties. Thermal weights could also be introduced at this stage but it is not done in this paper since these would only affect the peak heights and not the frequencies.

x full 000 111

x full 001 110

x upper 00

x upper 01

x upper 10

x upper 11

x lower 00

x lower 01

x lower 10

x lower 11

theta full 000 111

theta upper 00

theta upper 01

theta lower 01

theta lower 10

theta lower 11

power vs freq x upper

power vs freq x lower

power vs freq x full

power vs freq theta upper

power vs freq theta lower

power vs freq theta lower

§ EXPERIMENTAL IMPLEMENTATION ON A DISTRIBUTED SET OF ION-TRAP SYSTEMS

Experiments are performed on IonQ's 11-qubit trapped-ion quantum computer, Harmony, and accessed through the cloud via Amazon Braket [Wright et al., 2019, Amazon Web Services, 2020]. Harmony's laser beam systems can separately address each ion, or pairs of ions, enabling arbitrary single-qubit gate rotations and all-to-all two-qubit gate connectivity. The wavepacket dynamics are simulated and measured on Harmony with 1000 repetitions for each time point. Quantum circuits are executed in parallel such that the maximum number of qubits are in use. Finally, a classical Fourier transform of the data is used to reveal the frequency components of the vibrational motion using the formalism mentioned in <ref>.

Quantum circuits are compiled directly from the target unitary operators. This approach reduces error accumulation due to continuous analog simulation or due to Trotterization errors which may degrade the spectral signal in the frequency domain. Quantum circuits have the same gate layout for a given Hamiltonian, but vary the single-qubit rotation angles to simulate evolution for a given time. The largest contribution to errors comes from two-qubit gate operations (largest by an order of magnitude) across all time points. Variations in error due to single-qubit gate rotations are negligible. Since we aim to measure vibrational spectra, time-varying errors in the dynamics are more impactful than constant-amplitude errors. We observe that in our compilation strategy, two-qubit gate errors lead to a notable loss in oscillation amplitudes, but the frequency domain spectra remain robust.

\begin{quantikz}
\lstick{} & \ctrl{1} & \gate{R_Z(\pi/2) \; \sqrt{X} \; R_Z(\alpha)} & \ctrl{1} & \gate{R_Z(\pi/2) \; \sqrt{X} \; R_Z(\pi/2)} & \ctrl{1} & \gate{\sqrt{X}} & \qw \\
\lstick{} & \targ{} & \gate{R_Z(\gamma)} & \targ{} & \gate{R_Z(\beta)} & \targ{} & \gate{R_Z(\pi) \; \sqrt{X} \; R_Z(\pi)} & \qw\\
\end{quantikz}
General quantum circuit diagram for simulating the dynamics of the two-qubit block Hamiltonians ($x_1$-upper/lower and $x_2$-upper/lower). The parameters $\gamma$, $\beta$, and $\alpha$ are varied per time point.

The two-qubit unitary time evolution blocks of the $x_1$-upper/lower and $x_2$-upper/lower Hamiltonians are compiled into quantum circuits of CNOT gates and single-qubit rotations using Cartan's $KAK$ decomposition, which guarantees a maximum of three two-qubit gates given our system size [Tucci, 2005, Earnest et al., 2021]. As mentioned previously, the single-qubit rotations encode the time point values. A generalized quantum circuit diagram is shown in Fig. <ref>. Lastly, the circuit is compiled into the native gate basis supported by Harmony as described in Ref. wright2019benchmarking. A two-qubit evolution for any two-qubit block may look like the following circuit where single-qubit gate parameters $\gamma$, $\beta$, and $\alpha$ depend on the evolution time $t$[Earnest et al., 2021].

[]        $
\begin{quantikz}
\lstick{} & \gate[wires=2][1.5cm]{U(t_i)} & \gate[wires=3][1.5cm]{D(t_i)} & \gate[wires=2][1.5cm]{V(t_i)} & \qw \\
\lstick{} & & & & \qw \\
\lstick{} & \qw & & \gate{H} & \qw \\
\end{quantikz}
[]        $
\begin{quantikz}
\lstick{} & \qw & \qw & \qw & \ctrl{2} & \qw & \qw & \qw & \ctrl{2} & \qw \\
\lstick{} & \qw & \ctrl{1} & \qw & \qw & \qw & \ctrl{1} & \qw & \qw & \qw \\
\lstick{} & \gate{R_Z(\theta_0)} & \targ{} & \gate{R_Z(\theta_1)} & \targ{} & \gate{R_Z(\theta_2)} & \targ{} & \gate{R_Z(\theta_3)} & \targ{} & \qw \\
\end{quantikz}
(a) The general quantum circuit diagram for simulating the full three-qubit block Hamiltonian dynamics ($x_1$-full and $x_2$-full). State preparation and measurement are not shown. The additional Hadamard gate maps the basis back to the grid basis as described in the text. (b) A generalized quantum circuit diagram for the diagonal unitary $D(t_i)$.

The three-qubit unitary time evolution of the $x_1$-full and $x_2$-full Hamiltonians is compiled as follows, with a general quantum circuit diagram shown in <ref>. First, the Single Value Decomposition (SVD) separates the unitary operator of the given Hamiltonian at a specified time point $t_i$ into three parts $U(t_i)$, $D(t_i)$, and $V(t_i)$ as given by

\begin{equation}
\begin{pmatrix}U_0 & 0 \\ 0 & U_1\end{pmatrix} = \begin{pmatrix} U & 0 \\ 0 & U \end{pmatrix} \begin{pmatrix} D & 0 \\ 0 & D^{\dagger} \end{pmatrix} \begin{pmatrix} V & 0 \\ 0 & V \end{pmatrix}
\end{equation}

where $U_0$ and $U_1$ are the known two-qubit upper and lower blocks of the unitary evolution and $0$'s are $4 \times 4$ blocks. Since $U_0 = UDV$ and $U_1 =  UD^{\dagger}V$,

\begin{equation}
U_0 U_1^{\dagger} = U D^2 U^{\dagger} \;.
\end{equation}

From the eigenvalue decomposition of $U_0 U_1^{\dagger}$, we can solve for $U$ and $D$. We find $V$ using the relation $V = D^{-1} U^{-1} U_0$ defined earlier. Next, these matrices are complied as elements of the quantum circuit. Cartan's $KAK$ decomposition is used for the two-qubit $U(t_i)$ and $V(t_i)$ unitaries where each contributes at most three two-qubit gates and some single-qubit rotations. The remaining unitary $D(t_i)$ is diagonal which is decomposed into four two-qubit gates and four single-qubit $R_Z(\theta)$ rotations as shown in Fig. <ref>(b) [Shende et al., 2005].

Computation         Givens-transformed grid       Shuffled

$\ket{000}$     $(\vert x^0\rangle + \ket{x^7})/\sqrt{2}$       $(\vert x^0\rangle + \ket{x^7})/\sqrt{2}$

$\ket{001}$     $(\vert x^1\rangle + \ket{x^6})/\sqrt{2}$       $(\vert x^1\rangle + \ket{x^6})/\sqrt{2}$

$\ket{010}$     $(\vert x^2\rangle + \ket{x^5})/\sqrt{2}$       $(\vert x^2\rangle + \ket{x^5})/\sqrt{2}$

$\ket{011}$     $(\vert x^3\rangle + \ket{x^4})/\sqrt{2}$       $(\vert x^3\rangle + \ket{x^4})/\sqrt{2}$

$\ket{100}$     $(\vert x^4\rangle - \ket{x^3})/\sqrt{2}$       $(\vert x^7\rangle - \ket{x^0})/\sqrt{2}$

$\ket{101}$     $(\vert x^5\rangle - \ket{x^2})/\sqrt{2}$       $(\vert x^6\rangle - \ket{x^1})/\sqrt{2}$

$\ket{110}$     $(\vert x^6\rangle - \ket{x^1})/\sqrt{2}$       $(\vert x^5\rangle - \ket{x^2})/\sqrt{2}$

$\ket{111}$     $(\vert x^7\rangle - \ket{x^0})/\sqrt{2}$       $(\vert x^4\rangle - \ket{x^3})/\sqrt{2}$

Enumeration of the different bases. The computation basis are the bitstrings measured in the experiment. The Givens-transformed grid basis is in terms of the grid basis states $\ket{x^i}$. The shuffled basis is used experimentally to reduce the quantum circuit. It is similar to the Givens-transformed grid basis with the eighth and fifth states swapped and the seventh and sixth states swapped.

With this procedure of generating a quantum circuit given one of the three-qubit Hamiltonians ($x_1$-full or $x_2$-full), the dynamics can be simulated in the grid basis on the trapped-ion quantum computer. The final step is to append a series of gates which map the wavepacket onto the grid basis. We hand-optimize this circuit by adding just a single Hadamard gate by shuffling the lower block of states: the eighth and fifth states are switched, as are the seventh and sixth states, as shown in the third column of Table <ref>.

For example, a Hadamard acting on the second qubit of the $\ket{000}$ state, in the shuffled basis, gives a one-to-one mapping between the computation basis and the grid basis $\ket{x^i}$:

\begin{align}
H \ket{000} &= (\vert 000\rangle + \ket{100})/\sqrt{2} \\
& = (\vert x^0\rangle + \ket{x^7} + \ket{x^7} - \ket{x^0})) / 2 \\
& = \ket{x^7}
\end{align}

This generalizes to all states. Working in the shuffled basis has an additional simplification for state preparation. The states of interest in the transformed basis are entangled states which require two-qubit gates. However, in the shuffled basis, the initial states are no longer entangled states and single-qubit rotations can prepare them.

§ COMPARISON OF WAVEPACKET TIME-EVOLUTION FOR ALL INITIAL WAVEPACKETS FROM QUANTUM AND CLASSICAL SIMULATIONS

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^1_1)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^2_1)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^3_1)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^4_1)$]
Illustrating the coherent evolution of the wavepacket pertaining to the upper block of the Hamiltonian along the $x_1$ direction, encapsulating the vibrational dynamics along the internuclear axis connecting the two nitrogen atoms within the molecular structure depicted in <ref>. The temporal evolutions depicted in each subfigure is initialized at distinct grid points as delineated in the subfigure captions, with each time-traces projected onto specific grid points denoted in the panel titles. Solid yellow lines denote the exact numerical solution to the Schrödinger equation, calculated classically, while the blue curve represents outcomes derived from ion-trap dynamics.

This section compares the quantum dynamics trajectories obtained using the classical and quantum hardware for

the initial wavepackets provided in <ref>. As shown in Figs. <ref>(c) and <ref>, the quantum dynamics simulation for each of the effective one-dimensional Hamiltonians, $H^{[j]}_{\gamma,\beta}$, has been conducted independently. Since these one-dimensional Hamiltonians can be further expressed in a block-diagonal form, as discussed in <ref>, each of these diagonal sub-blocks is simulated independently on both quantum and classical computers. The information regarding simulation time-steps and total propagation time is delineated in <ref>.

<ref> complement figure 5 in the paper and provide a complete picture of our study. The <ref> shows the dynamics for the upper-block of the block diagonal Hamiltonian, $\tilde{H}^{[1]}_{\gamma,\beta}$, for the $x_1$ dimension, encoding the proton vibration along the donor-acceptor axis connecting the two nitrogen atoms in the structure shown in <ref>(a). Specifically, <ref>(a) presents the projection of the time-traces at distinct grid points ($x^1_1, x^2_1, x^3_1, x^4_1$) in 4 distinct panels. Note on notation: while the subscripts here represents dimensions, the superscripts refers to grid points. Hence, the four time traces in <ref> refer to the time evolution of various nuclear wavepacket on each of the four grid points. The blue time dots in <ref> are obtained from measurements on the ion-trap quantum computer and the yellow traces are from classical computation. Here, the initial wavepacket is localized at grid point `1' and has a zero value at all other grid points, that is, $\Psi_0(x_1) = \delta(x_1-x^1_1)$. Similarly, <ref>(b-d) shows the time-traces when the initial wavepackets are $\Psi_0(x_1) = \delta(x_1-x^2_1)$, $\Psi_0(x_1) = \delta(x_1-x^3_1)$ and $\Psi_0(x_1) = \delta(x_1-x^4_1)$ respectively.

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^5_1)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^6_1)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^7_1)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^8_1)$]

Showcasing the coherent evolution of the wavepacket associated with the lower block of the Hamiltonian along the $x_1$ direction, capturing the vibrational dynamics along the internuclear axis linking the two nitrogen atoms within the molecular structure portrayed in <ref>. The temporal evolution depicted in each subfigure is initiated from distinct grid points, as specified in the subfigure captions, with each time-trace projected onto specific grid points indicated in the panel titles.

The dynamics of the lower block of the block diagonal Hamiltonian $\tilde{H}^{[1]}_{\gamma,\beta}$, for the $x_1$ dimension, are illustrated in fig. <ref>. Specifically, <ref>(a) presents the projection of the time-traces at distinct grid points ($x^5_1, x^6_1, x^7_1, x^8_1$) in four distinct panels. Note: Grid points 5-8 appear on the lower block of the Hamiltonian. Here, the initial wavepacket is localized at grid point '5' and has a zero value at all other grid points, that is, $\Psi_0(x_1) = \delta(x_1-x^5_1)$. Similarly, <ref>(b-d) show the time-traces when the initial wavepackets are $\Psi_0(x_1) = \delta(x_1-x^6_1)$, $\Psi_0(x_1) = \delta(x_1-x^7_1)$ and $\Psi_0(x_1) = \delta(x_1-x^8_1)$ respectively. Once again, solid yellow lines denote the exact numerical solution to the Schrödinger equation, calculated classically, while the blue dots represents outcomes derived from ion-trap quantum computation. These $4\times4$ upper and lower block of the one-dimensional Hamiltonian $\tilde{H}^{[1]}_{\gamma,\beta}$ are simulated on a 2-qubit ion-trap systems. Detailed information regarding the time-steps and total propagation time for each of the time-traces can be found in <ref>.

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^1_1)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_1}$ for initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^2_1)$]

The coherent evolution of the wavepacket associated with the full $8\times8$ Hamiltonian along the $x_1$ direction as shown in the molecular structure in <ref>. The temporal evolution depicted in Figure (a) is initialized at the first grid point while the wavepacket in Figure (b) is initialized at the second grid point. The time-traces in each panel correspond to the time-evolution of the wavepacket projected onto specific grid points indicated in the panel titles.

As noted above, the symmetry of the potential allows the block diagonalization of the effective one-dimensional Hamiltonians. Time-evolution of the separate blocks along with the their respective Fourier transforms (discussed in the next section), yield information regarding the vibrational states for the system being studied. However, one of the critical features of Fourier analysis provided in <ref> is that only eigenvalue differences are obtained from the Fourier transform. This in general does not present a problem since most chemical spectroscopy experiments also yield eigenvalue differences. However, here, obtaining eigenvalue differences for each of the Hamiltonian blocks separately implies that the overall spectrum of the system is not accessible. Hence, to obtain the relative placements of the eigenvalues of the separate blocks, we also compute here the time-evolution of the system arising from the non-block-diagonalized Hamiltonian. The resultant quantum dynamics simulation of the $8\times8$ full-Hamiltonian, $\tilde{H}^{[1]}_{\gamma,\beta}$, has been conducted on a 3-qubit ion-trap system. As noted, this simulation aims to determine the relative energy differences between the eigenstates obtained from the Fourier transform of the dynamics simulation of the upper and lower blocks of the block-diagonal Hamiltonian. In <ref>, we illustrate the quantum dynamics of the Hamiltonian, $\tilde{H}^{[1]}_{\gamma,\beta}$, for two different initial wavepackets: $\Psi_0(x_1) = \delta(x_1-x^1_1)$ and $\Psi_0(x_1) = \delta(x_1-x^2_1)$,
as shown in <ref>(a-b) respectively. Each of these figures presents the projection of the time-traces at distinct grid points ($x^1_1, x^2_1, x^3_1, x^4_1, x^5_1, x^6_1, x^7_1, x^8_1$) in eight distinct panels.

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_2}$ for initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^1_2)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_2}$ for initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^2_2)$]

The figures illustrate the temporal evolution of the wavepacket associated with the upper block of the Hamiltonian, focusing on the $x_2$ direction, capturing the torsional degree of freedom arising from planar rotations of the pyridyl rings about the C-C bond within the molecular architecture (refer to <ref>). The temporal evolutions depicted in each subfigure is initialized at distinct grid points as delineated in the subfigure captions, with each time-traces projected onto specific grid points indicated in the panel titles.

<ref> illustrates the dynamics for the upper block of the block-diagonal Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$, along the $x_2$ direction. This encodes the torsional degree of freedom that acts here as a gating mode and regulates proton transfer (larger angles have larger donor-acceptor distances). The gating mode originates here from planar rotations of the pyridyl rings about the C-C bond within the molecular structure depicted in <ref>(a). Specifically, <ref>(a) displays the projection of the time-traces at distinct grid points ($x^1_2, x^2_2, x^3_2, x^4_2$) in four distinct panels. In these panels, the initial wavepacket is localized at grid point '1' and has a zero value at all other grid points, defined as $\Psi_0(x_2) = \delta(x_2-x^1_2)$. Similarly, <ref>(b) presents the time-traces when the initial wavepacket is $\Psi_0(x_2) = \delta(x_2-x^1_2)$. In both figures, solid yellow lines represent the dynamics simulated classically through the exact numerical solution to the Schrödinger equation, while the blue dots indicate outcomes derived from ion-trap dynamics.

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_2}$ for initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^6_2)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_2}$ for initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^7_2)$]

[Quantum (blue) and classical (yellow) time-traces for $\ket{x_2}$ for initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^8_2)$]
The time-evolution of the wavepacket linked with the lower block of the Hamiltonian, emphasizing the torsional degree of freedom along the $x_2$ direction within the molecular configuration illustrated in <ref>. The temporal evolution depicted in each subfigure is initiated from distinct grid points, as specified in the subfigure captions, Additionally, each time-trace is projected onto specific grid points as indicated in the panel titles.

The dynamics of the $4\times4$ lower block of the block diagonal Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$ is simulated on a two-qubit ion-trap system.

Specifically, <ref>(a) presents the projection of the time-traces at distinct grid points ($x^5_1, x^6_1, x^7_1, x^8_1$) in four distinct panels. Here, the initial wavepacket is localized at grid point '6' and has a zero value at all other grid points, that is, $\Psi_0(x_2) = \delta(x_2-x^6_2)$. Similarly, <ref>(b-c) depict the time-traces when the initial wavepackets are $\Psi_0(x_2) = \delta(x_2-x^7_2)$ and $\Psi_0(x_2) = \delta(x_2-x^8_2)$ respectively.

Simulation type     Initial Wavepacket     Wavepacket Error

4*$x_1$-upper block     $\delta(x_1-x_1^1)$     0.278

$\delta(x_1-x_1^2)$     0.200

$\delta(x_1-x_1^3)$     0.192

$\delta(x_1-x_1^4)$     0.264

4*$x_1$-lower block     $\delta(x_1-x_1^5)$     0.221

$\delta(x_1-x_1^6)$     0.162

$\delta(x_1-x_1^7)$     0.236

$\delta(x_1-x_1^8)$     0.232

2*$x_2$-upper block     $\delta(x_2-x_2^1)$     0.199

$\delta(x_2-x_2^2)$     0.233

3*$x_2$-lower block     $\delta(x_2-x_2^6)$     0.226

$\delta(x_2-x_2^7)$     0.188

$\delta(x_2-x_2^8)$     0.222

2*$x_1$-full     $\delta(x_1-x_1^1)$     0.080

$\delta(x_1-x_1^2)$     0.090

$x_2$-full     $\delta(x_2-x_2^1)$     0.077

Wavepacket Error

The figure illustrates the quantum dynamics associated with the torsional degree of freedom ($x_2$), arising from planar rotations of the pyridyl rings around the C-C bond. The initial wavepacket is $\Psi_0(x_2) = \delta(x_2-x^1_2)$, and the projection onto individual grid points is presented separately in each panel.

The quantum dynamics simulation of the $8\times8$ full Hamiltonian, $\tilde{H}^{[2]}_{\gamma,\beta}$, has been executed on a 3-qubit ion-trap system. This simulation is aimed at determining the relative energy differences between the eigenstates derived from the Fourier transform of the dynamics simulations of the upper and lower blocks of the block-diagonal Hamiltonian. In <ref>, we present the quantum dynamics of the Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$ for the initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^1_2)$. It showcases the projection of the time-traces at distinct grid points ($x^1_2, x^2_2, x^3_2, x^4_2, x^5_2, x^6_2, x^7_2, x^8_2$) in eight distinct panels.
While the dynamics' trends are effectively conveyed in the figures presented in <ref>, we further quantify the numerical difference between the quantum and classical simulations as follows:

\begin{align}
\Delta\Psi = \sqrt{\frac{1}{NT}\int_0^{T}\dd{t}\int_0^{N} \dd{x_j} \|\psi_{\text{Classical}}(x_j;t)- \psi_{\text{Quantum}}(x_j;t)\|^2},
\label{Eq:T-Ave_Error}
\end{align}

Here, $\psi_{\text{Quantum}}$ represents the wavepacket at time, $t$, obtained from the ion-trap, while $\psi_{\text{Classical}}$ signifies the corresponding wavepacket derived from exact diagonalization results. The variables $T$ and $N$ denote the total propagation time and the size of the Hamiltonian involved in simulating each Hamiltonian block, respectively. These errors between $\psi_{\text{Quantum}}$ and $\psi_{\text{Classical}}$ for each simulation are summarized in <ref>.
As noted, the error is of the order of $10-20\%$ loss in probability and this is consistent with that seen in Ref. Sandia-1D-3Qubits. However, as will be seen below,  and as already noted in Ref. Sandia-1D-3Qubits, the spectral features are in greater agreement since the error in quantum computation appears to not be frequency dependent.

Hamiltonian      $\Delta \omega$ (THz)     Max. frequency (THz)

$x_1$-upper/lower block     $1.250$     $200.0$

$x_2$-upper/lower block     $2.125$     $340.4$

$x_1$-full Hamiltonian    $0.031$     $5.0$

$x_2$-full Hamiltonian    $0.445$     $71.4$

Fourier transformation characteristics.

§ COMPARISON OF SPECTRAL FEATURES OBTAINED FROM THE WAVEPACKET TIME-EVOLUTION

The figure depicts the experimentally determined relative energy separation among the four eigenstates associated with the upper block of the Hamiltonian, which corresponds to the vibrational degree of freedom involving proton transfer along the internuclear axis ($x_1$) connecting the two nitrogen atoms in protonated 2,2'-bipyridine. These separations are obtained through the Fourier transform of the time evolution of the upper block of the Hamiltonian, considering four distinct initial states as detailed in <ref>. Panel `a' showcases ${\cal P}_{\chi_0}(\omega)$, mentioned in <ref>, resulting from the quantum propagation of the upper block, as depicted in <ref>(a), with the initial wavepacket $\Psi_0^1=\delta(x_1-x_1^1)$ projected onto four different grid points $\left\{x_1^1,x_1^2,x_1^3,x_1^4\right\}$. Similarly, panels `b', `c', and `d' display the Fourier transforms derived from the quantum propagation (see <ref>(b-d)) of the upper block of the Hamiltonian along the $x_1$-direction, each initialized with wavepackets $\Psi_0^2=\delta(x_1-x_1^2)$, $\Psi_0^3=\delta(x_1-x_1^3)$, $\Psi_0^4=\delta(x_1-x_1^4)$, respectively.
Finally, panel 'e' illustrates the ${\cal P}(\omega)$, mentioned in <ref>, obtained from summing the quantities in panels `a', `b', `c', and `d'. Dashed gray lines and labels show predicted frequencies from exact diagonalization of the nuclear Hamiltonian.

The Fourier transform of the quantum dynamics shown in <ref>, is computed using the formalism outlined in <ref>. Specifically, we calculated ${\cal I}(\omega, x_1)$ from <ref> for each subfigure (that is each initial wavepacket) in the mentioned figures. Subsequently, the ${\cal P} (\omega)$ values for each of the initial wavepackets is cumulated and the results are presented in <ref>.

But before we delve into the details of the analysis, an important factor to consider while inspecting the Fourier transform figures is the range of vertical axes, ${\cal P}_{\chi_0} (\omega)$ and ${\cal P} (\omega)$. To understand this, we use the Cauchy-Schwarz inequality[Riesz and Sz.-Nagy, 1990] in Eq. (<ref>) to obtain $\|{\cal I}_{\chi_0} (\omega; x)\| < T_{\chi_0}$, where $T_{\chi_0}$ is the total length of simulation (see <ref>). As a result, $\| {{\cal P}_{\chi_0} (\omega)} \| < {T_{\chi_0}}^2$ and these inequalities, as we will see, are consistent with the Fourier transform figures here.

absolute square of the quantity in <ref>, ${\left\vert {\cal I} (\omega; x) \right\vert}^2$ obtained for each grid point ${x}$ is normalized and integrated over the grid space to yield:

\begin{align}
{\cal P} (\omega) = \int dx \; \frac{{\left\vert {\cal I} (\omega; x) \right\vert}^2} {{\cal I}_{\text{max}} (x)}
\label{Eq:integ-Density-timecorrelation-FT-4}
\end{align}

where ${\cal I}_{\text{max}} (x) = \max({\left\vert {\cal I} (\omega; x) \right\vert}^2)$.

Given that the dynamics simulation has been performed for different initial wavepackets as specified in <ref>, the Fourier transforms, more specifically ${\cal P} (\omega)$, have been calculated for each wavepacket dynamics initiated with `r' different initial wavepackets and thus, cumulated to obtain:

\begin{align}
{\cal P}_{avg} (\omega) = \sum_i^r {\cal P}_{i} (\omega)/r
\label{Eq:integ-Density-timecorrelation-FT-5}
\end{align}

This analysis aims to extract the relative energy differences between the corresponding eigenstates. Subsequently, these energy differences are utilized to construct the energy spectrum of the molecule as illustrated in <ref> and in Figure 8 of the main paper. The resulting energy spectrum is then compared with the energies obtained through the exact diagonalization of the nuclear Hamiltonian. The frequency spacing, $\Delta\omega$, and the maximum frequency for each computed Fourier transform are presented in <ref>.

The Fourier transform of the quantum dynamics depicted in <ref>, associated with the upper block of the block-diagonal Hamiltonian $\tilde{H}^{[1]}_{\gamma,\beta}$, is computed and presented in <ref>. This analysis aims to extract the relative energy differences between the respective four eigenstates associated with it. Specifically, the Fourier transform of each panel in <ref>(a) is conducted to obtain ${\cal I}(\omega, x_1)$ and integrated over grid space to obtain ${\cal P}_{\chi_{0}} (\omega)$
and the resultant Fourier transform shown in panel `a' in <ref>. The resultant Fourier transform corresponds to the upper block of the Hamiltonian $\tilde{H}^{[1]}_{\gamma,\beta}$, with an initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^1_1)$. Similarly, panels `b', `c', and `d' show the resultant Fourier transform obtained from
each subfigure in <ref>(b-d), respectively
associated with the corresponding distinct initial wavepackets summarized in <ref> under the `$x_1$-upper block' row. Finally, panel 'e' illustrates the cumulative
Fourier transforms,
obtained from panels `a', `b', `c', and `d'. Dashed gray lines and labels indicate predicted frequencies from the exact diagonalization of the nuclear Hamiltonian and the dots represent the data obtained from the quantum computation. Visually the agreement looks near perfect and quantitative analysis of errors in eigenvalues obtained from this analysis is presented later.

The figure illustrates the experimentally determined relative energy separation among the four eigenstates associated with the lower block of the Hamiltonian, corresponding to the vibrational degree of freedom involving proton transfer along the internuclear axis ($x_1$) connecting the two nitrogen atoms in protonated 2,2'-bipyridine. Panel `a' showcases ${\cal P}_{\chi_{0}} (\omega)$, as defined in <ref>, derived from the Fourier transforms ${\cal I} (\omega, x_1)$ (refer to <ref>) of the individual panels in <ref>(a), with the initial wavepacket $\Psi_0^1=\delta(x_1-x_1^5)$. Similarly, panels `b', `c', and `d' display the Fourier transforms obtained from the quantum propagation (see <ref>(b-d)) of the lower block of the Hamiltonian along the $x_1$-direction, each initialized with wavepackets $\Psi_0^2=\delta(x_1-x_1^6)$, $\Psi_0^3=\delta(x_1-x_1^7)$, $\Psi_0^4=\delta(x_1-x_1^8)$, respectively. These panels present  ${\cal P}_{\chi_{0}} (\omega)$ obtained from ${\left\vert {\cal I} (\omega; x_1) \right\vert}^2$ summed over the four grid points $\left\{x_1^5,x_1^6,x_1^7,x_1^8\right\}$ as mentioned in <ref>. Finally, panel e' illustrates \({\cal P}(\omega)\), as defined in <ref>, obtained by summing the quantities in panels a', b', c', and `d'. Dashed gray lines and labels indicate the predicted frequencies from the exact diagonalization of the nuclear Hamiltonian.

Similarly, the Fourier transform of the quantum dynamics portrayed in <ref>, linked with the lower block of the block-diagonal Hamiltonian $\tilde{H}^{[1]}_{\gamma,\beta}$, is computed and illustrated in <ref>. This examination aims to discern the relative energy differences between the respective four eigenstates associated with it. Specifically, the Fourier transform of each panel in <ref>(a) is performed to obtain ${\cal I}(\omega, x)$ and then integrated over grid space to obtain ${\cal P}_{\chi_{0}} (\omega)$.

The resultant Fourier transform shown in panel `a' of <ref>. This resultant Fourier transform corresponds to the lower block of the Hamiltonian $\tilde{H}^{[1]}_{\gamma,\beta}$, with an initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^5_1)$. Similarly, panels `b', `c', and `d' display the resultant Fourier transform obtained by integrating ${\cal I}(\omega, x_1)$
for each subfigure in <ref>(b-d), respectively, associated with the corresponding distinct initial wavepackets summarized in <ref> under the `$x_1$-lower block' row. Finally, panel 'e' showcases the accumulation of Fourier transforms from all initial wavepackets
obtained from panels `a', `b', `c', and `d'. Once again, dashed gray lines and labels indicate predicted frequencies from the exact diagonalization of the nuclear Hamiltonian and the dots represent the data obtained from the quantum computation.

The figure illustrates the experimentally determined relative energy separation among the eight eigenstates of the Hamiltonian, corresponding to the vibrational degree of freedom involving proton transfer along the internuclear axis ($x_1$) connecting the two nitrogen atoms in protonated 2,2'-bipyridine. Each of the time-traces shown in <ref>(a) with the initial wavepacket $\Psi_0^1=\delta(x_1-x_1^1)$ projected onto eight different grid points $\left\{x_1^1, x_1^2, x_1^3, x_1^4, x_1^5, x_1^6, x_1^7, x_1^8\right\}$

undergoes Fourier transformation to compute ${\cal I}(\omega, x_1)$. These transformations are then summed to obtain ${\cal P}_{\chi_{0}} (\omega)$ as described in <ref>, resulting in the energy separation presented in panel 'a' of this figure. Similarly, panel `b' displays the average Fourier transforms derived from the quantum propagation (see <ref>(b)) of the Hamiltonian along the $x_1$-direction initialized with wavepacket $\psi_0^2=\delta(x_1-x_1^2)$. Finally, panel `c' illustrates the sum of the Fourier transforms, ${\cal P}(\omega)$, obtained from panels `a' and `b'. Dashed gray lines and labels indicate predicted frequencies from exact diagonalization of the nuclear Hamiltonian.

The Fourier transforms obtained from the quantum dynamics illustrated in <ref> of the Hamiltonian $\tilde{H}^{[1]}_{\gamma,\beta}$, which represents the vibrational degree of freedom of the proton along the internuclear axis between the two nitrogen nuclei in <ref>, are depicted in <ref>. Each panel in <ref>(a) undergoes Fourier transformation to compute ${\cal I}(\omega, x_1)$, as mentioned in <ref>. This is further integrated to derive the resultant Fourier transform ${\cal P}_{\chi_{0}} (\omega)$.
Panel `a' of <ref> shows the resultant Fourier transform corresponding to the quantum dynamics of the Hamiltonian $\tilde{H}^{[1]}_{\gamma,\beta}$ with an initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^1_1)$. Similarly, panel `b' corresponds to the Fourier transform associated with the initial wavepacket $\Psi_0(x_1) = \delta(x_1-x^2_1)$. Panel `c' illustrates the accumulation of Fourier transforms for all initial wavepackets
obtained from
the Fourier transforms in panels a' and `b'. The purpose of the Fourier transform of the full Hamiltonian $\tilde{H}^{[1]}_{\gamma,\beta}$ is to compute the relative energy separations between the eigenstates associated with the upper and lower blocks of the Hamiltonian.

The figure illustrates the experimentally determined relative energy separation among the four eigenstates associated with the upper block of the Hamiltonian, corresponding to the torsional degree of freedom arising from planar rotations of the pyridyl rings about the C-C bond ($x_2$) in protonated 2,2'-bipyridine. These separations are obtained through the Fourier transform of the time evolution of the upper block of the Hamiltonian, considering two distinct initial states as detailed in <ref>. Specifically, panel `a' showcases the Fourier transform, ${\cal P}_{\chi_{0}} (\omega)$, resulting from the quantum propagation of the upper block, as depicted in <ref>(a), with the initial wavepacket $\Psi_0^1=\delta(x_2-x_2^1)$ projected onto four different grid points $\left\{x_2^1,x_2^2,x_2^3,x_2^4\right\}$. Panel `b' displays the Fourier transform derived from the quantum propagation (see <ref>(b)) of the upper block of the Hamiltonian along the $x_2$-direction, initialized with wavepackets $\Psi_0^2=\delta(x_2-x_2^2)$. Finally, panel `c' illustrates the sum of the Fourier transforms obtained from panels `a' and `b'. Dashed gray lines and labels indicate predicted frequencies from the exact diagonalization of the nuclear Hamiltonian.

The Fourier transform of the quantum dynamics depicted in <ref>, corresponding to the upper block of the block-diagonal Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$, is calculated and presented in <ref>. Each panel in <ref>(a) undergoes Fourier transformation to compute ${\cal I}(\omega, x_2)$ as described in <ref>. This is then integrated over the grid space to produce the resultant Fourier transform ${\cal P}_{\chi_{0}} (\omega)$, shown in panel `a' of <ref>. This transform corresponds to the upper block of the Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$ with the initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^1_2)$. Similarly, panel b' displays the Fourier transform ${\cal P}_{\chi_{0}} (\omega)$, obtained from the integration of Fourier transforms ${\cal I}(\omega, x_2)$ for each panel in <ref>(b), associated with the initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^2_2)$. Finally, panel c' illustrates the accumulated Fourier transform, ${\cal P}(\omega)$,
derived from the transforms shown in panels a' and b'.

The figure illustrates the experimentally determined relative energy separation among the four eigenstates associated with the lower block of the Hamiltonian, corresponding to the torsional degree of freedom arising from planar rotations of the pyridyl rings about the C-C bond ($x_2$) in protonated 2,2'-bipyridine. Panel `a' showcases the average Fourier transform resulting from the quantum propagation of the lower block, as depicted in <ref>(a), with the initial wavepacket $\Psi_0^1=\delta(x_2-x_2^6)$ projected onto four different grid points $\left\{x_2^5,x_2^6,x_2^7,x_2^8\right\}$. Similarly, panels `b' and `c' display the Fourier transform derived from the quantum propagation (see <ref>(b-c)) of the lower block of the Hamiltonian along the $x_2$-direction, initialized with wavepackets $\Psi_0^2=\delta(x_2-x_2^7)$ and $\Psi_0^3=\delta(x_2-x_2^8)$ respectively. Finally, panel `d' illustrates the sum of the Fourier transforms obtained from panels `a', `b', and `c'. Dashed gray lines and labels indicate predicted frequencies from the exact diagonalization of the nuclear Hamiltonian.

Similarly, the Fourier transform of the quantum dynamics depicted in <ref>, associated with the lower block of the block-diagonal Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$, is computed and presented in <ref>. Specifically, the Fourier transform of each panel in <ref>(a) is performed to calculate ${\cal I}(\omega, x_2)$ in <ref> and then integrated to obtain the resultant Fourier transform, ${\cal P}_{\chi_{0}} (\omega)$,  shown in panel `a' of <ref>. This resultant Fourier transform corresponds to the lower block of the Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$, with an initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^6_1)$. Similarly, panels `b' and `c' display the resultant Fourier transform obtained by integrating the quantity ${\cal I}(\omega, x_2)$ over grid space, computed for each panel in <ref>(b-c), respectively, associated with the corresponding distinct initial wavepackets summarized in <ref> under the `$x_2$-lower block' row. Finally, panel `d' showcases the sum of the Fourier transforms obtained from panels `a', `b', and `c'. The peaks are compared with the predicted frequencies from the exact diagonalization of the nuclear Hamiltonian and shown as gray lines and labels in the figure.

The figure illustrates the experimentally determined relative energy separation among the eight eigenstates of the Hamiltonian, corresponding to the torsional degree of freedom arising from planar rotations of the pyridyl rings about the C-C bond ($x_2$) in protonated 2,2'-bipyridine. Each of the time-traces shown in <ref> with the initial wavepacket $\Psi_0^1=\delta(x_2-x_2^1)$ projected onto eight different grid points $\left\{x_2^1, x_2^2, x_2^3, x_2^4, x_2^5, x_2^6, x_2^7, x_2^8\right\}$ is Fourier transformed and summed to produce the energy separation in the figure here. Dashed gray lines and labels indicate predicted frequencies from the exact diagonalization of the nuclear Hamiltonian.

The Fourier transforms derived from the quantum dynamics depicted in <ref> of the Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$, which represents the torsional degree of freedom of the proton arising from the planar rotation of the planes containing the pyridyl rings in <ref>, are showcased in <ref>. Each panel in <ref> undergoes Fourier transformation to calculate ${\cal I}(\omega, x_2)$ and subsequent integration over grid space to yield the resultant Fourier transform, ${\cal P}_{\chi_{0}} (\omega)$ in <ref>,  as shown in <ref>. This resultant Fourier transform corresponds to the quantum dynamics of the Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$, with an initial wavepacket $\Psi_0(x_2) = \delta(x_2-x^1_2)$. The Fourier transform of the full Hamiltonian $\tilde{H}^{[2]}_{\gamma,\beta}$ aims to compute the relative energy separations between the eigenstates associated with the upper and lower blocks of the Hamiltonian.

The figure presents the experimentally determined frequency and energy spectra of protonated 2,2'-bipyridine corresponding to the vibrational degree of freedom of transferring proton along the internuclear axis ($x_1$) joining the two nitrogen atoms. Due to the block diagonal nature of the Hamiltonian (see <ref>), each block of the Hamiltonian is simulated separately on two separate two-qubit ion-trap systems. The Fourier transform of the time evolution of the full-Hamiltonian allows for the determination of relative energy separations between the eigenstates corresponding to the respective blocks. The panel `a' displays the average Fourier transform obtained from the quantum propagation of the upper block with different initial wavepacket (see detailed plot in <ref>) , while panel `b' illustrates the average Fourier transform derived from the quantum propagation of the lower block of the Hamiltonian along the $x_1$-direction initialized with different wavepackets (see <ref>). The panel `c' shows the average Fourier transform derived from the quantum dynamics of the full-Hamiltonian with two distinct initial wavepackets, also shown in <ref>, providing the relative separation between the two sets of eigenstates corresponding to the lower and upper blocks. The <ref>d presents the extracted frequencies from panels `a' to `c', color-coded according to their parent spectrum. This presentation allows for the experimental determination of the relative energies of all eigenstates. Additionally, <ref>e compares the quantum-computed energy eigenstates of the nuclear Hamiltonian (depicted as blue dots) with the results of exact diagonalization (shown as dashed gray lines).

Finally, the Fourier transforms
obtained from the <ref> are utilized to get the energy spectrum of the molecule, as depicted in <ref>, corresponding to the vibrational degree of freedom associated with the transferring hydrogen along the internuclear axis connecting two nitrogen atoms. The panel `a' in <ref>, which is imported from <ref>, displays the cumulative Fourier transform
obtained from the quantum propagation of the upper block with different initial wavepacket, while panel `b' (see <ref> for detailed plot) illustrates the cumulative Fourier transform derived from the quantum propagation of the lower block of the Hamiltonian along the $x_1$-direction initialized with different wavepackets. The panel `c', also shown in <ref>, shows the cumulative Fourier transform derived from the quantum dynamics of the full-Hamiltonian with two distinct initial wavepackets providing the relative separation between the two sets of eigenstates corresponding to the lower and upper blocks. In <ref>d, the extracted frequencies from panels `a' to `c' are presented, color-coded according to their parent spectrum, facilitating the experimental determination of the relative energies of all eigenstates. Additionally, <ref>e compares the quantum hardware-computed energy eigenstates of the nuclear Hamiltonian (depicted as blue dots) with the results of exact diagonalization (shown as dashed gray lines).

The figure presents the experimentally derived frequency and energy spectra of protonated 2,2'-bipyridine, focusing on the torsional degree of freedom (represented as $x_2$) associated with planar rotations of the pyridyl rings about the C-C bond. Leveraging the block diagonal structure of the Hamiltonian, each block is individually simulated using separate two-qubit ion-trap systems. The panel `a' illustrates the average Fourier transform resulting from quantum propagation of the upper block with 4 distict initial wavepackets (see <ref>), while panel `b' (see <ref> for more details) depicts the average Fourier transform obtained from the quantum propagation of the lower block along the $x_2$-direction starting with different initial wavepackets. Panel `c' exhibits the Fourier transform derived from the quantum dynamics of the full-Hamiltonian, elucidating the relative separation between the eigenstates of the lower and upper blocks. Furthermore, in <ref>d, frequencies extracted from panels `a' to `c' are color-coded by their respective spectra, facilitating the experimental determination of relative energies for all eigenstates. Additionally, <ref>e presents a comparison between the quantum-computed energy eigenstates of the nuclear Hamiltonian (represented by blue dots) and exact diagonalization results (depicted by dashed gray lines).

The energy spectra from Figures <ref>e and <ref>e are combined to form the full energy spectra of protonated 2,2'-bipyridine, depicted as blue dots in bottom panel. The bottom panel shows first thirty energy levels for simplicity. These spectra are compared to the exact-diagonalization results of the nuclear Hamiltonian, shown as gray lines. The bottom panel shows in the agreement in units of THz, whereas the top panel illustrates the absolute energy difference between the quantum-computed energies and the exact-diagonalization results in units of $cm^{-1}$. The root-mean square error between the quantum-computed energies and the exact-diagonalization results is $41.66\;\unit{cm^{-1}}$ for first thirty energy levels.

Similarly, the Fourier transforms
obtained from <ref> are utilized to derive the energy spectrum of the molecule, as depicted in <ref>. This figure encapsulates the experimentally derived frequency and energy spectra of protonated 2,2'-bipyridine, with a focus on the torsional degree of freedom ($x_2$) associated with planar rotations of the pyridyl rings about the C-C bond. Employing the block diagonal structure of the Hamiltonian, each block undergoes individual simulation using separate two-qubit ion-trap systems. Panel `a', also featured in <ref>, illustrates the cumulative Fourier transform
resulting from the quantum propagation of the upper block with 4 distinct initial wavepackets. Meanwhile, panel `b' (referenced in <ref>) portrays the cumulative Fourier transform obtained from the quantum propagation of the lower block along the $x_2$-direction, initiated with different wavepackets. Panel `c', imported from <ref>, showcases the Fourier transform derived from the quantum dynamics of the full Hamiltonian, providing insights into the relative separation between the eigenstates of the lower and upper blocks. Further insight is provided in <ref>d, where frequencies extracted from panels `a' to `c' are color-coded according to their respective spectra, enabling the experimental determination of relative energies for all eigenstates. Additionally, <ref>e offers a comparison between the quantum hardware-computed energy eigenstates of the nuclear Hamiltonian (depicted by blue dots) and exact diagonalization results (represented by dashed gray lines).

The energy spectra from Figures <ref>e and <ref>e are combined to obtain the vibrational eigenenergies corresponding to the two coupled modes for protonated 2,2'-bipyridine.
These states are compared with exact-diagonalization of the nuclear Hamiltonian with mean absolute error between the computational results is found to be 0.14 kcal/mol.

Next, the energy spectra in <ref> are combined to obtain the vibrational energy levels corresponding to the two coupled degrees of freedom.
The accuracy of these energies is compared with that of exact diagonalization results. As shown in <ref>, the mean absolute error is well within chemical accuracy.

[Wright et al., 2019]
authorK. Wright,
authorK. M. Beck,
authorS. Debnath,
authorJ. Amini,
authorY. Nam,
authorN. Grzesiak,
authorJ.-S. Chen,
authorN. Pisenti,
authorM. Chmielewski,
authorC. Collins,
et al., journalNature communications
volume10, pages5464 (year2019).

[Amazon Web Services, 2020]
authorAmazon Web Services,
titleAmazon Braket (year2020),

[Saha et al., 2021]
authorD. Saha,
authorS. S. Iyengar,
authorP. Richerme,
authorJ. M. Smith, and
authorA. Sabry, journalJ.
Chem. Theory Comput. volume17, pages6713
(year2021), ISSN issn1549-9618, 1549-9626,

[Orús, 2014]
authorR. Orús,
journalAnn. Physics volume349,
pages117  (year2014).

[Pathak, 2013]
authorA. Pathak,
titleElements of quantum computation and quantum
communication (publisherTaylor & Francis,

[Huang et al., 1994]
authorY. Huang,
authorD. J. Kouri,
authorM. Arnold,
authorI. Thomas L. Marchioro,
and authorD. K.
Hoffman, journalComput. Phys. Comm.
volume80, pages1 (year1994).

[Hoffman et al., 1991]
authorD. K. Hoffman,
authorN. Nayar,
authorO. A. Sharafeddin,
and authorD. J. Kouri,
journalJ. Phys. Chem. volume95,
pages8299 (year1991).

[Iyengar, 2006]
authorS. S. Iyengar,
journalTheo. Chem. Accts. volume116,
pages326 (year2006).

[Campbell, 1897]
authorJ. E. Campbell,
journalProc. Lond. Math. Soc. volumes1-29,
pages14 (year1897), ISSN issn0024-6115.

[Baker, 1901]
authorH. F. Baker,
journalProc. Lond. Math. Soc. volumes1-34,
pages347 (year1901), ISSN issn0024-6115,

[[Trotter, 1959](1959)]Trotter
authorM. F. Trotter,
journalProc. Am. Math. Soc. volume10,
pages545 (year1959).

[Suzuki, 1976]
authorM. Suzuki,
journalCommun. Math. Phys. volume51,
pages183 (year1976), ISSN issn1432-0916.

[Nelson, 1964]
authorE. Nelson, journalJ.
Math. Phys. volume5, pages332

[Iyengar and Jakowski, 2005]
authorS. S. Iyengar and
authorJ. Jakowski,
journalJ. Chem. Phys. volume122,
pages114105 (year2005).

[Kouri et al., 1995]
authorD. J. Kouri,
authorY. Huang, and
authorD. K. Hoffman,
journalPhys. Rev. Lett. volume75,
pages49 (year1995).

[Sumner and Iyengar, 2007]
authorI. Sumner and
authorS. S. Iyengar,
journalJ. Phys. Chem. A volume111,
pages10313 (year2007).

[Vendrell et al., 2007]
authorO. Vendrell,
authorF. Gatti, and
authorH.-D. Meyer,
journalAngew. Chemie Intl. Ed. volume46,
pages6918 (year2007).

[Diken et al., 2005]
authorE. G. Diken,
authorJ. M. Headrick,
authorJ. R. Roscioli,
authorJ. C. Bopp,
authorM. A. Johnson,
and authorA. B. McCoy,
journalThe Journal of Physical Chemistry A
volume109, pages1487 (year2005),
notepMID: 16833469.

[Press et al., 1992]
authorW. H. Press,
authorS. A. Teukolsky,
authorW. T. Vetterling,
and authorB. P.
Flannery, titleNumerical Recipes in C
(publisherCambridge University Press, New York,

[Richerme et al., 2023]
authorP. Richerme,
authorM. C. Revelle,
authorC. G. Yale,
authorD. Lobser,
authorA. D. Burch,
authorS. M. Clark,
authorD. Saha,
authorM. A. Lopez-Ruiz,
authorA. Dwivedi,
authorJ. M. Smith,
et al., journalThe Journal of Physical Chemistry
Letters volume14, pages7256

[Tucci, 2005]
authorR. R. Tucci,
journalarXiv preprint quant-ph/0507171  (year2005).

[Earnest et al., 2021]
authorN. Earnest,
authorC. Tornow, and
authorD. J. Egger,
journalPhysical Review Research volume3,
pages043088 (year2021).

[Shende et al., 2005]
authorV. V. Shende,
authorS. S. Bullock,
and authorI. L.
Markov, in booktitleProceedings of the 2005
Asia and South Pacific Design Automation Conference (year2005),