Converged quantum simulations of reactive solutes in superfluid helium: The Bochum perspective

Experimental techniques such as vibrational photodissociation spectroscopy using messenger tagging^1–4 or helium nanodroplet isolation (HENDI) experiments^5–10 have become increasingly popular in recent years for the investigation of complex molecules, clusters, and chemical reactions. Both classes of experimental techniques use helium either as a very gentle tagging agent or to fully embed the chemical system of interest in a superfluid helium environment. They have led to great insights into the properties of a wealth of different systems such as neutral^11–17 and protonated water clusters,^18,19 HCl⋯water aggregates,^20,21 reactive intermediates²² such as $C6H62+$ or metal clusters,^23–28 and amino acid peptides or even small proteins,^29–32 to name but a few. Moreover, such experiments turned out to be fundamentally important for understanding the bosonic nature of few to many interacting ⁴He atoms (as well as of para-H₂ molecules) and the so-called molecular superfluidity^33,34 via the superfluid response of bosonic helium around molecular impurities.^5,35–43

While spanning the full range from small sized helium aggregates to bulk like superfluid helium nanodroplets, the experimental techniques operate under unique temperature conditions on the order of a few kelvin in the case of photodissociation tagging spectroscopy or even at sub-kelvin temperatures³⁵ as established in HENDI experiments. This grants access to the system of interest close to its rovibrational ground state essentially in the absence of thermal fluctuations. For a detailed understanding of these systems under such ultra-cold conditions, quantum simulation techniques that take nuclear quantum effects (NQEs) into account are key to provide insights into the underlying physical and chemical processes from a fully atomistic point of view. They can therefore assist, complement, and guide experimental studies to shed light on many fundamental questions. Pioneering such quantum simulations were, for example, able to predict⁴⁴ that finite clusters consisting of only 64 ⁴He atoms already feature “manifestations of superfluid behavior” below about 2 K, which has been confirmed experimentally ten years later.³⁷ Many of these simulation techniques are based on the Feynman–Kac formulation of quantum statistical mechanics via discretized imaginary time path integrals.^45–48 This formalism is able to provide, in principle, the exact quantum structural and thermodynamic properties of chemically complex and even reactive many-body systems at finite temperatures, notably also including Bose–Einstein statistics⁴⁹ and thus superfluid response.⁴⁸ 

Indeed, atomistic quantum simulation approaches mainly based on finite-temperature path integral Monte Carlo (PIMC) techniques, but also using various quantum Monte Carlo approaches to access the quantum ground state, have greatly expanded our knowledge of single atoms and ions^50–61 as well as of chemically inert non-reactive molecules.^{38,41,42,62–86} These vary from systems interacting with a few up to many helium atoms with a particular focus on the quantum (micro-)solvation shell structure around these solute species and on the molecular superfluidity of that bosonic quantum solvent. Most of the pioneering studies have used force fields for the description of the molecular impurities in superfluid helium, where vibrational degrees of freedom were usually neglected, while rigid-body translations and rigid-body rotations were increasingly taken into account. There also exist dedicated hybrid schemes where exclusively inter-molecular degrees of freedom within molecular complexes or a single coordinate describing large-amplitude motion within a molecule are explicitly considered in the sampling together with helium,^65,73,86 while other degrees of freedom such as the intra-molecular vibrations are kept frozen. Yet, a large body of literature exists where the molecular solutes are even fully clamped in space in the sense of space-fixed classical point particles, contributing only a fixed external potential that acts on the helium atoms.

Although these PIMC approaches provide deep insights into the quantum solvation of chemically complex molecular systems at ultra-low temperatures, they allow one to study neither chemical reactions nor intra-molecular vibrations in the presence of bosonic helium, the latter being of overriding importance for flexible, floppy, or fluxional molecules. In contrast, quantum simulation techniques based on path integral molecular dynamics (PIMD)^87–89 and, in particular, ab initio path integral (AIPI) simulations^90–94 are ideally suited to study NQEs on fully flexible molecules, clusters, or molecular complexes including even covalent chemical reactions via the latter approach.^95–99 Such molecular dynamics based sampling, unfortunately, is not readily applicable to superfluid solutions since bosonic exchange cannot be efficiently included in PIMD for systems consisting of many bosonic particles despite promising recent progress to be addressed later.

The recently developed hybrid path integral molecular dynamics/bosonic path integral Monte Carlo (HPIMD/MC) method¹⁰⁰ combines above-mentioned advantages by efficiently sampling the bosonic solvent using path integral Monte Carlo while treating the fully flexible and reactive solute using ab initio path integral molecular dynamics. This has, for the first time, opened the door to investigate quantum solutes including their full reactivity inside superfluid helium environments. The pioneering simulations^100,101 of reactive HCl⋯water clusters with this hybrid technique relied on “on-the-fly” evaluations of the interactions within the solute species based on concurrent electronic structure calculations in the spirit of AIPI simulations.^90–94 However, this approach imposes rather high computational costs and is usually limited to a density functional theory based description of the electronic structure. In the end, this severely limits the convergence of these simulations in many respects, namely, in terms of (i) the statistical error of the computed properties, (ii) the systematic error due to the finite path integral discretization, (iii) the accuracy of the underlying electronic structure method for describing the solute itself, and (iv) the solute⋯helium interactions. In certain cases, computationally very efficient but specialized force fields can also be applied in the case of molecular fluxionality, such as developed for the highly flexible protonated methane molecule,¹⁰² for which, indeed, intriguing couplings between the large-amplitude motion and the bosonic exchange in the helium surrounding were disclosed only recently.¹⁰³ However, this approach has only poor capability for generalization to other interesting and larger chemical systems and, in addition, is limited in terms of accuracy by the underlying functional form, in particular when it comes to describing complex chemical reactions. An additional complication of essentially all simulation approaches for solute⋯helium complexes is that carefully developed, system-specific interaction potentials for the coupling of the solute and the He solvent are essential.^104–106 This has the disadvantage that for each system of interest, new interaction potentials with helium have to be computed and analytically fitted, which turns out to be a daunting task if molecular flexibility and chemical reactivity need to be taken into account.^105,106

In the following, we review in detail how essentially fully converged quantum simulations of fluxional and reactive solutes in superfluid helium can be performed due to most recent methodological advances while also presenting some unpublished developments that contribute to solving the above-mentioned limitations (i)–(iv). For this purpose, we first introduce our HPIMD/MC method¹⁰⁰ and also allude to various other simulation approaches. Subsequently, we address the issue of converging path integral sampling at finite but ultra-low temperatures within our hybrid approach by showing how the challenges that arise from ultra-cold conditions^107,108 in conjunction with the bosonic quantum nature of helium can be met. Thereafter, we outline a general strategy that provides computationally efficient and largely automated access to accurate descriptions of both inter- and intra-molecular interactions. Our particular focus is on recent developments based on machine learning techniques that rely on artificial neural network representations of both solute⋯helium¹⁰⁹ and intra-solute¹¹⁰ full-dimensional potential energy surfaces (PES) at the level of essentially converged coupled cluster electronic structure calculations. We finally conclude with an outlook on the simulation of reactive impurities in bosonic quantum solvents and spell out remaining challenges.

