Tunneling splittings observed in molecular rovibrational spectra are significant evidence for tunneling motion of hydrogen nuclei in water clusters. Accurate calculations of the splitting sizes from first principles require a combination of high-quality inter-atomic interactions and rigorous methods to treat the nuclei with quantum mechanics. Many theoretical efforts have been made in recent decades. This Perspective focuses on two path-integral based tunneling splitting methods whose computational cost scales well with the system size, namely, the ring-polymer instanton method and the path-integral molecular dynamics (PIMD) method. From a simple derivation, we show that the former is a semiclassical approximation to the latter, despite that the two methods are derived very differently. Currently, the PIMD method is considered to be an ideal route to rigorously compute the ground-state tunneling splitting, while the instanton method sacrifices some accuracy for a significantly smaller computational cost. An application scenario of such a quantitatively rigorous calculation is to test and calibrate the potential energy surfaces of molecular systems by spectroscopic accuracy. Recent progress in water clusters is reviewed, and the current challenges are discussed.

## I. INTRODUCTION

Water is one of the most ubiquitous and vital substances on Earth. It is widely involved in chemical reactions and the cycles of matter and heat. The seemingly simple hydrogen atoms and oxygen atoms form relatively strong intermolecular hydrogen bonds and a complex hydrogen bond (HB) network. Therefore, water has many abnormal properties and abundant phases.^{1–5} Conventionally, these properties and phases were understood using the ball-and-stick model, i.e., electrons were treated using quantum mechanics, which results in chemical bonds, and the nuclei were treated as classical particles, which form the skeleton of the structures. In this scenario, difficulties for theoretical descriptions of water arise from the need to explore the large HB network and get accurate HB interactions. In recent years, advances in computer simulation methods mean that quantum mechanical descriptions of the nuclei are also possible,^{6–11} opening a door for studies of new exciting quantum effects of water. Assisted by high-resolution experimental techniques, such nuclear quantum effects (NQEs) had been found to be highly relevant to a series of nontrivial properties.^{12–21} The tunneling motion of hydrogen atoms, which involves the breaking and reforming of hydrogen bonds, is a prominent example.^{21–24} However, the condensed nature of bulk, which results in broadening of its vibrational and rotational spectra, makes it impossible for one to obtain high-resolution experimental data to decipher details of these processes in its bulk phases.

Water clusters provide a good entry point to circumvent this problem^{25–27} as its small size allows us to capture the specific hydrogen movement patterns without being obscured by the broadening caused by a large amount of chaotic movement like in the bulk phase. The high-resolution rovibrational spectrum of gas-phase water clusters can give much more precise and detailed information, where the splittings of the rovibrational energy levels reveal different tunneling behaviors of hydrogen atoms.^{28–35} These splittings have been investigated for decades. Permutation-inversion symmetry is useful for assigning the splittings.^{36} At very low temperatures, the geometry of the clusters can only be rearranged by tunneling motions, which results in equivalent conformations in different wells being connected through comparatively low barriers. Splitting of the energy levels happens because of couplings of these degenerate states.^{37} A molecular symmetry group can be generated by these permutation-inversion operation elements corresponding to the geometry rearrangements. Then, the splitting pattern can be described by the irreducible representations of this group, and the splitting sizes are determined by the ease of tunneling motions.^{38}

With the splitting sizes calculated, one can compare them with the experimental observations and analyze the atomic details of the tunneling paths. The accuracy for the value of the splitting sizes depends on two factors, i.e., the quality of the potential energy surface (PES) and the reliability of the method with which tunneling is described.^{39} Assuming that a rigorous and practical method is available to address the latter factor, experimental values of the splitting sizes can then be used as new calibrations for the improvement of the water potential.^{40,41} This is not only essential for simulations of bulk water in a more general sense but also in the cases when NQEs are of interest.^{41–46} Therefore, it seems that we have reached at a point in time when both the experiment and the first-principles based PES are almost ready, but there is still a lack of a rigorous and practical method to describe the quantum mechanical tunneling of the nuclei. With such a method, one can first calculate the ground-state splitting to understand the tunneling behavior more deeply and then analyze the finer structure of the infrared spectrum. In turn, the experimental splitting sizes serve as calibrations to further improve the quality of the water potential.^{40,41} Positive feedback from the quality of the PES to the theoretical method in describing tunneling, and vice versa, will push theoretical descriptions of water systems to the final stage with dreamed accuracy.

