where the sum is based on the bosonic modes. Here, the parameter  \(\epsilon\) denotes the energy difference between the ground-state energies of the states localized in the two wells without tunneling, and  \(\Delta_0\) indicates the matrix for tunneling between the two wells (Figure S2 (c)). The  \(\sigma_x, \sigma_y\) and \(\sigma_z\) are Pauli operators while \(\hat{b}_{\gamma}^{\dag}\) and \(\hat{b}_{\gamma}\) are the  creation and annihilation operators acting on the \(\gamma \ th\) bosonic mode of the heat bath and they obey the relationship \([ \hat{b}_{\gamma}, \hat{b}_{\delta}^{\dag}]= \delta_{\gamma \delta}\). The  \(c_{\gamma} \) is the coupling strength between the \(\gamma \ th\) bosonic mode and the wells. The dynamics of \(\sigma_z\) depends on the ratio of \(\hbar \Delta_0/\epsilon\). The eigenvalue of \(H_S\) is \(E_\pm=\pm \frac{1}{2} \sqrt{\epsilon^2 + (\hbar \Delta_0)^2 } \)\(\). Here,  \(|1\rangle \prod_{\gamma} |g_{\gamma1}\rangle \)and \(|0\rangle \prod_{\gamma}|g_{\gamma0}\rangle\)  denote the ground states in the separated right and left isolated local minimum, respectively. Note that \(\prod_{\gamma}|g_{\gamma i}\rangle\), i = 0 and 1, indicates the ground bosonic states in the well 0 and 1. To estimate the degree of population for each atomic eigenfunction \(|0\rangle\)and \(|1\rangle\) in the molecular eigenfunction \(\psi = e |0\rangle + f |1\rangle\), we adopt the ratio \((\frac{f}{e})\) , which satisfies the population in the ground state and the excited state as \((\frac{f}{e})_{ground}= \frac{\sqrt{\epsilon^2+(\hbar \Delta_0)^2}}{\hbar \Delta_0}\) and  \((\frac{f}{e})_{excited}= - \frac{\sqrt{\epsilon^2+(\hbar \Delta_0)^2}}{\hbar \Delta_0}\). Then, the excited state is \(\ \psi_+ = \cos(\frac{\theta}{2}) \ e^{- i \frac{\phi}{2}} \ |0\rangle \prod_{\gamma} \ |g_{\gamma0}\rangle + \sin(\frac{\theta}{2}) \ e^{i \frac{\phi}{2}} \ |1\rangle \prod_{\gamma} \ |g_{\gamma1}\rangle\) and the ground state is\(\ \psi_- = - \ \sin(\frac{\theta}{2}) \ e^{- i \frac{\phi}{2}} \ |0\rangle \prod_{\gamma} \ |g_{\gamma0}\rangle + \cos(\frac{\theta}{2}) \ e^{i \frac{\phi}{2}} \ |1\rangle \prod_{\gamma} \ |g_{\gamma1}\rangle\ \), where \(\tan \theta = \hbar\Delta_0 /\epsilon\) and \(\hbar\Delta_0 = |\hbar\Delta_0|e^{i\phi}\) or \(\phi=0\). By applying the spin system to the H-bond N-H...O in solution, the atomic orbitals of the N and O atoms are denoted as \(|0\rangle=\psi_N \) and \(|1\rangle=\psi_O\). Initially, the proton is located close to the N-atom. We here define the initial state as \(\Psi_S(0)=|0\rangle \prod_{\gamma}|g_{\gamma0}\rangle =[\cos(\frac{\theta}{2}) \ \psi_+ - \sin(\frac{\theta}{2}) \ \psi_-]\) at t = 0 and solve the Schrodinger equation of the time evolution of state  \(\Psi_S(t)= U_0(t,0) \ U_I(t,0) \ \Psi_S(0)\), where 'S' means Schrodinger picture and 'I' indicates the interaction picture. By using the exponential of Pauli identity \(\)\(\)\(\vec{\sigma}=(\sigma_x,\ \sigma_y,\ \sigma_z)\) and \(\exp(i \ \theta \ \vec{v} \cdot \vec{\sigma})= \cos(\theta) \ \sigma_0 + i \ \sin(\theta) \ \vec{v} \cdot \vec{\sigma}\), where \(\sigma_0\) = identity 2 x 2 matrix, we obtain
\(U_I(t,0)= \exp_+[- \frac{i}{\hbar} \int_0^t \ d\tau \ V_I(\tau)]\)
                    \(= \cos(\kappa) \ \sigma_0+ i \ \sin(\kappa) \ (\frac{A(t)}{\kappa} \sigma_x + \frac{B(t)}{\kappa} \sigma_y + \frac{C(t)}{\kappa} \sigma_z)\)
                    = evolution operator in the interaction picture
and
\(U_0(t,0)= exp_+[\frac{i}{\hbar} \ \int_0^t d\tau \ H_0(\tau)]\)
                     \(\)\(\)\(= \ e^{\frac{i}{\hbar} \ t \ \sum_{\gamma} \hbar \omega_\gamma \ \hat{b}_{\gamma}^{\dagger} \hat{b}_{\gamma} } \ [ \cos(\lambda t) \ \sigma_0 + i \ \sin(\lambda t) \ ( - \frac{\Delta_0}{2 \lambda} \ \sigma_x + \frac{\epsilon}{2 \ \lambda \ \hbar} \sigma_z ) ] \)
where \(\lambda = \frac{1}{2} \ \sqrt{\Delta_0^2 \ + \ (\frac{\epsilon}{\hbar})^2 }\). Then, we have the superposition state