In the following, we first briefly discuss different approaches to simulate molecular impurities in superfluid helium before reviewing the main aspects of our general hybrid strategy for reactive solutes in quantum solvents. Most finite temperature path integral approaches for the simulation of molecules in superfluid helium rely on path integral Monte Carlo (PIMC) sampling for which various algorithms have been specifically designed to describe the rotation of rigid molecular impurities in ⁴He, e.g., established in Refs. 34, 71, 77, 81, and 111–113. However, as already mentioned in the Introduction, PIMC and QMC approaches, in general, are not ideally suited for the efficient simulation of chemically complex molecular systems if molecular vibrations or even chemical reactions shall be included since a lot of the standard sampling moves would usually be rejected. Thus, algorithms introducing specific MC moves are needed to ensure efficient sampling in some of these cases, which have been devised repeatedly such as rigid-body DMC to sample inter-molecular interactions, while the interacting molecules themselves are treated as frozen entities.^114–118 Rigorous ground state QMC methods have also been tailored for the specific needs of simulating molecules in helium environments.^119–122 Although granting access to the correct bosonic ground state, these approaches prevent the analysis of finite temperature effects on properties such as the superfluid density.

Therefore, in order to sample the quantum partition function of a molecular system solvated by a bosonic solvent, it is advantageous to split the system into two coupled subsystems, which can be sampled individually using the most efficient method. This is the main idea behind the HPIMD/MC method introduced in Ref. 100, where the interested reader can find all details, while we focus here on a short summary of the fundamental concepts. Each subsystem is described with a sampling algorithm that is best suited for its individual challenges. Molecular impurities usually feature high flexibility and possibly reactivity and, thus, are best described by PIMD propagation. This offers great sampling advantages for chemically complex molecular systems that are governed by highly anisotropic and strong intra-molecular interactions, since the forces required in molecular dynamics propagation guide best the sampling of the partition function, while MC moves relying only on energy differences would be rejected most of the time (except in the limit of using tiny displacements where, however, sampling of configuration space becomes prohibitively slow or upon introducing system-specific tailored MC moves). Quantum solvents that are subject to relatively weak intermolecular interactions, such as ⁴He but also para-H₂, are efficiently sampled by PIMC, which additionally enables one to include their bosonic nature. These two methods, which will be described in detail in Sec. III, are combined using the rigorous so-called MDMC hybrid sampling scheme introduced originally by LaBerge and Tully for a distinctly different purpose, namely, for simulating classical fluids.¹²³ 

For the simulation of bosonic helium, one has to account for quantum-mechanical particle exchange, which is at the very heart of many quantum phenomena such as the Lambda transition from normal to superfluid helium that occurs at 2.17 K in bulk helium at ambient pressure. For ⁴He, Bose–Einstein statistics has to be obeyed, which means that the wavefunction and the corresponding density matrix need to be symmetric under exchange of any two atoms. This implies that in addition to the configuration space, the “permutation space” has to be sampled if quantum exchange statistics shall be enforced. In the path integral isomorphism of a classical system of P beads connected by harmonic springs, this means that originally closed paths (also called “ring polymers” in some literature) of different originally distinguishable atoms can now be connected to each other,⁴⁹ thus lifting the distinguishable nature of identical particles according to Maxwell–Boltzmann statistics and thereby establishing Bose–Einstein statistics. Such global (topological) changes are achieved in a non-continuous manner so that sampling via molecular dynamics is significantly hampered, whereas suitable global PIMC moves can be introduced to rigorously sample symmetric exchange permutations, thus generating superfluid helium at finite temperatures. For further background and details on the theory and implementation of bosonic path integrals, we refer to the classic review by Ceperley.⁴⁸ 

Two highly efficient sampling algorithms based on PIMC have been developed for the task of simultaneously sampling permutation and configuration space in an ergodic manner. We note only in passing that in the realm of PIMD sampling, there have been several attempts to also generalize PIMD to bosonic systems.^125–128 This also includes largely¹²⁹ unpublished efforts by two of us, which in the end were not useful for simulating superfluid helium in view of the strongly repulsive He⋯He interactions that result in a radial distribution function that is essentially zero up to a distance of about 2 Å, thereby strongly suppressing the particle overlap and thus exchange permutations. Very recently,¹³⁰ a PIMD method has been proposed that avoids explicit enumeration of permutations of bosons subject to a very weakly repulsive model interaction potential resulting in a pair correlation function that is still sizable even for zero distance. In stark contrast to PIMD, successful PIMC methods to sample symmetric permutations are known since long. The first bosonic PIMC algorithm, called the bisection algorithm,^48,131 relies on resampling sections of a ring polymer by a Lévy construction that guarantees acceptance of the move for a free particle. For weakly interacting particles, this provides a very effective move with rather high acceptance ratios. Permutation space sampling is incorporated by resampling more than one section of different ring polymers at the same time, where the interconnection is swapped. Note that this includes parts of imaginary time paths of different atoms, thus making these moves rather unfavorable in the general case. In order to increase the probability of accepting such a change of the permutation state, the free particle transition probability of pair permutations between neighboring atoms can be used as a bias in the selection of appropriate permutation changes.⁴⁸ This algorithm has been implemented quite some time ago by HF in the CP2k program package^132,133 and underlies our original ab initio HPIMD/MC method.^100,101

Alternatively, the so-called worm algorithm has been introduced for efficient sampling of bosonic systems.^134,135 This technique opens the cyclic ring polymers, which transforms the sampling of the diagonal density matrix into off-diagonal elements [cf. open move in Fig. 1(d)]. Afterward, the opened ring polymer, called worm, can connect with that of other atoms via similar MC moves as the ones used for the bisection algorithm, thus effectively changing the permutation state of the system [cf. swap move in Fig. 1(g)]. However, in contrast to the bisection algorithm, only individual imaginary time sections need to be sampled and accepted, thus increasing the overall acceptance probability of permutation changes. After closing [cf. Fig. 1(d)], which moves the sampling back to the diagonal density matrix elements, the system ends up in a different permutation state. Originally introduced in the grand canonical ensemble,^134,135 the algorithm can be reduced to a canonical version by a suitable combination of the original moves and omission of some others.¹³⁶ Our canonical version of the worm algorithm (as introduced in the  Appendix and Fig. 1) has recently been implemented by us in the CP2k program package.^132,133 As usual in MC sampling, specialized moves can be designed for the sake of improving efficiency, such as the centroid move [cf. Fig. 1(b)] and other moves to resample different parts of the ring polymer [cf. Figs. 1(c), 1(e), and 1(f)] in order to enhance permutation sampling also in the limit of small finite systems where a molecule is microsolvated by only a few helium atoms, which, however, fairly strongly interact with the solute.^103,137 The different moves used in our CP2k implementation of the canonical worm algorithm are described in detail in the  Appendix.