Along this route, several methods have been developed to find the tunneling paths and to calculate the splitting sizes. To date, only dimer splittings can be calculated by directly solving the Schrödinger equation in full dimensions^{47} due to the unfavorable scaling of the cost with system size. Diffusion Monte Carlo (DMC) is a promising method^{48,49} but suffers from the fixed-node approximation especially in the case of larger clusters and high-excited levels. The Wentzel–Kramers–Brillouin (WKB) method combined with the group theory is convenient in describing the splitting pattern^{50–53} but lacks accuracy in predicting the splitting sizes due to the approximations made in the method and the *a priori* specified tunneling path. The ring-polymer instanton method can avoid the latter problem by optimizing a minimum-action path as the optimal tunneling path in full dimensionality.^{54–57} It is also applicable to large water clusters, thanks to the relatively inexpensive computational cost and favorable scaling with system size. Now, it is becoming a mainstream method and has been used to calculate the finer structure of the splitting spectrum.^{58,59} However, the steepest-descent approximation made in the instanton method makes it short of spectroscopic accuracy for systems with anharmonic soft modes.^{35,56} For example, previous studies have exhibited departures from the experimental value within a factor of two and little dependency on the different water potentials used.^{35,56} To cure this problem, a path-integral molecular dynamics (PIMD) method making use of density matrix elements has recently been developed by Althorpe and co-workers.^{60–63} Their pioneering work showed that it can give an accurate splitting size of the acceptor path in the water dimer compared to benchmark results and showed promising results for larger water clusters as well.^{63} Yet, there was still room for improvement, e.g., resolving the relatively large sampling errors and other issues. Our recent theoretical work has reproduced the experimental values of the torsional tunneling splitting in a water trimer, which further demonstrated the potential of the PIMD method in providing the most accurate tunneling splitting results to date for large water clusters.^{64}

This Perspective focuses on these two path-integral based methods, i.e., ring-polymer instanton and PIMD. In Sec. II, we review the two methods in detail and then in Sec. III show that a connection exists between them, i.e., the instanton is exactly a steepest-descent approximation of the PIMD method although they originate from the quantum partition function and density matrix element, respectively. In Sec. IV, we introduce the progress on accurate calculations of the tunneling splittings in water clusters from the dimer to hexamer one by one and discuss the challenges in deciphering more finer structures of the spectra. Section V concludes our Perspective.

## II. PATH-INTEGRAL BASED METHODS ON GROUND-STATE TUNNELING SPLITTING

### A. Ring-polymer instanton method

Here, we briefly introduce the tunneling splitting derived in the formalism of the ring-polymer instanton by Richardson and Althorpe, first in a simple two-degenerate-well system^{55} and then in a multiple-well system.^{56}

*E*

_{0}split by the mixing of the two states due to tunneling. In the low temperature limit, the contribution of higher energy levels can be neglected. This leads to

*Q*(

*β*) (2

*Q*

_{0}(

*β*)) is the partition function of the system with (without) the inclusion of tunneling,

*E*

_{0}±Δ/2 are the energy levels split by the tunneling, and

*β*is 1/

*k*

_{B}

*T*. This expression provides a link between the tunneling splitting and the partition functions of the system.

*Q*

_{0}(

*β*) is easy to evaluate. Applying the steepest-descent approximation, it corresponds to the situation of a collapsed ring-polymer in the bottom of one of the wells. This results in a simple harmonic vibrational partition function. With

*β*

_{N}=

*β*/

*N*,

*N*being the number of beads, one has

*ω*

_{s}is the harmonic frequency of the well. The formulas in this subsection are for one-dimensional systems but can be easily generalized to realistic systems.

Unlike *Q*_{0}, *Q* corresponds to the partition function of the entire double-well region with the barrier included due to the existence of tunneling. Again, applying the method of steepest descent to *Q*(*β*), one gets a saddle point in the ring-polymer space. This results in the instanton, i.e., a closed trajectory in the imaginary time that connects the two wells over a barrier with the least action. In the low temperature limit, with infinitely long imaginary time *βℏ*, the instanton trajectory spends most of its time staying in the wells and passes the barrier occasionally. Such a pass is called a “kink.” Each kink process requires a finite time, which can be ignored compared to the infinite total time. As will be seen in Sec. III, this finite time is of fixed size and can be defined as an imaginary “tunneling time.” The most simple trajectory involves two kinks to form a closed orbit. If no kink occurs, it will fall back to the situation of *Q*_{0}. In the ring-polymer formalism where a large number of *N* beads are used to uniformly discretize the total imaginary time *β*ℏ, most beads would stay in the wells and just a few on the barrier describing kink processes. One example is shown in Fig. 1. A finite number of *M* beads is needed for a “one kink” process, marked with hollow circles. The instanton trajectory in the figure contains four kinks, and each kink is exactly the same. In the gap between the kinks, it will stay for any length of imaginary time, which is equivalent to an arbitrary number of *N*_{1}, *N*_{2}, *N*_{3}, and *N*_{4} beads in well ** a** or

**, represented by the dashed lines. This leads to**

*b**N*

_{1}+

*N*

_{2}+

*N*

_{3}+

*N*

_{4}+ 4

*M*=

*N*, where

*N*≫

*M*.

*Q*

_{n}(

*β*) represents the contribution from a trajectory containing

*n*kinks. The prefactor is the number of ways of arranging

*n*kinks in a ring-polymer of

*N*beads. Given that

*N*is large, $CNn\u2243Nn/n!$. A trajectory can start at well

**or well**

*a***, which contributes a factor of two. Only even numbers of**

*b**n*are included in the sum since otherwise the trajectory would not be closed. In the

*β*→

*∞*limit, the kink processes will be fully separated by the long stay in wells, which allows us to treat these kinks as isolated and equivalent. Meanwhile,

*Q*

_{n}(

*β*) are to be divided by

*Q*

_{0}as shown in Eq. (1); the contribution of the beads in wells shall be offset. At last,

*Q*

_{n}(

*β*) is only related to

*n*(the number of kinks) and the contribution of one kink (donated as

*θ*) through

*θ*can be calculated by a

*linear*polymer composed of finite

*M*beads connecting the two wells from

**to**

*a***. The potential energy for this linear polymer is in a familiar form,**

*b**r*

_{i}to shorten the notation. A one-dimensional model with mass

*m*is considered for simplicity. With these treatments, the action reads

*T*→ 0. Each kink has such a mode since they are assumed isolated.

*Q*

_{n}, these zero-frequency modes should not be integrated out by the Gaussian integral. They donate a factor of $(\beta N\u210fSkink)n/2$ to the partition function, and the resulting

*Q*

_{n}reads

^{55}

**G**of the instanton. The prime indicates that the

*n*zero frequencies are excluded from the product and replaced by $(\beta N\u210fSkink)n/2$. This

*Q*

_{n}also reads

**G**

_{0}through

**J**is the mass-weighed Hessian of the linear polymer and

**J**

_{0}is the equivalent Hessian for the non-tunneling system with

*M*beads collapsed in a well.

*β*→

*∞*is always required in the derivation. Comparing with Eq. (1), one can get the final expression for tunneling splitting of a double-well system,

### B. Tunneling splitting in a multi-well system

*G*degenerate wells. Similar to the double-well system, mixing of these

*G*states allowed by tunneling makes the ground-state energy

*E*

_{0}split into a set of

*G*levels {

*E*

_{ν}}. Equation (1) generalizes to

*Q*(

*β*) in the low temperature limit can be evaluated by the sum of the contribution of various instanton trajectories with different numbers of kinks. However, this time, the orbit does not locate on a single barrier. Any two of the

*G*wells are connected by a unique tunneling path unless the barrier is very high to allow tunneling to occur. These tunneling paths can be marked as different kinds of kinks. Contribution of a closed instanton trajectory with

*n*kinks must include all possible arrangements starting from a well and returning to this well through

*n*kinks.

*Q*

_{n,ν}as the contribution of the closed trajectories of

*n*kinks that start and end at well

*ν*, the partition function can be written as

*Q*

_{n,ν}, a concept of connected graph,

*adjacency matrix*

**A**, is introduced. A kink connects two wells through a barrier. The element

*A*

_{λμ}is the number of such tunneling paths that connect well

*λ*and well

*μ*directly. In general,

*A*

_{λμ}is equal to 1 or 0; however, there are some special cases where there are multiple symmetrical paths connecting

*λ*and

*μ*.

^{56}An example of such path is shown in Fig. 3. The red circles represent the degenerate wells. They are connected by solid lines if tunneling paths exist between them; then, the corresponding

*A*

_{λμ}would be non-zero. $\u2211\kappa =1GA\lambda \kappa A\kappa \mu $ provides the number of connections composed by two kinks that start at

*λ*and end at

*μ*via an intermediate

*κ*. Similarly, $(An)\lambda \mu $ contains all the arrangement of

*n*kinks connecting

*λ*and

*μ*. With the help of

*N*

^{n}/

*n*!, all the closed

*n*-kink trajectories that start and end at the same well

*ν*can be found by $(An)\nu \nu $. The example in Fig. 3 shows a six-kink trajectory of the instanton. $(A6)\nu \nu $ can help to traverse all possible connections of the kinks.

*n*kinks. Here, it is assumed that there is no coupling between different kinds of kinks. One can take wells

*λ*and

*μ*as a double-well system, and

*θ*

_{λμ}has exactly the same form as

*θ*in Eq. (12). Then, the contribution of a trajectory is the product of corresponding

*θ*

_{λμ}. With the help of

*adjacency matrix*

**A**, one can enumerate all possible sequences to obtain

*Q*

_{n,ν}. A Hückel-type

*tunneling matrix*

**W**is defined for this purpose,

*h*

_{λμ}= −

*θ*

_{λμ}/

*β*

_{N}. Similar to Eq. (5), one has

**W**are the splittings

*E*

_{ν}−

*E*