Due to the simple additive two-body nature of the He⋯He potential, the pair density approximation has been developed to efficiently converge the path integral discretization.^48,138 We note that the pair potential is not only a good approximation in the bulk phase of helium, but also in the helium microsolvation shell around even fairly strongly interacting molecules as demonstrated in Sec. III C of Ref. 106, where the three-body contribution to the interaction energy has been shown to not exceed one percent as estimated by explicit coupled cluster electronic structure calculations. The pair density approximation allows one to efficiently perform PIMC simulations of bosonic helium even at temperatures on the order of 1 K or less. In this approach, the many-body density matrix is expanded in one-body and two-body contributions to the many-body action, where the nontrivial two-body contribution to the total action is approximated by the sum over the two-body contributions of the exact pair action for all pairs of helium atoms. These two-body contributions can be efficiently calculated using the matrix squaring method combined with the partial wave expansion.⁴⁸ In practice, they are precomputed and tabulated on a regular grid, and further fitted to a few one-dimensional spline tables for quick evaluation during the simulation as implemented long ago by HF in the CP2k program package.^132,133 Overall, the number of replicas P required to converge the path integral describing the system can be significantly reduced by this approach (see Ref. 48 for more details).

Other work has been devoted to defining locally resolved superfluid response,^{63,67,140–142} where pragmatic proxies for superfluidity based on simple cutoffs for sufficiently long exchange cycle lengths were also proposed.^63,142,143 Clearly, superfluidity is a genuinely macroscopic concept, and any local perspective or analysis of finite bosonic systems introduces necessarily specific viewpoints beyond what is a physical observable in the strict sense. It is for these reasons wise to speak about “manifestations of superfluid behavior”⁴⁴ or to use the notion of “molecular superfluidity”^33,34 instead of simply calling significant bosonic exchange phenomena “superfluidity” in such a context.

As already explained, PIMD is used in our hybrid HPIMD/MC approach in order to sample the configuration space of the solute molecule that is embedded in a bosonic helium environment because MD is well suited for the solute subsystem and has several advantages over MC. In particular, the configuration space is sampled according to Newton’s equations of motion due to the classical isomorphism, and thus, no specialized moves need to be designed. In addition, thermostats to sample the canonical (NVT) ensemble, such as massive Nosé–Hoover chain thermostats,^92,150 or stochastic friction-based techniques, such as the Path Integral Langevin Equation (PILE)¹⁵¹ thermostat, further improve the sampling efficiency, as explicitly shown for doped helium clusters at low temperatures.¹⁵² Moreover, time dependent properties can, in principle, be approximately computed in the framework of PIMD using variants such as Centroid MD (CMD)¹⁵³ or Ring Polymer MD (RPMD)¹⁵⁴ if the temperature is not too low. Moreover, both methods are plagued with serious deficiencies (namely, the curvature problem in CMD and spurious resonances for RPMD) when it comes to computing vibrational densities of states or infrared spectra,^155,156 which can be partially cured in the case of RPMD.^157–160 

However, performing path integral simulations at very low, but finite, temperatures is a computational challenge as the number of beads P required to converge the discretized path integral numerically becomes prohibitively high. For this reason, most simulations of clusters or molecules embedded in ⁴He have been performed using either a completely frozen solute molecule or a rigid-rotor approximation for the description of the impurity. These simulations are thus completely neglecting the vibrational degrees of freedom, which requires the highest number of beads to converge. While this enormous simplification is usually a valid approximation for quasi-rigid molecules, such as SF₆, OCS, or CH₄, to name but a few such cases, it breaks down for more flexible molecules, such as the hydronium (H₃O⁺) cation subject to pronounced umbrella tunneling¹⁶¹ or the protonated methane (⁠$CH5+$⁠) molecule^103,137 in view of its hydrogen scrambling. In such cases, one needs to perform converged path integral simulations while accounting for all degrees of freedom, which is a nearly impossible task at ultra-low temperatures on the order of 1 K, if standard techniques based on the so-called primitive (lowest-order Trotter) approximation of the thermal path integral action are used. As we explicitly demonstrated by converged benchmark calculations using the PILE thermostat¹⁵¹ for some representative showcase molecules,^107,108 namely, $CH5+$⁠, CH₄ and $H5O2+$⁠, Trotter numbers of about P ≈ 10 000 are required to converge intra-molecular observables at about 1 K.

Various methods have been designed in order to accelerate the convergence of path integrals with respect to the number of discretization replica P. This includes the mentioned pair density approximation used to describe the He⋯He interactions, which is based on approximating the high temperature N-body density matrix as a product of pair density matrices.⁴⁸ While this approach is perfectly suited for the pairwise interaction of point-like particles used to describe the He⋯He interactions based on the tabulated or splined pair functions as explained earlier, it is not easily applicable for more complex cases, such as flexible or even reactive solute species and the corresponding solute⋯He interactions due to the curse of dimensionality. Fortunately, other techniques based on different approaches such as higher-order Trotter discretization schemes to go beyond the standard primitive approximation,^162–167 perturbation theory,¹⁶⁸ or colored noise thermostats^169,170 can also be used to effectively enhance the convergence of the path integral as a function of increasing discretization number P. We note in passing that in PIMD (or PIMC), fully converged values of different quantities, such as energies or structural properties, can, in general, be estimated by extrapolation to the P → ∞ limit^71,152,171 by making use of the asymptotic nature of the Trotter discretization.^172–174 This approach however requires us to perform calculations for different numbers of replica in order to extrapolate.

The striking computational advantage when using colored noise thermostats is illustrated with the help of Fig. 2, where the convergence with respect to the number of beads of the distribution functions of H–H distances for the quasi-rigid methane (CH₄) molecule and its highly fluxional cousin protonated methane (⁠$CH5+$⁠) is presented. Clearly, the use of colored noise thermostats, here PIQTB,¹⁷⁷ allows for a significant acceleration of the convergence with respect to P compared to standard PIMD simulations (using PILE¹⁵¹ in the present case). In particular, it should be noted that the distributions obtained using a small number of beads (such as P = 32) are already very close to the fully converged ones even in the difficult case of the highly fluxional protonated methane molecule. We note in passing that the PIQTB¹⁷⁷ and PIGLET^175,176 approaches yield quite similar convergence acceleration at ultra-low temperatures, as demonstrated recently in Ref. 108.

Having shown that colored noise thermostatting greatly improves the path integral convergence of structural properties by adding artificial fluctuations according to quantum harmonic oscillators, one might wonder if these do not distort the partly anisotropic quantum delocalization properties of nuclei within molecules. As demonstrated in Fig. 3 of Ref. 108, it is reassuring to find that this is not the case, although colored noise thermostats tend to overestimate the particle delocalization if too few beads are used, yet they convincingly reproduce anisotropic NQEs.

Thus, the use of colored noise thermostats allows one to perform converged PIMD simulations of molecular systems even down to ultra-low temperatures on the order of one K while accounting for the full flexibility of the molecule including notably large-amplitude motion and fluxionality.

As introduced in Secs. III A and III B, ultra-low temperatures on the order of 1 K require many thousands of replicas to yield a converged representation of the path integral in the primitive approximation. At the same time, the pair density approximation,^48,138 utilized for PIMC simulations of helium, provides the converged results using one to two orders of magnitude fewer beads, while Trotter convergence in PIMD simulations can be accelerated via colored noise thermostats by roughly the same factor of 10–100 as we have just seen. However, it is still possible that the solute part, which features usually very high frequency vibrations that are harder to converge, requires a higher discretization of the path integral than the helium environment when using the pair density approximation. If these distinctly different acceleration approaches shall be combined, as required in our HPIMD/MC scheme, one possibly has to couple different numbers of path integral replica for the helium environment (using the pair density approximation in PIMC) and the hosted impurity (using some colored noise thermostatting in PIMD). In recent years, different approaches for such a task have been developed, ranging from ring polymer contraction^190–192 to multiple-time stepping in imaginary time.¹⁹³ In the following, such ideas are transferred to our coupled HPIMD/MC approach in order to achieve a converged quantum description of the solute while coupling it efficiently and most accurately to the quantum solvent at ultra-low temperatures. Our approach relies mainly on the multiple-time stepping technique introduced in Ref. 193.

In order to derive a meaningful solute⋯solvent bead coupling, let us start by focusing on the Hamiltonian of the MC and MD subsystem individually. During the PIMD propagation of the solute in the fixed environment of the solvent, the solute is fully described by the solute Hamiltonian Ĥ_mol and its coupling via the interaction potential with the solvent Û_int so that Ĥ_PIMD = Ĥ_mol + Û_int. Here, it is important to realize that the He⋯solute interaction is relatively weak compared to the intramolecular solute interactions. Such energetically low-lying modes of a system exhibit less zero point energy than higher modes. As a consequence, it is expected that the He⋯solute interaction can be converged with fewer path integral beads than the solute. Having this in mind, one can derive a two-step Trotter splitting of the density matrix as shown in the following.

Consequently, the He⋯solute interaction is only applied to every P_strideth bead of the molecular solute; a scheme of the resulting coupling approach is presented in Fig. 3. As apparent from Eq. (11), these beads, which are directly coupled to the quantum solvent, experience this interaction P_stride-times stronger and distribute it via the harmonic couplings within the solute’s ring polymer to the adjacent beads. This also indicates the limitations of this approach: Once the instantaneous forces due to the He–solute interaction times P_stride are too large, the Trotter splitting will break down. As outlined above, the approximation is therefore only valid in the limit of sufficiently weak interactions and thus small U_int.

Let us next focus on the advantages of this coupling via striding compared to other approaches, especially at very low temperatures. In particular, ring polymer contraction and an averaging scheme are considered as other approaches to couple two systems with different numbers of beads. These methods are not derived here, but details for the contraction approach can, for example, be found in Ref. 190. The main idea of the first approach is to reduce the number of replicas of the path integral ring polymer by discarding the highest frequencies of the system in the normal mode representation. Averaging achieves a reduction in the number of beads by a simple geometric average of adjacent beads. In particular, the normal mode contraction approach has been shown to be a powerful tool to boost the performance of path integral simulations under ambient conditions^191,192 where rather low bead numbers suffice to ensure convergence.

However, for the purpose of studying systems at ultra-low temperatures using the HPIMD/MC method, multiple-time stepping in imaginary time has some advantages. To highlight this, it is useful to consider the PIMC subsystem in more detail. During the MC phase of our hybrid approach, the system is described by the Hamiltonian Ĥ_PIMC that is comprised of the solvent Hamiltonian Ĥ_He and the coupling via the interaction potential with the solute Û_int so that Ĥ_PIMC = Ĥ_He + Û_int. In this case, it is important that the solvent interacts with the solute in representative configurations of its quantum ensemble since it is kept frozen during MC sampling. However, the averaging and contraction based methods significantly alter the structural and delocalization properties of the solute. At the same time, multiple-time stepping in imaginary time does not modify the structural properties at all because it comes down to essentially skipping a certain amount of replica. In the path integral formalism, all replica are equivalent and thus provide the same structural information. The structural properties are of major interest for the correct coupling of the PIMC and PIMD subsystems within our hybrid sampling approach, but the delocalization of the PIMD systems also needs to be transferred correctly to the helium environment. Explicit tests by us have shown that multiple-time stepping in imaginary time is able to conserve the correct quantum delocalization even for large striding lengths P_stride, while both the averaging and contraction methods provide always smaller particle delocalization, in particular, at ultra-low temperatures. This clearly highlights that imaginary time striding provides the best representation of a reduced ring polymer for the purpose of coupling it with the helium environment under ultra-cold conditions.

After having confirmed that multiple-time stepping in imaginary time is a valid option for the coupling of two differently discretized subsystems, we will focus next on quantifying the remaining uncertainties of this approach. For this purpose, conventional PIMD techniques can be used to describe solute⋯He complexes with only a single He atom since this treatment is known to straightforwardly converge to the numerically exact result for P → ∞, although a very high number of beads are needed to provide such fully converged benchmarks at ultra-low temperatures. In direct comparison to this path integral sampling of the solute⋯He complex, our HPIMD/MC method is used to describe the same complex using a combination of all above-introduced techniques in terms of solute and solvent sampling, thus using a substantially lower number of replicas compared to the conventional simulations. For the sake of one-to-one comparison of the PIMD benchmark to the HPIMD/MC scheme, the identical additive intra-molecular solute and intermolecular solute⋯He coupled cluster potentials as introduced in Secs. IV B and IV C, respectively, are employed, although a global surface for the entire solute ·He_N complex could be easily generated for N = 1. This outlined validation procedure has been performed for the CH₅⁺ · He and CH₄ · He complexes at 1.67 K in order to assess fluxional and quasi-rigid solute molecules, respectively.

In Fig. 4, the H⋯He and O⋯He distance distribution functions of these two representative solute⋯He complexes are compared for the two simulation approaches. In order to investigate if the hybrid coupling between the PIMC and PIMD regimes is valid, these two distance distributions were chosen as the most important features. For both structural properties and both systems, our hybrid HPIMD/MC approach, using all introduced methods for efficient sampling at ultra-low temperatures, reproduces the conventional results with many more replica in the frame of standard PIMD simulations extremely well (still, some minor differences can be detected). The shoulder and main peak in the H⋯He distance distribution functions at smaller and larger distances, respectively, result from the different equivalent hydrogen atoms in the two systems with respect to the single He atom, as resolved in terms of the two sets of dashed lines (see the caption). These subtle features are nicely reproduced by the hybrid approach with only slight shifts in the distribution. Note that the corresponding average values differ only by about 1%, as visualized by the vertical lines.