_{0}.

### C. PIMD method

The multi-dimensional integral in the partition function can be solved analytically through semiclassical approximation,^{65–67} i.e., the method of steepest descent. The computational cost is normally low. Errors of harmonic approximation, however, exist. It is natural to think of using molecular dynamics sampling to avoid errors of this harmonic approximation.

^{68,69}following the derivation of Mátyus

*et al.*

^{60}For a molecular Hamiltonian, the coordinate representation of its density matrix has a clear definition. By inserting a complete set of eigenfunctions

*ψ*

_{n}, it has the form

*E*

_{1}−

*E*

_{0}is the tunneling splitting to be determined.

**is an arbitrary geometry as the molecular wavefunctions**

*r**ψ*

_{n}spread in the coordinate space. The splitting size Δ, however, is independent of

**and**

*r**β*. This means that one particular geometry is enough to determine the splitting in Eq. (21). The reactant

**=**

*r***is a convenient and practical choice as $P\u0302a=b$, where**

*a***and**

*a***are the reactant and product geometry, respectively. They are stationary points in the two wells. A concise relationship between density matrix elements and the splitting size exists,**

*b**ρ*(

**,**

*a***;**

*b**β*)/

*ρ*(

**,**

*a***;**

*a**β*) is labeled as

*I*(

*β*). The splitting size Δ can be calculated from

*I*(

*β*) with at least two different

*β*as $\beta \u0304$ is also an unknown quantity. The extension to a multi-well system is easy by making use of the permutation-inversion symmetry group.

^{61}Next, we shall concentrate on the calculation of

*I*(

*β*) using PIMD.

*ρ*(

**,**

*a***;**

*b**β*) is transformed into a multi-dimensional integral in a classical linear polymer phase space,

*M*beads are inserted between fixed end points

**and**

*a***, with the coordinates**

*b*

*r*_{1}, …,

*r*_{M}, each bead being a replica of the molecular system with

*f*dimensions. Different from the convention that

*N*is used to label the number of time slices in path-integral simulations, we use

*M*deliberately since in the instanton method

*M*is the number of time slices for the linear polymer. As will be shown later, the equivalence between the steepest-descent approximation of the path for

*ρ*(

**,**

*a***;**

*b**β*) in Eq. (24) and the linear polymer for one kink process in the instanton method is crucial to demonstrate that the instanton method is a steepest-descent approximation of the PIMD method. The imaginary time

*β*is cut into

*M*+ 1 parts,

*β*

_{M}=

*β*/(

*M*+ 1) and

*ω*

_{M}= 1/(

*β*

_{M}

*ℏ*). The

*p*

_{i,j}s are introduced by the

*M*×

*f*inverse Gauss integrals.

*μ*

_{i,j}can be arbitrarily chosen while a judicious choice can improve efficiency of sampling here as proposed in Ref. 62. Their contribution to the prefactor

*A*will later be offset by the same prefactor in

*ρ*(

**,**

*a***;**

*a**β*).

*I*(

*β*) can be written as

*M*→

*∞*. Performing molecular dynamics sampling directly on this formula is not practical because

**and**

*a***are far apart. It is natural then to think of dividing the space between**

*b***and**

*a***into sufficiently small segments. This strategy is equivalent to thermodynamic integration. A free energy function is defined as**

*b**λ*here denotes a reaction coordinate, ranging [0, 1], and geometry

**(**

*c**λ*) connects

**(0) =**

*c***and**

*a***(1) =**

*c***along a smooth path. Customarily, we take the instanton as the smooth path as it provides sensible orientations of the degenerate conformations by finding a minimum-action pathway. Then, one has**

*b***(**

*c**λ*) from

**to**

*a***,**

*b*_{λ}represents a thermodynamic ensemble average of the linear polymer with an end point fixed at

**(**

*c**λ*). The Hamiltonian $Ha,c(\lambda )M$ is obtained by replacing

**in Eq. (25) by**

*b***(**

*c**λ*).

**(**

*c**λ*) is the only variable in the Hamiltonian that explicitly includes

*λ*; therefore,

*V*(

**) =**

*a**V*(

**). Then, one has**

*b***(**

*c**λ*) depends on the specific path chosen from

**to**

*a***. The integral of**

*b**λ*from 0 to 1 can be implemented by numerical methods, such as Gauss–Legendre quadrature. The expectation values of

*r*_{M}at several

*λ*s are needed to evaluate Δ

*F*in Eq. (31), as determined by

## III. CONNECTION BETWEEN INSTANTON AND PIMD METHOD

In the case of tunneling splitting calculations, the relationship between instanton theory and the PIMD method may be unclear because the former is derived from the partition function in Eq. (1) while the latter starts from the density matrix element in Eq. (21). In this section, through a simple derivation, it is shown that the instanton method is also a semiclassical approximation of the PIMD method in the calculation of tunneling splitting.