In conclusion, this analysis demonstrates that our HPIMD/MC method for efficient treatment of the quantum convergence is able to provide an accurate description of solutes in helium for the given interaction potentials. While for a single He atom, it is, in principle, possible to resort to the global potential energy surface of the complex and thus to conventional PIMD techniques (such as PILE or massive Nosé–Hoover chain thermostatting), albeit at the price of employing very large Trotter numbers, this approach is impossible, in practice, for the simulation of larger solute⋯helium clusters including bosonic exchange involving many ⁴He atoms. The reviewed and introduced acceleration techniques therefore open the door to efficiently describe complexes solvated by superfluid helium and thus shed new light on these intriguing systems at ultra-low temperatures.

The van der Waals interactions between helium atoms, which are among the weakest intermolecular interactions known, lead to the truly unique properties of small helium clusters^194–197 and are deeply connected to the formation of the superfluid phase in liquid helium.⁴⁸ For the simulation of most of these properties, a pair potential fitted to highly accurate ab initio calculations¹⁹⁸ is sufficient to provide reliable results. This also enables the usage of the pair density matrix approximation,⁴⁸ described in Sec. III A, and therefore provides a computationally efficient and accurate approach for the simulation of many particle helium systems. Note that a lot of work has been devoted to the development of more sophisticated helium⋯helium potentials based on ab initio methods,^199–209 which include adiabatic, relativistic, or higher order terms in the many body expansion. However, for the most properties of interest in the present context, a sufficient agreement between computer simulation and experiment is already provided when using the long established pair potentials due to the work of Aziz and co-workers including the original empirical 1979 variant.²¹⁰ 

In order to illustrate the quality of the results obtained with the 1995 pair potential¹⁹⁸ that we usually employ in our simulations (combined with the pair density matrix approach using a Trotter discretization P such that P × T = 80 K and the bosonic PIMC sampling using our implementation of these techniques in the CP2k code, in particular the canonical worm sampling as defined in the  Appendix), different properties of simulated bulk helium are compared to experiment in Fig. 5. First, the specific heat C_V, computed by numerically differentiating the total energy with respect to the temperature, is compared to the experimental values as a function of temperature. The well-known divergence of the specific heat capacity of bulk ⁴He, called Lambda transition,²¹² when moving from the classical liquid phase to the superfluid phase of helium under ambient pressure conditions, is recovered at around 2.1–2.2 K, although finite size effects broaden the peak considerably compared to experiment even when using 512 ⁴He atoms (see Ref. 131 for the pioneering work and Ref. 145 for a finite size scaling study of the Lambda transition). Last but not least, the superfluid fraction f_s is also shown as a function of temperature in the bottom panel of Fig. 5. Again, the experimental behavior is well described in the simulations, and the phase transition close to the experimental value, 2.17 K, is indeed recovered, if the number of atoms is sufficiently large. We stress that we use a truncated octahedral supercell in our implementation in order to minimize the number of atoms located beyond the largest non-periodic spherical cutoff and also compute the winding number and thus the superfluid density accordingly, as detailed in Sec. II E of Ref. 100.

As demonstrated based on the physical observables collected in Fig. 5, the weak intermolecular interactions are faithfully described by the standard He⋯He pair potential used¹⁹⁸ and the bosonic quantum nature of ⁴He is also nicely reproduced by our PIMC simulations in agreement with the pioneering work.^131,144

What remains challenging is the description of the intra-molecular solute interactions, i.e., the full-dimensional global PES of the solvated species, as well as the solute–helium interactions, where small errors are sufficient to tip the balance with respect to the delicate helium–helium interactions. This demands both very accurate correlated electronic structure calculations and their accurate representation in terms of computationally efficient potentials that provide analytical forces in addition. In order to address these remaining challenges, our most recent developments based on machine learning potentials are presented in Secs. IV B and IV C.

In most theoretical studies of helium (micro-)solvation, the molecular solute species has been treated either as a classical point-like particle system^{62–64,67,83–85,213} or as a quantum rigid rotor,^{41,42,66,68–72,74,75,77–82} thus excluding vibrational degrees of freedom. These studies have led to a wealth of pioneering insights into quantum solvation as well as the general understanding of superfluidity and its manifestations at the microscopic scale. In particular, for quasi-rigid molecules, such as SF₆, OCS, or CO₂, to name but a few, it is expected to be an excellent approximation to restrict the description to rotational and translational degrees of freedom by freezing the molecular skeleton, thus neglecting any small-amplitude intra-molecular motion. At the same time, these approaches prevent a more detailed analysis of the influence of helium solvation on the solute molecule beyond rotational dynamics and are severe simplifications in the case of more flexible species, in particular so for those undergoing large-amplitude motion such as fluxional molecules. Other cases are those with soft inter-molecular degrees of freedom such as hydrogen bonded dimers where an interesting impact of helium microsolvation on tunneling splittings has been discovered,^65,73 albeit at 0 K and while keeping the interacting monomers frozen.

Including intra-solute motion and thus the internal flexibility of molecules, complexes, or clusters forces one to deal with full-dimensional potential energy surfaces of these species (which finally must also be coupled to helium via solute⋯He interactions for all possible solute configurations and, thus, at each point of the respective PES, as detailed in Sec. IV C). The situation gets even more complicated if the solute species can undergo chemical reactions. For small gas phase molecules, vibrational degrees of freedom have usually been included by means of specialized parameterizations of physical or purely mathematical functional forms to represent the PES.^214–218 This task can quickly become cumbersome and difficult, which easily hampers any routine exploration of chemically complex species including their reactions. In the realm of quantum simulations at finite temperatures, this problem has been solved long ago by introducing ab initio path integrals (AIPI)^90–94 where the electronic structure problem is solved on-the-fly while sampling the path integral using molecular dynamics techniques. This avoids the pre-parameterization of typically high-dimensional potential energy surfaces or interaction potentials, yet the price to pay is using computationally efficient electronic structure theory, which means resorting to density functional theory in practice. As is well known, the quality of these density-based methods very much depends on the specific case and on the choice of the functional, although excellent descriptions can be achieved as demonstrated repeatedly, for instance, for protonated methane^219–222 based on extensive benchmarking.

In the last few years, new approaches based on machine learning and thus on general functional forms have been established as an alternative to the more traditional approaches, as reviewed, for instance, in Refs. 223–227. Among many distinct approaches,^228–239 high-dimensional neural network potentials (NNPs)²²⁸ combined with atom centered symmetry functions²⁴⁰ have been shown to reliably and accurately reproduce electronic structure based reference calculations. Moreover, they allow for a high level of automation of the fitting process, as reviewed in detail in Refs. 224, 225, 241, and 242. This enables the rapid development of highly accurate potential energy surfaces for inter- and intra-molecular interactions.^243–248 Recently, such approaches have been demonstrated to even grant efficient access to custom-made high-dimensional interaction potentials at a quality level^{108–110,249,250} that is broadly considered to be the “gold standard” in quantum chemistry, namely, coupled-cluster theory using singles, doubles and perturbative triple excitations [CCSD(T)].²⁵¹ 

In order to illustrate in the present perspective the power of these modern machine-learned approaches as applied to quantum simulations, the first machine-learned CCSD(T) PIMD simulation¹⁰⁸ being published in 2018, our recently established adaptive automated fitting procedure¹¹⁰ (that has been employed earlier to generate the NNP of bare H₅O₂⁺, which has been used previously in Ref. 108) is applied here to four different systems. All four species are interesting targets for helium solvation, namely, the methane molecule, its fluxional cousin protonated methane, as well as the hydronium and Zundel cations since they represent a most challenging mix of quasi-rigidity, low-dimensional large-amplitude motion, and high-dimensional fluxionality. Their global PES are described by high-dimensional NNPs,^{224,225,241,242} which are trained to energies obtained from explicitly correlated coupled cluster calculations. These include single, double, and perturbative triple excitations in conjunction with the F12a correlation factor²⁵² and scaled triple contributions²⁵² in the aug-cc-pVTZ basis set [i.e., CCSD(T^*)-F12a/AVTZ], as implemented in the Molpro package.²⁵³ This reference method has been demonstrated to yield results very close to the complete basis set (CBS) limit²⁵² due to using this explicitly correlated triple-zeta basis.

The generation of our CCSD(T^*)-F12a/AVTZ quality NNPs is performed using an automated fitting procedure that identifies by itself the most representative training points¹¹⁰ by drawing on the established properties of neural network potentials.²²⁴ Key to our procedure is to perform the least number of computationally demanding CCSD(T^*)-F12a/AVTZ electronic structure calculations by an adaptive sampling strategy [while not using forces that are increasingly demanding to compute numerically or analytically for large N-atom systems when using consistently the CCSD(T^*)-F12a/AVTZ methodology]. Our strategy targets specifically these parts of the configuration space of the specific system, which are relevant for the envisaged usage of the NNP in computer simulations. We stress that this subspace is certainly different for low-temperature quantum simulations close to the quantum ground state compared to thermally activated dissociation or reactive scattering calculations, which must be carefully considered when using a certain NNP. In other words, no NNP should ever be used for purposes for which it has not been explicitly trained.

The main steps of our automated iterative procedure¹¹⁰ are the following:

1. A large number of configurations are generated by classical MD or PIMD simulations at different temperatures, which are initially performed on-the-fly using general purpose, efficient, and easy-to-handle DFT-based methods. Once a dynamically stable NNP becomes available, this NNP is used to continue these MD or PIMD simulations, thus greatly improving sampling efficiency.

2. Relevant configurations are selected specifically either by comparing two (or more) different NNPs, which feature large differences in regions that are still under-represented in the MD or PIMD sampling procedure, or by detecting situations where the configuration space of the existing training set is left. Note that for the very first iteration, a small number of randomly chosen configurations on the order of as few as 10–20 configurations are used to seed the procedure.

3. Reference energies (at the desired coupled cluster level) are computed exclusively for those specifically selected configurations (rather than using random or otherwise sampled configurations from the above-mentioned MD or PIMD runs), which are then added to the dataset.

4. New NNPs are fitted to the improved dataset, and the cycle starts over at the stage of using NNP MD or PIMD simulations.

In order to explicitly validate the accuracy of the resulting NNPs, short quantum PIMD trajectories at 1.67 K have been generated using the NNPs. Afterward, single-point coupled cluster calculations have been carried out along these trajectories at each and every individual PIMD sampling step. The resulting comparison of the energies along the trajectories is depicted in Fig. 6 for our four representative solute species, namely, CH₄, $CH5+$⁠, H₃O⁺, and $H5O2+$ from top to bottom. For all four systems, almost perfect agreement between the NNP and the explicit coupled cluster reference is observed, thus highlighting the quality of our approach (note that only a few red circles are added to visualize the CC lines that perfectly overlay the NNP lines on the scale set by the respective energy fluctuations). At the same time, the evaluation of the NNPs is many orders of magnitude faster than the coupled cluster reference, granting access to the converged electronic structure results at negligible cost.

Moreover, the individual NNPs can be developed very efficiently due to the automated procedure (so far successfully applied up to the protonated water tetramer in Ref. 110). This procedure iteratively improves the training set while keeping the demanding reference calculations to a minimum. This opens the door for routine exploration of medium sized molecules and molecular complexes under ultra-cold conditions including all their degrees of freedom as well as the quantum nature of the nuclei while employing an essentially converged description of the interactions.

Transcending this specific class of applications, we are confident that our adaptive automated NNP fitting procedure¹¹⁰ can also be used to generate global potential energy surfaces, including analytical forces and Hessians, for distinctly different purposes than (PI)MD, such as semiclassical wavepacket propagation or full quantum dynamics.

In order to investigate the impact of molecular motion on helium solvation and superfluidity (and vice versa), an accurate description of the intermolecular interactions is essential. Indeed, these solute⋯helium interactions are usually rather weak and can be comparable to the He⋯He interactions, and the competition between these two different interactions gives rise to complex solvation patterns that can change drastically with the number of ⁴He atoms. While wavefunction based ab initio methods such as coupled cluster theory are able to provide the required level of accuracy, they usually are computationally too demanding to be used “on-the-fly” in path integral simulations with notable exceptions.²⁵⁴ Indeed, most studies of doped helium clusters have thus used physically motivated analytical interaction potentials fitted to reproduce reference data usually obtained by accurate ab initio methods.^{41,62,64,66–68,70,72–76,78–85,103,137,213} However, this requires the design of specific interaction potentials for each system under study, which can quickly become a time consuming and complex task. This is especially complicated if one wants to rigorously reproduce the fine structure of helium solvation from small clusters composed of a few ⁴He atoms up to much larger ones, including possibly periodic bulk while still accounting for the full large-amplitude flexibility of the solute, including chemical reactions.^105,106,255 For this reason, previous studies have been focused on a relatively small set of dopant molecules for which accurate potentials are available, such as SF₆, OCS, or CO₂, and among others.¹⁰⁴ 

As for the above-mentioned intra-solute interactions, machine learning techniques such as neural networks can also be used to accurately describe the solute⋯helium interactions in an efficient way. Here, we use the same family of high dimensional NNPs^224,228 combined with atom centered symmetry functions²⁴⁰ to describe the atomic environment, as introduced before for the description of the solute potential energy surface. Because the typical values of the solute⋯helium interaction energies are usually completely different from the energy range associated with the inter-molecular interactions of the solute, it is beneficial to separate between inter- and intra-molecular interactions. Consequently, we use two different NNPs to describe these two types of interactions. Moreover, this separation is ideally suited for our hybrid HPIMD/MC sampling approach presented earlier, which relies on a different sampling strategy for the solute and the solvent. Now, the NNP is trained to reproduce the interaction energy between one helium atom and the solute molecule, thus accounting for two-body terms in a many-body expansion. The total interaction energy is then obtained in a pairwise fashion by adding all solute⋯helium two-body terms in the system. Thus, three-body and higher order terms are neglected in our approach. This has been shown to be an excellent approximation in the case of the $CH5+⋯$He interactions for which three-body terms have been shown by extrapolated coupled cluster calculations to account for less than 1% of the total interaction energy.¹⁰⁶ 