*ρ*(

**,**

*a***;**

*b**β*), we can get a stationary point in the

*Mf*-dimensional linear polymer space corresponding to a fixed-end linear polymer connecting

**and**

*a***over the barrier that minimizes the Euclidean action. It will be exactly the same as the kink trajectory. This equivalence is crucial in establishing the connection between the instanton and PIMD methods. After integrating out the momentum**

*b*

*p*_{1}, …,

*p*_{M}, one has

*M*is assumed. The prefactor would be offset by the same term in

*ρ*(

**,**

*a***;**

*a**β*). $Ua,bM$ can be obtained by excluding the kinetic energy terms in Eq. (25). It is numerically equivalent to the potential energy in the instanton method in Eq. (6) with the only difference that fixed ends are used here. The Hessian matrix determined by $Ua,bM$, therefore, equals that of the instanton method.

**and**

*a***as shown in Fig. 2. In the instanton method, such permutation motion takes**

*b**N*steps to complete a cycle (

*N*is the number of time slices of the trajectory). In the density-matrix based PIMD simulation, however, the ends of the linear polymer are fixed at

**and**

*a***. If we assume that there are**

*b**l*beads on the barrier under the requirement that the action of the linear polymer is unchanged, only

*M*−

*l*steps of such permutation are allowed. Therefore, integration over this zero-frequency mode contributes a factor of $(M\u2212l)\beta M\u210fSkink$. This is the only difference with the instanton method, while other terms such as action

*S*

_{kink}and fluctuation Φ are all the same. Finally, we label the semiclassical approximation of the ratio of the density matrix elements

*I*(

*β*) as

*I*

_{SC},

**J**of the linear polymer that stands for

*ρ*(

**,**

*a***;**

*b**β*). $\omega k2$s are those of

**J**

_{0}from

*ρ*(

**,**

*a***;**

*a**β*). The prime indicates that the zero-frequency mode is omitted from the product.

*β*

_{M}, we further define the imaginary time experienced by these

*l*beads on the barrier as the imaginary tunneling time, $\beta \u0304inst=l\beta M$. Then,

*x*) ≃

*x*when

*x*is small (one can arrive at the same conclusion without making the tanh(

*x*) ≃

*x*approximation), we have

_{inst}in Eq. (14), we can get

## IV. TUNNELING SPLITTINGS IN WATER CLUSTERS

In the rovibrational spectra of gas-phase water clusters from microwave and far-infrared experiments, splittings of the rovibrational levels were detected, which revealed profoundly the wave nature and quantum dynamics of nuclei in water. These splittings, coming from the couplings between degenerate wells as introduced in Sec. I, revealed a more complex molecular symmetry^{36} than an ordinary system with a rigid structure. In the literature, the study of tunneling splitting was accompanied by but more challenging than the resolution of the structure of water clusters.^{25,70–72} From a structure in the ball-and-stick form, a complete permutation-inversion group^{36,73} could be constructed, which traverses all potential tunneling rearrangements. The tunneling motions that are important to the splittings can be mainly categorized as flipping (torsional tunneling), geared (anti-geared) interchange, and bifurcation tunneling.^{74} “Interchange” refers to when two water molecules in an H-bond change their roles as the donor and acceptor. Early calculations use *ab initio* electronic structures to estimate the barriers of different tunneling pathways, which were competent to identify the possible tunneling rearrangements responsible for the observed splitting. The estimations of splitting sizes were based on the WKB approximation or low-dimensional models.^{50,75,76} Consistent with intuition, breaking fewer H-bonds leads to a lower barrier and hence a larger splitting size. Flipping tunnelings do not break H-bonds and tend to produce the largest splittings. Although numerically rough, the splitting patterns could be explained successfully by the permutation-inversion group generated from the most possible rearrangements.

With the advancement in computing power, the instanton, as one of the semiclassical methods, has achieved more precise tunneling pathways and, hence, numerical results that are more convincing.^{56} The complex tunneling dynamics and splitting patterns in water clusters have been reviewed by Cvitaš and Richardson.^{74} The instanton not only represents the most appropriate path but also gives results accurate in the order of magnitude. The sources of errors are then gradually clarified. There is now an opportunity to build on the instanton to enable more rigorous calculations and decipher the finer structure of the vibration–rotation–tunneling (VRT) spectra. As one of the most successful attempts, the torsional tunneling splitting in water trimer was the first to be reproduced rigorously with an unprecedentedly high accuracy using the PIMD method.^{64} However, for other clusters, various problems still exist and block the way of accurate calculation, such as the rotation– and vibration–tunneling couplings. The current two theories, i.e., instanton and PIMD, focus on calculation of the ground-state tunneling splittings. Some attempts have been made to include these couplings as will be introduced later. We review recent progress and discuss the challenges by enumerating the water clusters from the dimer to hexamer.