The interaction NNPs are obtained in a very similar fashion as just described for the intra-solute potential energy surfaces using the adaptive automated fitting procedure.¹¹⁰ However, small modifications of the procedure have been introduced¹⁰⁹ for the present task as follows. The initialization of the procedure is extended,¹⁰⁹ since it is not possible to perform reasonably accurate on-the-fly electronic structure based PIMD simulations of solute⋯helium complexes in view of the significant dispersion (also known as London or van der Waals) contributions to these interactions, which are difficult to faithfully include in DFT based AIPI simulations. First of all, the solute configurations are not sampled from previous MD or PIMD trajectories, but randomly extracted from the full dataset of reference coupled cluster calculations, which has been generated previously to fit the solute NNP. Second, a large number of positions for the helium atom are generated on randomly perturbed grids around these solute configurations as elaborated in the original article,¹⁰⁹ where an unbiased simple sampling MC procedure is applied on the grid to stochastically select a computationally manageable subset of points. At this stage, a short range cutoff is used in addition to avoid generating helium positions that are too close to the solute molecule, which would result in very high (repulsive) energies that could distort the NNP fit with no gains. Similarly, a long range cutoff is used to avoid generating the helium position too far away from the solute species for which the interaction energy approaches zero.

Akin to what is done in the case of the intra-solute interactions, as explained in Sec. IV B, two different NNPs are used to predict the interaction energies associated with the generated configurations and only the most relevant ones are selected based on the difference of predicted energies. However, the absolute values of the interaction energies in the repulsive (short range) region are significantly higher than the ones for helium positions far away from the molecule. Consequently, the absolute difference between the two NNP predictions, ΔE, is also significantly higher in the repulsive region, and these short range configurations are over-represented if only the configurations with the highest ΔE are selected. To circumvent this problem, a large number of configurations associated with the highest ΔE are extracted and then a Metropolis importance sampling MC procedure (instead of simple sampling MC) is used to select configurations with a probability distribution p(E) ∝ exp(−βE). The value of the β coefficient is not the inverse thermal energy 1/k_BT as usual but is also randomly chosen within meaningful bounds at each new iteration to automatically ensure a satisfactory sampling of helium configurations, which includes both short range and long range configurations. The β parameter bounds are adaptively chosen up to four times some average interaction strength of the solute species and the helium atom. Only after new configurations have been selected by this pre-screening technique, the associated reference interaction energies are computed using the standard “supermolecule” approach based on coupled cluster calculations. Here, the same electronic structure method as before for the solute molecule is used, namely, CCSD(T^*)-F12a/AVTZ, but with the counterpoise correction applied to correct for the basis set superposition error²⁵⁶ that remains even when computing interaction energies using explicitly correlated (F12) wavefunctions with a triple-zeta basis set.

This iterative fitting procedure, which is critically based on random sampling protocols to maximize ignorance as well as automation, has been shown to yield very accurate NNPs to describe the helium solvation of small protonated water clusters using datasets composed of only about 50 000 configurations for which the coupled cluster calculations need to be performed.¹⁰⁹ In the meantime, we have fitted several NNPs using this procedure for different dopants, namely, methane (CH₄), protonated methane (⁠$CH5+$⁠), hydronium (H₃O⁺), and the Zundel complex (⁠$H5O2+$⁠) with training and test RMSEs that are lower than 0.3 kJ mol⁻¹ even for the complex case of the highly flexible $CH5+$ molecule (while always using datasets containing about 50 000 configurations, which already includes the sampling of sufficiently many intra-molecular configurations of the solute species due to its vibrational displacements in view of the choice of the solute configurations as described above).

The ability of these NNPs to describe the solute⋯helium interactions is illustrated in Fig. 7, which shows helium densities in the three-dimensional Cartesian space as obtained from path integral simulations at T = 1.67 K for different doped helium clusters, namely, CH₄ · He₁₅, $CH5+⋅He8$⁠, H₃O⁺ · He₆, and $H5O2+⋅He10$ from top to bottom. For these particular simulations designed to benchmark the accuracy of the solute⋯helium interactions as given by the corresponding interaction NNP compared to the coupled cluster reference, the molecule is kept frozen in space (specified in the caption) and the solute⋯helium interaction energies needed to perform PIMC in the Cartesian continuum are interpolated from values that have been precomputed on rather dense grids using either the NNP or the reference method (being the counterpoise-corrected CCSD(T^*)-F12a/AVTZ calculations for these interactions). We note in passing that these grids contain typically 250 000 points at which the CC calculations have been carried out (which is one order of magnitude more than what is required to parameterize the interaction NNP), but solely for the purpose of carefully benchmarking the solute⋯helium NNPs “at work.” This procedure allows for visual comparison of helium densities obtained using the NNP and the coupled cluster method at the level of the relevant orientationally resolved spatial modulations due to the building up of quantum solvation shells. From the densities shown in Fig. 7, it is clear that the NNPs are able to very accurately reproduce the helium densities in all cases. Overall, this technique permits one to conduct path integral simulations of helium solvation at coupled cluster accuracy and without restriction on the number of ⁴He atoms.

In summary, the fully atomistic quantum simulation of molecules as well as molecular complexes and clusters solvated by superfluid helium has made significant steps forward due to many methodological improvements in the last few years. One of these developments is the hybrid (ab initio) path integral molecular dynamics/bosonic path integral Monte Carlo approach (HPIMD/MC)¹⁰⁰ that has been reviewed in detail in this perspective. This approach enables the description of fluxional and even reactive systems in bosonic helium while accounting for their quantum nature.¹⁰¹ The combination of (i) colored noise path integral thermostats^175–177 for the simulation of quantum solutes at ultra-low temperatures^107,108 with (ii) the pair density approximation and the sampling of the exchange permutations^{131,134,135,138} for the bosonic quantum solvent in the framework of (iii) the HPIMD/MC technique¹⁰⁰ provides an essentially converged quantum statistical description of dopants in superfluid helium. A striding approach to couple the solute with the solvent more accurately is introduced herein and shown to allow for the convergence of the molecular dynamics and Monte Carlo processes using different path integral methods and thus different levels of Trotter discretization to sample the entire path integral that represents the full solute–solvent system in the context of hybrid PIMD/MC sampling. Moreover, a canonical version of the so-called worm algorithm to efficiently sample bosonic permutations even for strongly interacting solute⋯helium systems containing only a few ⁴He atoms is presented herein and validated using as many as 512 bosonic atoms.