### A. Rotation–tunneling coupling of the acceptor tunneling splitting in a water dimer

A water dimer is the smallest and earliest studied water cluster.^{28,77} It has one hydrogen bond with the two water molecules as donor and acceptor. A full 12-dimensional wavefunction method has been performed recently where the complete molecular Hamiltonian (including permutation-inversion symmetry) was solved by a kind of discrete variable representation method, and the whole VRT levels can be obtained.^{47} Using the ring-polymer instanton, up to five different tunneling pathways have been found between eight degenerate configurations.^{56} The corresponding five permutation symmetry operations can be generated from two basic ones, i.e., the acceptor tunneling and geared tunneling. A combination of these two operations can generate a doublet-of-triplet splitting.

Among them, the acceptor tunneling contributes a large doublet splitting as it involves no breaking of hydrogen bonds and thus has a small barrier.^{56,73,77,78} This tunneling is special in the sense that it is constituted by a tunneling path with barrier and a barrierless rotation operation, as demonstrated in Fig. 4. The swap of the two hydrogens of the acceptor [Fig. 4(a)] has to be coupled with the overall rotation around the inertia spindle [Fig. 4(b)] in order to finish the permutation [Fig. 4(c)]. A recent work has calculated the acceptor tunneling splitting using the PIMD method with the thermodynamics integration along an instanton path and obtained a value in excellent agreement with the variation method (diffusion Monte Carlo method).^{63}

We noticed that a symmetry factor of 2 is applied to the PIMD result in the previous work based on the reasoning that there are two equivalent instantons in the system.^{62} Yet, according to our experience, such symmetry factor may be inappropriate for the PIMD method since the ratio of density matrix elements does not depend on the path of TI. Hence, the splitting size could be erroneously magnified twice, depending on the ergodicity of the sampling in the PIMD simulations. Due to the neglect of the rotation–tunneling coupling in the current PIMD method, we suspect that the PIMD result might not agree spot on with the benchmark for the water dimer. An attempt has been made to evaluate this coupling.^{58} They proposed an extension to instanton theory, inspired by the derivation of Mátyus and Althorpe in a complete SO(3) symmetry,^{61} and qualitatively evaluated the acceptor splitting of the first several rotational-excited levels in the water dimer. More theoretical efforts on deciphering this finer structure are desired.

### B. Accurate calculation of the torsional tunneling splitting in a water trimer and the nonadjacent interactions

Quantitative reproduction of experimental trimer splittings has long been viewed as an important criterion for testing the quality of the 3-body interactions of a water potential.^{79–81} To achieve this, a rigorous full-dimensional treatment is required, and on top of this, the accuracy depends on the quality of the water potential. The torsional tunneling splitting in a water trimer is probably the first one in all clusters that have been calculated rigorously from first principles with spectroscopic accuracy.^{64} It marks significant progress of the water potential and demonstrates the high predictive power of the PIMD + TI method.

A water trimer has a ring structure with three hydrogen bonds in the oxygen plane and three unbonded hydrogens up or down the oxygen plane.^{29,82} Flipping of one of the unbonded hydrogens to the other side of the plane leads to an equivalent conformation while experiencing a very low barrier. Flipping the three unbonded hydrogens one by one results in six permutationally equivalent conformations [Fig. 5(a)]. The flip (or torsional) tunneling operation generates a C_{6} permutation-inversion symmetry group, which then gives a quartet splitting with a 1, 2, 2, 1 degeneracy. Due to the low barrier, the splitting size is large (up to tens of wavenumbers). The rotational excitations are several orders of magnitude smaller than the torsional splitting and display as a fine structure in the spectrum. Fortunately, little rotational–tunneling coupling and vibrational–tunneling coupling have been observed, therefore, the ground-state splitting could then be inferred from the transitions observed in the VRT spectra. Extensive and detailed experimental measurements gave complete level information and a precise splitting size.^{29,33,82–84}

From the theoretical perspective, however, with nine more dimensions than dimer, the water trimer is very challenging to be treated with a full 21-dimensional wavefunction method.^{85} Semiclassical methods on the other hand are not quantitative enough for the task. The PIMD method here plays an important role to enable a full-dimensional and rigorous calculation.

Our recent theoretical work accomplished a reproduction of the experimental values with an unprecedented accuracy of just one wavenumber from first principles^{64} [Figs. 5(b) and 5(c)]. The Hückel approximation (considering only adjacent interactions) was shown to be no longer valid, which means that sticking to the tunneling matrix in Sec. II B of the present Perspective would encounter inaccuracy in a quantitative calculation. To achieve high accuracy, the nonadjacent terms that represent a delocalization of the Wannier-style nuclear wave packets were included in the torsional Hamiltonian as performed in Ref. 64. These nonadjacent interactions contribute the shift of the split levels from an even 1:2:1 spacing. Such effect may also exist in other multi-well systems with low barriers, e.g., cyclic water pentamer.

### C. The tunneling splitting of vibrational-excited states in a water tetramer

The water tetramer is an ideal platform to test the four-body term in a many-body water potential. It has a ring structure like the water trimer and four unbonded hydrogens up or down the oxygen plane alternatively^{34} (Fig. 6). To construct an equivalent conformation, the four unbonded hydrogens must flip together to opposite sides of the plane. The tunneling motion is between only two symmetric wells, which is consistent with the doublet splitting observed in the VRT spectra.^{34,86} Apart from this, no other type of tunneling splitting has been observed. This doublet splitting is much smaller in size than the torsional one in the water trimer (2.26 GHz for $(H2O)4$ and 5.6 MHz for $(D2O)4$)^{34,86} due to the fact that although no hydrogen bonds are broken, all the four unbonded hydrogen atoms are involved in the tunneling. The tunneling mechanism is still undetermined due to complexity in the multi-proton transfer process.^{49,86} There are questions on whether the four hydrogen atoms go through a concerted tunneling pathway or a step-wise one and whether the intermediate state contributes or not. Fortunately, the PIMD + TI method should not be affected by uncertainties in the path since the result of TI is path-independent, and an ongoing work of our group is aimed to answer this question.

Another problem of water tetramer is the possible vibration–tunneling couplings. Because the splittings were observed in low-energy vibrational bands, there has been a possibility that the observed doublets are the sum or difference of splittings of the ground state and a vibrational-excited state^{86} as shown in Fig. 6. The vibration–tunneling couplings were considered to be weak because the vibration was supposed as a ring-deformation motion of the bonded hydrogens, which would have little influence on the flip motion.^{86} To achieve the goal of taking spectra data as the calibrator, methods that can accurately calculate the tunneling splittings of vibrational-excited states are desired. Early theoretical attempts were made in model systems.^{87} For molecular systems, an effort has been made recently to build a semiclassical method based on the instanton from a Hamilton–Jacobi formalism.^{59} Besides this, deciding whether the doublets are ground-state splittings or the sum still requires more evidence.

### D. Larger water clusters

The cyclic water pentamer was supposed to have a similar torsional tunneling motion like the water trimer.^{88,89} Flipping the five unbonded hydrogens in sequence could result in ten equivalent conformations that split the rovibrational state into six levels with 1, 2, 2, 2, 2, 1 degeneracy. In addition, the splitting size was considered to be large up to tens of wavenumber as the water trimer because of the low barrier of torsional motion. Several bands in the VRT spectra have been identified as transitions between six levels,^{90} but the complete information of all the splitting levels is still lacking to our knowledge. Theoretically, Cvitaš and Richardson computed the ground-state splittings using instanton theory^{91} with a tunneling matrix of 320 equivalent minima, including the torsional and bifurcation tunneling motions.^{88} This effort might be helpful in the continuing effort to assign the experimental spectrum. However, the nonadjacent interactions should be presented in the torsional tunneling splitting of the water pentamer as well,^{64} which were not included in the tunneling matrix of the instanton method and will impact the ratio of the spacings between the torsional splitting levels. Therefore, a prediction from a high-quality calculation using the PIMD + TI method is desirable, yet such calculation would require an *ab initio* water potential with highly accurate four-body interactions (and perhaps accurate 5-body interactions may also be required).

The water hexamer has several possible structures that were considered to have similar energy. During the identification of the structures, ticks of tunneling splitting have been found in the pure rotational spectra,^{35,71} showing a unique doublet-of-triplet splitting pattern. In the prism structure, the interchange tunneling pathways have been identified to be responsible for the splittings; the anti-geared mechanism breaks one H-bond, while the geared mechanism concertedly breaks two. Using ring-polymer instanton theory, Richardson *et al.* discovered that the geared mechanism plays a crucial role in producing the unique splitting pattern. The absolute size of the splittings predicted by instanton theory is larger than the experimental result due to the semi-quantitative nature of the method and the accuracy of the water potentials available at the time.^{35} Vaillant *et al.* made a first PIMD calculation and obtained splitting sizes closer to the experimental results with the same water potential.^{63} However, the sampling uncertainties remain large, and from our experience, fully converging the PIMD + TI for the hexamer is very challenging due to the relatively high barrier and the complexity of the tunneling pathway. Improving the sampling efficiency is a prior task for accurate calculations of the tunneling splittings in large-size water clusters. On the other hand, the experiment in Ref. 35 observed that the transitions between different pairs of rotational levels exhibit different splitting sizes. This indicates that rotation–tunneling coupling is non-negligible in the hexamer, and thus, the exact ground-state splitting size has not been precisely derived. The rotation–tunneling coupling may be the next key point to be considered in future theoretical calculations.