As recently demonstrated for the protonated cousin of methane, $CH5+$⁠, the HPIMD/MC technique can provide unique insights into the delicate couplings between large-amplitude motion leading to molecular fluxionality¹³⁷ and bosonic exchange statistics in the helium environment.¹⁰³ However, these recent studies were restricted to only a few helium atoms in the microsolvation shell of the solute due to the delicate nature of the solute⋯helium interactions that render the fitting of traditional pairwise interaction potentials difficult at best.^105,106

Recent developments based on neural network potentials (NNPs)^109,110 allow us to systematically overcome these previous shortcomings. By means of adaptive automated improvement, these NNPs provide a straightforward way for the investigation of complex molecular systems including clusters in superfluid helium. This approach gives access to interactions described at the essentially converged coupled cluster level, namely, CCSD(T^*)-F12a/aug-cc-pVTZ including the counterpoise correction when computing interaction energies. With these techniques now at hand, it will be possible to study the entire progression from small helium clusters to helium droplets to periodic bulk superfluid solutions that host fluxional and even reactive molecular species using a quantum simulation approach that is essentially converged at all levels. We hope that this will stimulate new experimental work toward applying helium messenger and nanodroplet isolation spectroscopies to those molecular species that are subject to large-amplitude motion, such as fluxional or floppy molecules, and to address cryochemical reactions beyond what has been achieved in the last decade.^20,21

In the context of rotating impurities in superfluid helium, as amply referred to in the previous text, it will be especially interesting to understand the salient couplings between the superfluid response and the dynamics of fluctuating charged molecules that feature much stronger solute⋯helium interactions compared to neutral species. Another possible direction for future work is the study of other bosonic quantum solvents that can also potentially feature “manifestations of superfluid behavior,” such as para-H₂ for which both pure^{140–142,257–259} and doped clusters^{33,139,260–263} have already been the subject of many insightful theoretical studies as comprehensively reviewed not long ago.³⁴ 

Let us now address some of the major computational challenges awaiting progress in the future. Accounting for the bosonic quantum nature of helium is already a very challenging task, but fermionic exchange can also become important at ultra-low temperatures, if light particles such as hydrogen atoms are present in the solute—not to mention the use of ³He as a fermionic quantum solvent.³⁷ Although it has recently been demonstrated to be of minor importance for the quantum structural properties specifically of protonated methane (see the detailed analysis in the supplementary material of Ref. 103), it remains unclear more generally speaking how the fermionic nature of the nuclei would influence different properties such as the entangled²⁶⁴ rovibrational dynamics being intricately coupled to the nuclear spin statistics within the solute species.^264–266 However, accounting for fermionic exchange in finite-temperature path integral simulations of sufficiently large and strongly interacting many-body systems is notoriously difficult and is even today one of the unsolved challenges in the field of computer simulations, although many attempts have been made for decades trying to overcome the difficulty. More importantly, for a more direct comparison with spectroscopic experiments, it would be beneficial to compute time dependent observables, which unfortunately is not directly accessible in imaginary time path integral simulations. Specifically in the context of rigid-body rotation of molecules in bosonic ⁴He or para-H₂ environments, several dedicated methods have been designed to extract spectroscopic information from imaginary time correlation functions.^{33,34,66,68,70,71,77,261,263} More direct approaches have also been designed in the framework of path integral molecular dynamics to generate approximate quantum dynamics at finite temperatures, thus providing access to time-correlation functions as noted above in passing. Yet, they seem to break down when reaching cryogenic temperatures due to different intrinsic problems. Clearly, computing the accurate quantum dynamical properties of fluxional and possibly reactive molecular solutes in bosonic solvents at cryogenic temperatures, in particular rotation-resolved infrared spectra as commonly recorded in helium tagging photodissociation and HENDI action spectroscopy, remains another grand challenge.

F.B. and C.S. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This research was partially supported by Grant No. MA 1547/19 to D.M. and is also part of the Cluster of Excellence “RESOLV”: funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy (Grant No. EXC 2033 – 390677874). In addition, we would like to thank Professor Joachim Werner for providing the auxiliary basis set functions of He needed to evaluate the F12a correlation factor within the Molpro package. C.S. acknowledges partial financial support from the Studienstiftung des Deutschen Volkes and the Verband der Chemischen Industrie. The computational resources were provided by HPC@ZEMOS, HPC-RESOLV, and BOVILAB@RUB.

There are no conflicts to declare.

In the following, the details of our implementation of the canonical worm algorithm are presented (see the main text for general background and panels B to G of Fig. 1 for its pictorial presentation). In order to formulate a canonical version of the worm algorithm that is well integrated into our hybrid HPIMD/MC simulation approach, some restrictions on the set of allowed MC moves are required. First of all, the moves obviously need to conserve the number of particles. Therefore, the “remove” and “insert” moves of the original grand canonical implementation^134,135 are entirely excluded. Second, as a result thereof, the chemical potential μ has to be removed from all acceptance probabilities of the different moves. To achieve this, the sets of moves from Refs. 134 and 135 are combined in such a way that the number of ring polymer springs is conserved when proposing a new configuration, thereby removing the chemical potential from the acceptance probability. The precise nature of the combined moves used in our canonical worm algorithm is described later in detail, and we refer to Ref. 136 for a previous version of canonical worm sampling.

Our canonical worm algorithm as described in this appendix and implemented in the CP2k program package^132,133 has been validated for bulk ⁴He by one-to-one comparison to data available in the literature^131,135,144 as compiled in Fig. 5.

The traditional centroid move [Fig. 1(b)] is the first update, which is additionally introduced. To perform a centroid update, an atom is selected randomly and the exchange cycle it belongs to (in open configurations also the head/tail bead) is translated by a random vector, where the components are drawn independently from an uniform distribution [−Δ_max, Δ_max]. As shown in Fig. 1(b), the change in action depends on all beads of all atoms from the selected permutation cycle, and this update is therefore computationally demanding. Note that the change in the kinetic energy contribution to the action is zero because the kinetic and potential energies can be evaluated independently.

The staging update is able to resample any imaginary time section but relies on two anchor beads that stay fixed during the update. As a consequence, it is not possible to move the head or tail bead with the staging update since these beads can exclusively be anchor beads. Therefore, specialized moves to resample the ends of the worm are introduced as follows. The move head and move tail update [Fig. 1(e)] are identical but act on the different ends of the worm. From the head or the tail end of the worm, l beads are removed like the first step in the close move and then resampled equivalently to the second step in the open move. This effectively moves the ends of the worm, which is accepted with a probability according to Eq. (A1).

The crawl move [Fig. 1(f)] is similar to the movement of the head or tail. However, instead of reconstructing the imaginary time section that has just been deleted, the other end of the worm is extended by the same number of beads. This shifts the worm in imaginary time by ±l depending on the direction of the crawl and effectively leads to the worm crawling through the configuration space. Although this move is not strictly required, it greatly helps in preventing situations where the worm’s head and tail are trapped in separate local minima, thus making it difficult to close the worm. These are situations that are easily encountered when using a rather small number of ⁴He atoms around solutes that have strongly anisotropic interactions with deep potential minima that are well separated in space.