## V. DISCUSSION AND SUMMARY

Before concluding, we discuss a key issue, which is the impact of the water potential energy surface on the tunneling splittings. Previous studies have shown that the use of different water potential energy surfaces yields quantitatively different results. Specifically, in the case of the water trimer, a comparison was made between the MB-pol and “Our-pol” water potential energy surfaces.^{64} PIMD simulations using the former yielded *I*_{1}(*β* = 8000) = 0.58(2) and *h*_{1} = 26(2) cm^{−1} while using the latter produced *I*_{1}(*β* = 8000) = 0.520(3) and *h*_{1} = 22.1(2) cm^{−1}. The tunneling splitting value from “Our-pol” was much closer to the benchmark obtained from spectra data, which was 21.7 cm^{−1}. In a study of the water hexamer prism, Richardson *et al.* compared the tunneling splittings obtained using the HBB2-pol^{92} and MB-pol water potential energy surfaces.^{35} They found that the results obtained using the two water potentials could differ by 15%–30%. This disparity implies that the accuracy of the potential energy surfaces can be effectively evaluated by calculating tunneling splittings, provided that a rigorous treatment of the nucleus is achieved. Tunneling splittings can be used as the most stringent benchmarks to evaluate the accuracy of the potential energy surfaces.

In this Perspective, we reviewed the instanton and the PIMD methods on calculating ground-state tunneling splittings in water clusters. A clear connection is established between the two methods through a simple derivation. Instanton theory provides a computationally efficient way to predict tunneling splitting for larger systems, yet the splitting size it predicts is semi-quantitative due to the steepest-descent approximation made. The PIMD method improves the accuracy of calculation results by correcting this error, especially for systems with low barriers and strongly anharmonic modes.

The ability of the PIMD method to calculate rigorously and accurately the ground-state tunneling splitting has been illustrated in a water trimer. However, for other water clusters, before one can compare with the spectra data, more efforts are required. Each rovibrational level of an isolated molecular system should be split by the tunneling motions with the splitting size reflecting the ease of tunneling. The transitions observed in spectra are between rotational- or vibrational-excited levels. If the system exhibits a strong rotation– or vibration–tunneling coupling, it would be hard to figure out the splitting size of the ground state experimentally. We listed the current issues including the vibration–tunneling coupling in a water tetramer and rotation–tunneling coupling in a water dimer and hexamer prism. Some recent progress was summarized. The target of taking the spectra data of these clusters as calibrators of an *ab initio* water potential is challenging, and more theoretical efforts on solving rigorously the tunneling splittings of rovibrational-excited states are required. Besides, for large-size clusters and complex tunneling motions, there is still room for improving the sampling efficiency to enhance reliability and reduce computational cost. With the support of the two path-integral based methods and further developments in the methodology, we believe that accurate calculation of tunneling splittings in water clusters is possible, and we hope this also will facilitate constructions of water potentials with spectroscopic accuracy in the near future.

## ACKNOWLEDGMENTS

We thank Dr. Christophe Vaillant and Professor Stuart Althorpe for their guiding help on the PIMD method. The authors are supported by the National Science Foundation of China under Grant Nos. 11934003, 12234001, and 22288201; the National Basic Research Programs of China under Grand Nos. 2021YFA1400503 and 2022YFA1403500; the Beijing Natural Science Foundation under Grant No. Z200004; and the Strategic Priority Research Program of the Chinese Academy of Sciences under Grant No. XDB33010400. The computational resources were provided by the supercomputer center at Peking University, China.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Yu-Cheng Zhu**: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Writing – original draft (lead). **Shuo Yang**: Data curation (equal); Software (equal). **Jia-Xi Zeng**: Methodology (equal); Validation (equal). **Wei Fang**: Funding acquisition (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). **Ling Jiang**: Writing – review & editing (equal). **Dong H. Zhang**: Funding acquisition (equal); Project administration (equal). **Xin-Zheng Li**: Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

## REFERENCES

*Water: A Comprehensive Treatise*

_{2}O: A biography of water

*Ab initio*path integral molecular dynamics: Basic ideas

*ab initio*molecular dynamics

*Chemical Dynamics at Low Temperatures*

*ab initio*potential energy surface and the vibration-rotation-tunneling levels of (H

_{2}O)

_{2}and (D

_{2}O)

_{2}

*ab initio*and ‘hybrid’ potential energy surfaces, intramolecular vibrational energies, and classical IR spectrum of the water dimer

*ab initio*study of tunneling splittings in the water dimer

*ab initio*study of tunneling splittings in the water trimer

*The Whys of Subnuclear Physics*

*Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets*

^{3}He with the path-integral Monte Carlo method

*Molecular Spectroscopy and Quantum Dynamics*

*ab initio*calculations

*ab initio*and semiempirical potentials