Quantum microsolvation and size-resolved energetics of He-solvated H₂⁺ and D₂⁺ cations

Abstract

Helium nanoclusters provide a highly quantum solvent environment in which light ions experience pronounced zero-point motion and delicate competition between ion–solvent and solvent–solvent interactions. We report quantum diffusion Monte Carlo (DMC) calculations for protonated hydrogen cations embedded in helium clusters, H₂⁺@He_N and D₂⁺@He_N, quantifying ground-state energetics, evaporation energies, and microsolvation motifs up to medium-size N = 30 systems. Using CCSD(T)/CBS-quality interaction potentials augmented by short-range nonadditivity and correct long-range behavior, we analyze size-resolved stabilization patterns, aiming to identify possible shell-closure signatures governed by the balance of anisotropic cation–He forces and weak He–He binding. The resulting DMC probability densities reveal heterogeneous quantum microsolvation patterns, exhibiting both solidlike and liquidlike characteristics. Quantum effects are substantial, with the first two helium atoms found well localized, and as the cluster size increases, inner helium layers retain partial solidlike character, while pronounced quantum delocalization and He-exchange dominate the outer shells. As a consequence, no clear, well-defined shell-closure (“magic” N) could be identified, although the D₂⁺He_N clusters show a more compact, well-formed central ring around the [He–H–H–He]⁺ core for N ≥ 10 compared to the H₂⁺He_N ones. The computed evaporation energy trends rationalize experimentally observed ion–yield anomalies and delineate where semiclassical Feynman–Hibbs corrections deviate from exact ground-state sampling. Overall, this work establishes a benchmark quantum description for protonated He clusters and offers microscopic insight into the interplay of nuclear quantum effects, many-body interactions, and isotopic substitution in He-solvated ionic complexes, with implications for ion spectroscopy and low-temperature cluster physics.

---

1 Introduction

Helium droplets and finite He clusters are unique quantum media for isolating molecules and ions at ultralow temperatures, enabling precision spectroscopy and incisive studies of microsolvation.^1–17 Their defining features – weak He–He binding, near-frictionless environments, and large nuclear delocalization – create a sensitive balance between dopant–He attraction and the solvent's gentle structure. Charged impurities, in particular, can reorganize helium locally and promote distinct solvation shells whose rigidity and mobility vary with distance from the ion. Prior work has documented solidlike, localized first shells around strongly interacting cations, transitioning to liquidlike, delocalized second/third shells in larger clusters;^18–21 these structural motifs correlate with energetic “magic numbers” in ion–yield distributions.^18–30

From a computational perspective, a central challenge is to achieve a theoretical description that is both quantitatively accurate and fully quantum mechanical. On the potential-energy side, high-level ab initio data (e.g., CCSD(T)/CBS) have been condensed into analytical or machine-learned sum-of-potential models – combining ion–He and He–He contributions with selective three-body corrections – that accurately capture long-range interactions, as well as the anisotropy in both linear and T-shaped configurations.^31–43 On the nuclear motion side, fully quantum finite-temperature path-integral molecular dynamics (PIMD) has demonstrated pronounced zero-point effects, consolidated solidlike first shells and enhanced delocalization in outer shells,^19,20 as well as semiclassical approaches, such as the second-order Feynman–Hibbs (FH2).^20,21

Against this backdrop, protonated hydrogen cations embedded in He constitute compelling paradigms. Their light masses and anisotropic ion–He interactions strongly enhance nuclear quantum effects, while their simple electronic structure enables the construction of high-quality interaction potentials and facilitates stringent quantum tests. Moreover, experiments on ionic dopants in He droplets have revealed size-dependent abundance anomalies, commonly interpreted as signatures of shell completion and stability patterns that are extremely sensitive to minute energy differences between neighboring cluster sizes.

However, only a few high-resolution mass spectrometry studies of protonated He clusters (up to 50 atoms) and nanodroplets are available.^4,15–17 Early experiments detected only a very weak signal for He_NH₂⁺ clusters (3 ≤ N ≤ 10),¹⁵ while in another study, the limited mass resolution of doped ultracold superfluid He-nanodroplets prevented the identification of He_ND₂⁺ contributions.⁴ More recent measurements on hydrogen- and deuterium-doped helium nanodroplets ionized by electron impact reported weak anomalies in the ion–yield distributions for the He_ND₂⁺ and He_NH₂⁺ series at N = 13, 14, and 19.¹⁶ Owing to contamination effects, however, only the feature at N = 19 appears to provide reliable evidence of magic structures.

Thus, a theoretical approach that can deliver ground-state energetics with sub-cm⁻¹ accuracy for a given interaction model is essential for disentangling deficiencies in the potential energy surface (PES) from effects arising purely from nuclear quantum sampling. Here we employ diffusion Monte Carlo (DMC) to determine the zero-temperature quantum ground states of H₂⁺@He_N (D₂⁺@He_N) clusters over a broad size range. DMC provides essentially exact ground-state energies and probability densities for distinguishable nuclei given a specified PES, capturing the large zero-point delocalization that dominates He microsolvation. Building on CCSD(T)/CBS-quality ion–He interactions with carefully treated anisotropy and short-range nonadditivity, we (i) map total and single-atom evaporation energies as a function of N, (ii) identify and characterize magic-number candidates associated with shell completion and subsequent shell growth, and (iii) visualize quantum densities that reveal compact, snowballlike inner shells and increasingly liquidlike outer regions. Where appropriate, we contrast DMC trends with semiclassical FH2 results for these cation–He systems to clarify how effective-potential softening or finite temperature broadening compares with exact T = 0 quantum sampling.

The present DMC analysis complements prior quantum studies on similar He-solvated systems to show fixed inner-shell populations, suppressed exchange across the first-shell boundary, and progressive outer-shell delocalization. It also aligns with the broader picture that energetics alone may underresolve stability at medium sizes, requiring joint analysis of structure (radial/angle distributions and shell populations) and dynamical proxies to diagnose solidlike versus liquidlike behavior. By providing an exact ground-state benchmark for protonated He clusters, our DMC results help disentangle PES-driven effects from nuclear quantum contributions to stability patterns, thereby guiding the interpretation of experimental ion-yield anomalies and informing future many-body modeling strategies for such He-solvated ions.

2 Computational details, results and discussion

We used Cartesian (x, y, z) and spherical (R, θ, ϕ) coordinates to describe the He_NH₂⁺/D₂⁺ clusters. The origin of the coordinate systems is placed at the center of mass of the H₂⁺ or D₂⁺ molecule, and the z axis of the body-fixed (BF) frame is chosen to be parallel to the r vector joining the two H or D atoms. The R_i vectors are connecting the center of mass of the H₂⁺ or D₂⁺ molecule to each He atom, with i = 1 − N and θ_i is the polar angles of the vectors R_i with the z-axis, and ϕ_i the azimuthal angles of rotation around the z-axis (see upper-right panel in Fig. 1).

Fig. 1 Potential curves (in meV) of the VHeH2+21,42 and VHeHe44 (solid lines) as a function of the R and RHeHe distances at the indicated θ values (left panel) and minimum energy path of the VHeH2+ along θ (lower right panel) with r = re = 1.05886 Å. The coordinates used for the H2+HeN clusters are also shown (upper-right panel).

2.1 Interaction potential models

Following our previous work on He-doped clusters,^{18,21,37,45–49} the sum-of-potentials (SOP) approach was used to represent the PESs for the He_NH₂⁺ clusters,where V_HeiH2⁺(R_i,r_e) are the 3-body (3B) contributions represented by the CCSD(T)⁵⁶ potentials of the HeH₂⁺ complex,^21,42 with r_e being the equilibrium H₂⁺ bond length, and V_HeHe(R_ij) is the 2-body (2B) potential function for the He₂ given in ref. 44, the R_ij being the vectors between distinct (i,j) He atoms.

Fig. 1 shows the 3B and 2B potential curves (see the left panel) as a function of R or R_HeHe and the minimum energy path (see lower-right panel) of the 3B HeH₂⁺ CCSD(T)/CBS potential surface as a function of θ angle, with H₂⁺ fixed at r_e distance. The 3B PES shows two symmetric minima at linear configurations with R = 1.54714 Å at an energy of −332.06 meV, connecting by a barrier at T-shaped geometry and with the He atom at R = 1.57219 Å from the H₂⁺ molecular ion and energy of −42.68 meV (see inset plot). The 2B potential curve corresponding to the He–He interactions presents a very weak well at an energy of just 0.95391 meV for R_HeHe = 2.95 Å.

2.2 Ground state quantum DMC calculations

The diffusion Monte Carlo (DMC) method (see ref. 50 and references therein) is a stochastic approach that projects out the ground state of a quantum system by evolving the time-dependent Schrödinger equation in imaginary time, τ = it/ħ,where H is the Hamiltonian, R is the 3N-dimensional configuration vector, and E_T is a reference energy. In the long-time limit (τ → ∞), all excited-state components decay exponentially, and the wave function converges to the ground state. The propagation is carried out using a short-time approximation of the Green's function, which, in the importance-sampled form, reads

G(R,R′,δτ) ≈ G_dd(R,R″,δτ/2)G_b(R,R′,δτ)G_dd(R″,R′,δτ/2) (3)

whereincluding diffusion, drift, G_dd, and branching G_b terms, with D = ħ²/(2M) and M is the mass of the particles. Each walker undergoes random displacements drawn from a Gaussian distribution (diffusion), drifts along the gradient of the trial wave function (drift velocity v_D = 2∇ ln|Ψ_T|), and replicates or dies according to the local energy,which controls the branching weight. This importance-sampling procedure dramatically improves convergence and stability, particularly for systems with strongly anisotropic or weakly bound interactions.^51,52

For doped helium clusters, the Rigid Body DMC (RBDMC) formulation treats the molecular impurity as a rigid rotor coupled to He atoms.^53,54 This standard approximation is routinely employed to decouple the high-energy intramolecular stretching vibrations from the low-energy, large-amplitude intermolecular motions of the surrounding He atoms around the cation. Thus, the Hamiltonian separates translational, rotational, and interaction contributions as

where m and M are the He and cations' masses, I is the molecular moment of inertia of the cation, and V is the total interaction potential including He–He and H₂⁺–He terms^21,42,44 (see eqn (1)). The rotational kinetic energy is expressed in terms of the angular ϕ_x and ϕ_y coordinates (see upper-right panel in Fig. 1), which parameterize rotations around the body-fixed x and y axes. Importance sampling is implemented with a Jastrow-type guiding wave function to reduce variance and stabilize branching. Given the angular anisotropy in the H₂⁺–He interaction, we allow an explicit angular dependence in the 3B term, while the 2B HeHe term captures short and medium range correlations, with

The 2B correlation factors are parameterized as

and the anisotropic dopant–He contribution is expanded in even-order Legendre polynomials,

The radial functions u_l(R) adopt the flexible form

with variational parameters {b_l,ν_l,α_l,μ_l,β_l} optimized to balance compactness and accuracy.^52,55–58 With this trial wave function, the derivatives appearing in the drift velocity and local energy are computed analytically. The generalization of the short-time Green's function, eqn (3), to the Hamiltonian (7) is straightforward. The sampling of the rotational Green's function in the RBDMC is carried out by rotating the body fixed frame by using the corresponding rotation matrices for the sampled angles ϕ_x and ϕ_y at each step. To reduce the time step error, half of a rotation is done first in one axis, then the full rotation around the other axis and finally the other half of the rotation around the first axis. The coordinates of the molecule in the space frame are calculated straightforwardly by using the position of the center of mass of the molecule and the orientation of the body frame. This RBDMC implementation yields statistically exact ground-state energies and reliable structural properties for He clusters doped with molecular or ionic impurities, overcoming the limitations of conventional DMC for highly anisotropic interactions.⁵⁴ Regarding the computation of the He densities, the DMC mixed estimator was employed. While these quantities are affected by the trial wave function bias, we find that, for the variational wave functions employed here, no qualitative differences are observed between mixed and extrapolated density estimators. In particular, the predicted radial positions of density maxima and the number of He atoms per solvation shell are the same within statistical uncertainty.

2.3 Results and discussion

We first perform a variational Monte Carlo (VMC) calculation, where we optimize the parameters of the wave function by minimizing the total energy of the corresponding N size cluster. The order of the expansion in Legendre polynomials (see eqn (10)), l_M is set to 4, and thus 18 parameters were used for all clusters studied. It is worth remembering that the guiding function is employed to guide the RBDMC calculation, in such a way that an involved expansion would lead to an excessively long computational time. For the RBDMC calculation, an initial population of 2000 walkers is propagated by using 101 independent blocks (the initial one is not employed for properties calculation) with 10⁶ steps per block, and δτ = 10⁻⁴ meV⁻¹.

In Tables S1 and S2 (see the SI) we list the total, potential, kinetic, rotational and translational average energies as obtained from the RBDMC calculations for the H₂⁺He_N and D₂⁺He_N clusters up to N = 30.

In the upper panel of Fig. 2, we show the ground state quantum DMC, E_N, energies and their average per added He atom, , for each H₂⁺He_N and D₂⁺He_N clusters up to N = 30. For comparison reasons, the corresponding potential BH optimized energies, including the harmonic [scr E, script letter E]_ZPE correction for the H₂⁺He_N, are also displayed. These values, together with the classical optimized structures, are given in ref. 21. Both the H₂⁺He_N and D₂⁺He_N energy curves exhibit a similar size dependence. Furthermore, the DMC-fixed energies, clamped molecule, are found to closely follow the quasiclassical BH + [scr E, script letter E]_ZPE values, particularly for small clusters with N < 7. We marked four changes in the slope of the potential BH optimized energy curve (see colored solid lines in the upper panel of Fig. 2) at N = 2, 7, 13 and around 20, which can be associated with the formation of particularly stable structures or completion of solvation shells. These changes are also apparent in the quantum DMC curves for H₂⁺He_N and D₂⁺He_N. However, except for the first change, the remaining slope variations are considerably less pronounced than those in the BH-optimized profile, reflecting the smoothing of the energetic signatures of structural transitions, likely due to quantum delocalization.

Fig. 2 (Upper panel) DMC energies (E_N in meV) and average energy per He atom (E_N/N in meV) values for the H₂⁺He_N and D₂⁺He_N clusters as a function of cluster size. Comparisons with potential minima energy as obtained from the BH optimizations,²¹ including the harmonic [scr E, script letter E]_ZPE corrections, are also displayed. (lower panel) Single-atom evaporation energies for the H₂⁺He_N and D₂⁺He_N clusters as computed by the present DMC calculations in comparison with previous BH, BH + [scr E, script letter E]_ZPE and semiclassical FH2 results,21 as well as with the experimentally recorded ion–yield distributions from ref. 16.

The lower panel of Fig. 2 shows the computed single-atom evaporation energies, E_N−1 − E_N, as a function of cluster size, as obtained by the present DMC calculations, in comparison with the previously reported BH, BH + [scr E, script letter E]_ZPE, and semiclassical FH2 simulations at T = 2 K.²¹ The DMC curves present pronounced steps at small cluster sizes, N = 2 and 6, in agreement with classical and semiclassical results. For larger clusters, however, the DMC curves for both H₂⁺He_N and D₂⁺He_N are relatively flat, with a slight decreasing trend up to N = 30. In contrast, both BH and FH2 curves show oscillatory peaks for the N = 10–20, and weak local maxima up to N = 30. The BH + [scr E, script letter E]_ZPE energies closely follow the DMC ones, although the corresponding curve exhibits well pronounced drops at N = 6, 10, 13, and 19, as well as smaller ones at 23, 26 and 28.

In view of the discrepancies observed between the distinct computational studies, and in order to further examine them, we compared our DMC single-atom evaporation energy data with the experimental ion yield distributions reported in the literature¹⁶ (see lower panel in Fig. 2). For both H₂⁺He_N and D₂⁺He_N ion series, the reported¹⁶ distributions show an abrupt drop from N = 13 to 14, and a peak at N = 19, whereas the D₂⁺He_N distribution exhibits a generally flat behavior with only a weak maximum at N = 19. However, according to the authors, due to contamination in the ion yield distributions of both H₂⁺He_N and D₂⁺He_N, only the peak at N = 19 in the D₂⁺He_N series can be considered a reliable indicator of a magic cluster size. One should also note that the presence of the cationic dopant and the anisotropic He–cation interaction significantly modify the even–odd oscillations characteristic of pure He clusters, leading to the size-dependent patterns observed in our results, which are also comparable with the experimental ion–yield distributions.¹⁶

By analyzing the energy trends obtained from both computations and experiment, no definitive conclusions can be drawn regarding the energetic stability of the medium-sized clusters studied here. We therefore turned to an examination of the computed DMC H₂⁺ and D₂⁺He_N cluster structures, with the aim of gaining further insights into their structuring and possible shell-closure features. Fig. 3 presents the 3D spatial probability distributions obtained from the quantum DMC calculations for selected cluster sizes. One clearly observes highly compact and well-defined He structures surrounding the H₂⁺ cation. Depending on the cluster size, the formation and shape of the structures vary, with the first two He atoms remaining strongly bound to the H₂⁺ at linear configurations, and as the size of clusters increases, this [He–H–H–He]⁺ structure remains as a structural core of the larger clusters up to N = 30. One can observe that in the N = 7 cluster, the He atoms are distributed in a central and almost closed ring around this core, and as the number of He atoms increases, the resulting 3D spatial distributions show increasingly well-defined and compact regions.

Fig. 3 3D spatial probability distributions of He atoms around the H₂⁺ cation as obtained from the DMC simulations for the indicated H₂⁺He_N clusters. Both side and top views are shown. The He densities are computed by the DMC simulations, with the isosurface values being 0.025, 0.015, 0.009, 0.005, 0.005 and 0.005N_He per ų for N = 7, 9, 13, 17, 19 and 27, respectively.

In particular, when comparing these quantum distributions with the classical and semiclassical FH2 results reported in ref. 21, the 3D DMC distributions' notable differences emerge, especially with respect to the classical distributions, while similarities appear with those from the FH2 calculations. As the number of He atoms increases, the 3D classical distributions develop an increasing number of rings, reaching up to three inner rings at each side of the cation for N = 20, with an additional external ring placed around the center of the cation appearing for N = 27. In contrast, the semiclassical FH2 distributions are less structured, exhibiting incomplete and fewer He rings up to N = 19, while for N ≥ 20, the He atoms are preferentially aligned along ten radial directions around the cation. Overall, the growth pattern predicted by the FH2 calculations for the H₂⁺He_N clusters shows certain resemblance to the 3D DMC distributions, provided that important quantum effects are not explicitly included in the FH2 treatment. This comparison also highlights the degree of He delocalization in the DMC results. One can observe in the N = 13 and 17 clusters that the two rings visible in the FH2 distributions merge into a single and more diffuse one in the DMC plots. Moreover, the transition from ring to radials starts for N ≥ 19 in the DMC calculations. In these larger clusters, the inner layer consists of five well-defined and compact density maxima around the H₂⁺ core, whereas in the outer layer, the substantial exchange of He atoms between adjacent radial directions produces five broad, diffuse radial regions in the DMC results, compared to ten in the FH2 distributions.

The observed structural features, as well as the discrepancies between the different methods, highlight the delicate balance between the He–He and He–H₂⁺ forces, and the strong influence of nuclear quantum effects. It is also worth emphasizing that quantum delocalization of He dominates over any thermal effects present in the MD and FH2 data.

Furthermore, the persistent presence of the He–[He–H–H–He]⁺ structural core in all cluster size distributions provides an important insight into the SOP approach, suggesting that higher-order many-body contributions beyond the 3B He–H₂⁺ terms may also need to be considered in the PES expansion. Given the strong dependence of the cluster structure of the underlying PES topology, in the case that such contributions are significant for accurately describing medium-size clusters, then refinement of the present approach would be warranted.

Next, we examined the isotopic dependence of the He atom structuring around the D₂⁺ cation. Fig. 4 presents the radial R_HeHe probability distributions obtained from the quantum DMC calculations for both H₂⁺He_N (see upper panel) and D₂⁺He_N (see lower panel) clusters at the indicated N up to 30. Although only subtle differences can be discerned from these 1D distributions, the isotopic effects become more apparent when inspecting the 3D spatial probability distributions. To illustrate this, Fig. 5 displays representative 3D plots corresponding to the D₂⁺He₇, D₂⁺He₁₇, D₂⁺He₁₉ and D₂⁺He₂₀ clusters.

Fig. 4 Radial R_HeHe probability distributions obtained from the quantum DMC calculations for H₂⁺He_N (upper panel) and D₂⁺He_N (lower panel) for the indicated N cluster size.

Fig. 5 3D spatial probability distributions of He atoms around the D₂⁺ cation as obtained from the DMC simulations for the indicated D₂⁺He_N clusters. Both side and top views are shown. The He densities are computed by the DMC simulations, with the isosurface values being 0.025N_He per ų for N = 7, and 0.005N_He per ų for N = 17, 19 and 20.

By comparing these with their H₂⁺He_N counterparts, one can see structural differences, which weren't appreciated in the 1D profiles shown in Fig. 4. For instance, the D₂⁺He₇ presents a complete He ring centered at the D₂⁺ cation, whereas H₂⁺He₇ (see the top panel in Fig. 3) shows a discontinuous one. Conversely, the views of the D₂⁺He₁₉ differ from those of the H₂⁺He₁₉, but they are close to those corresponding to the H₂⁺He₁₇, presenting a central ring accompanied by diffuse ones at each side of the cation.

Fig. 6 (middle panels) shows the radial probability distributions obtained from the quantum DMC calculations as a function of R for selected size H₂⁺He_N and D₂⁺He_N clusters. As the cluster size increases, new peaks gradually emerge in these 1D distributions, with the R values extending up to about 7 Å for N = 30. To quantify the helium probability density and determine how many atoms populate each spatial layer, we performed partial integrations over the selected areas surrounding the H₂⁺ or D₂⁺ cation that correspond to these features. The vertical dashed lines in the middle panels mark the three layers defined by the R values at which new peaks appear in the H₂⁺He_N and D₂⁺He_N distributions. The resulting fractions of He atoms in each layer are shown in the upper panel of Fig. 6 where one can see a decrease in population in the first and second layers for N > 2 and N > 9, respectively, in the H₂⁺He_N clusters. The lower panel further reports the number of He atoms assigned to each of the three layers as a function of N for both H₂⁺He_N and D₂⁺He_N clusters, together with the corresponding H₂⁺He_N FH2 predictions. Both the semiclassical FH2 and quantum DMC results indicate that the first layer contains a fixed number of two He atoms for all cluster sizes, with no exchange with the outer layers. This reflects the stronger He–H₂⁺ interaction relative to the He–He interaction, which keeps one He atom on each side of the H₂⁺ or D₂⁺ molecular ion. The second and third layers emerge for clusters with N ≥ 3 and N ≥ 10, respectively. For N > 7, the differences between FH2 and DMC seen in the full 3D spatial distributions are also apparent in the layer populations, with FH2 predicting a larger number of He atoms in the second layer, whereas DMC places more atoms in the third one. In particular, both methods exhibit similar dependencies on the cluster size, with the number of He atoms in the second layer reaching a plateau of nearly ten He atoms for N > 17 from the FH2 data, and between five to six He atoms for N > 10 in the DMC calculations. In contrast, the population of the third layer grows almost linearly with N for N > 9 and 10 in both approaches, up to the largest N = 30 size considered here.

Fig. 6 Radial probability distributions (in R) obtained from the quantum DMC simulations for the indicated H₂⁺He_N and D₂⁺He_N, with N = 1–30 (middle panels). The upper panel shows the fraction of He atoms in each solvation layer as a function of N, computed by integrating the density in each layer. The vertical dashed lines in the middle panel mark the separation of three peaks in the R distributions. The lower panel shows the number of He atoms in solvation layers as a function of cluster size for the indicated He_NH₂⁺ and He_ND₂⁺, with N = 1–30, obtained from the DMC calculations in comparison with semiclassical FH2 data for the He_NH₂⁺ case.

On the basis of our analysis of average and single-atom evaporation energies, as well as structural distributions of the He probability densities, we found that both H₂⁺He_N and D₂⁺He_N clusters exhibit the coexistence of compact and diffuse spatial regions. These features may be interpreted as signatures of solidlike or liquidlike behavior of the He atoms surrounding the cations. To quantify this tendency more rigorously, we employed robust structural indicators that can capture such distinctions, providing additional insights. In particular, we calculated the Lindemann index,^59–61 defined as, , where r_ij refers to R and R_HeHe distances. The employed definition corresponds to the pair-distance fluctuation measure, commonly used in finite cluster studies. These δ indices quantify atomic fluctuations relative to their time-averaged positions, and their magnitude serves as an indication of structural rigidity or flexibility. A solidlike behavior is associated with highly localized He atoms, whose relative positions fluctuate only weakly, leading to small δ values (typically <0.1). Here, the term “localized” refers solely to the reduced spatial spread of the probability density for the inner He atoms, as inferred from their narrower radial distributions. In contrast, liquidlike behavior corresponds to significant fluctuations of the He atoms around their mean positions, yielding larger δ values (typically >0.2). We should note that the threshold values of δ are system-dependent, and for He-containing systems, one expects higher fluctuations than in classical ones.

Fig. 7 displays the system-averaged δ values for specific H₂⁺He_N and D₂⁺He_N clusters from the DMC simulations, in comparison with those from the FH2 results (see upper panel). The figure also shows the individual He atom δ values obtained from DMC for both isotopologues up to N = 30 (middle and lower panels). For all studied H₂⁺He_N clusters, the averaged δ values for R_HeHe (dashed black line) and R distance fluctuations (solid black line), lie between 0.11 and 0.165, slightly above the typical 0.1 threshold. This indicates an overall mixed structural character of the He atoms surrounding the H₂⁺ cation, with N = 12 cluster showing an enhanced He atom mobility. These values support the presence of compact He arrangements that nonetheless allow for the exchange of He atoms between them in the outer layers. For the D₂⁺He_N clusters, the corresponding system-averaged δ values are generally lower, with the smaller N = 6 and 7 clusters showing a more pronounced solidlike behavior.

Fig. 7 (upper panel) Lindemann index, δ, average values for both R_HeHe (dotted-triangle lines) and R (solid-circle lines) distance fluctuations for the indicated H₂⁺He_N and D₂⁺He_N cluster, as obtained from the DMC simulations. Comparison with the FH2 results is also presented for the He_NH₂⁺ case. (middle and lower panels) Lindemann index, δ, values for the R_i distance fluctuations for each He atom in the indicated He_NH₂⁺ and He_ND₂⁺ clusters, respectively.

The individual He atom δ values shown in the middle and lower panels of Fig. 7 reinforce this picture. In both isotopologues, most He atoms in each cluster show 0.1 ≤ δ ≤ 0.2, highlighting the coexistence of solidlike and liquidlike behavior. At least two He atoms per cluster have δ ≤ 0.1, consistent with the presence of well-localized, solidlike structures in the inner shell (see also Fig. 6). Only one or two He atoms in the N = 12, 13, and 20 H₂⁺He_N clusters show larger fluctuations with δ values of around 0.2, indicating a more liquidlike character for those specific atoms. In contrast, for the D₂⁺He_N clusters, only one He atom in the N = 13, 19, and 20 systems exhibits intermediate fluctuations with 0.15 < δ < 0.2.

3 Summary and conclusions

We carried out quantum DMC simulations to investigate the microsolvation, structure, and stability of the protonated hydrogen cations H₂⁺ and D₂⁺ embedded in helium clusters. Using a SOP expansion including 2B and 3B CCSD(T)/CBS-quality interaction terms, we quantified size-dependent energetics and analyzed the spatial probability densities that govern the organization of helium around these light quantum ions. The DMC results highlight the dominant role of nuclear quantum effects in shaping both the inner and outer solvation layers, and delineate the limits of the semiclassical FH2 effective potential approach in capturing the full quantum picture.

Our simulations show that He atoms do not distribute uniformly around the cations. Instead, the clusters form distinct solvation layers whose shape and compactness depend on the number of He atoms and the isotopic cation. For both ions, the first two He atoms are strongly localized and form a compact inner core, whereas additional He atoms increasingly display a mixture of solidlike and liquidlike behavior. This interplay leads to highly delocalized outer regions with significant He-exchange, preventing the emergence of clear magic-number structures in both species. The observed isotopic differences demonstrate how even small variations (in mass in this case) can alter the balance between cation–He anisotropy and He–He interactions.

The predicted solvation patterns and energy shifts provide a theoretical foundation for interpreting spectroscopic measurements of He-solvated ions. The computed evaporation-energy trends reproduce qualitative features of experimental ion yield distributions, including weak ion abundance anomalies in the He_NH₂⁺ and He_ND₂⁺ series. Such comparisons also show that quantum simulations are essential for validating data-driven potential models, and for assessing how potentials influence experimentally recorded spectral signatures. Although enhanced relative stability is found for certain medium-sized clusters, a comparison with available experimental data still remains inconclusive. The discrepancies between theoretical predictions and experimental data suggest that further refinement of the potential model approaches is needed to fully capture the underlying stabilization mechanisms and microsolvation dynamics of He_NH₂⁺ clusters. In particular, our results indicate that the [He–H–H–He]⁺ structural motif acts as the dominant core for all cluster sizes studied, which may point to limitations in the current SOP representation. At the same time, more detailed experimental information, especially high-resolution, contamination-free mass spectra of these He-doped nanoclusters, would be especially valuable for identifying and improving the remaining sources of error.

Beyond clarifying size-dependent stability, the present study highlights the importance of higher-order nonadditive contributions within many-body potential expansions, offers guidance for the development of next-generation interaction models, and supports future investigations of quantum solvation, isotopic effects, and ion–He interactions in cold environments.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data are generally included in the article itself and from the corresponding author upon reasonable request.

Supplementary information (SI): kinetic, rotational, translational, potential and total energy values together with their corresponding errors from the DMC simulations. See DOI: https://doi.org/10.1039/d5cp04786b.

Acknowledgements

The authors thank the Centro de Calculo del IFF/SGAI-CSIC and CESGA-Supercomputing centre for allocation of computer time. We acknowledge financial support from the Comunidad de Madrid grants ref: IND2018/TIC-9467 and INVESTIGO ref: INV23-IFF-M2-09, MICIU grant no. PID2024-155666-NB-I00, CSIC-PEICT ref: 2024AEP119 and COST Actions CA21101(COSY) and CA21126(NanoSpace).

Notes and references

1. J. P. Toennies and A. F. Vilesov, Angew. Chem., Int. Ed., 2004, 43, 2622–2648.
2. Y. Xu, W. Jäger, J. Tang and A. R. W. McKellar, Phys. Rev. Lett., 2003, 91, 163401.
3. M. Rossi, M. Verona, D. E. Galli and L. Reatto, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 212510.
4. M. Fárník and J. Toennies, J. Chem. Phys., 2005, 122, 014307.
5. S. Moroni and S. Baroni, Comput. Phys. Commun., 2005, 169, 404–407.
6. C. Di Paola, F. Sebastianelli, E. Bodo, F. A. Gianturco and M. Yurtsever, J. Chem. Theory Comput., 2005, 1, 1045–1054.
7. M. Y. Choi, G. E. Douberly, T. M. Falconer, W. K. Lewis, C. M. Lindsay, J. M. Merritt, P. L. Stiles and R. E. Miller, Int. Rev. Phys. Chem., 2006, 25, 15–75.
8. E. Bodo, E. Yurtsever, M. Yurtsever and F. A. Gianturco, J. Chem. Phys., 2006, 124, 074320.
9. F. Stienkemeier and K. K. Lehmann, J. Phys. B: At., Mol. Opt. Phys., 2006, 39, R127–R166.
10. S. Paolini, F. Ancilotto and F. Toigo, J. Chem. Phys., 2007, 126, 124317.
11. D. E. Galli, D. M. Ceperley and L. Reatto, J. Phys. Chem. A, 2011, 115, 7300–7309.
12. D. Verma, R. M. P. Tanyag, S. M. O. O'Connell and A. F. Vilesov, Adv. Phys.: X, 2019, 4, 1553569.
13. M. Gatchell, P. Martini, F. Laimer, M. Goulart, F. Calvo and P. Scheier, Faraday Discuss., 2019, 217, 276–289.
14. O. Asvany, S. Schlemmer, T. Szidarovszy and A. Czáscár, J. Phys. Chem. Lett., 2019, 10, 5325–5330.
15. T. M. Kojima, N. Kobayashi and Y. Kaneko, Z. Phys. D:At., Mol. Clusters, 1992, 23, 181–185.
16. P. Bartl, C. Leidlmair, S. Denifl, P. Scheier and O. Echt, ChemPhysChem, 2013, 14, 227–232.
17. L. Lundberg, P. Bartl, C. Leidlmair, P. Scheier and M. Gatchell, Molecules, 2020, 25, 1066.
18. R. Yanes-Rodríguez, R. Rodríguez-Segundo, P. Villarreal and R. Prosmiti, Eur. Phys. J. D, 2023, 77, 116.
19. F. Brieue, C. Schran and D. Marx, Phys. Rev. Res., 2023, 5, 043083.
20. R. Yanes-Rodríguez, P. Villarreal and R. Prosmiti, Phys. Chem. Chem. Phys., 2025, 27, 8259–8266.
21. M. J. M. de Oca-Estévez, P. Villarreal, J. Hernández-Rojas, J. Bretón, A. Sarsa and R. Prosmiti, Chem. Phys. Lett., 2025, 877, 142242.
22. J. P. Toennies, A. F. Vilesov and K. B. Whaley, Phys. Today, 2001, 54, 31–37.
23. E. Coccia, E. Bodo, F. Marinetti, F. A. Gianturco, E. Yildrim, M. Yurtsever and E. Yurtsever, J. Chem. Phys., 2007, 126, 124319.
24. P. Moroshkin, V. Lebedev and A. Weis, Phys. Rev. Lett., 2009, 102, 115301.
25. S. Müller, M. Mudrich and F. Stienkemeier, J. Chem. Phys., 2009, 131, 044319.
26. R. Prosmiti, G. Delgado-Barrio, P. Villarreal, E. Yurtsever, E. Coccia and F. A. Gianturco, J. Phys. Chem. A, 2009, 113, 14718–14729.
27. C. Callegari and W. E. Ernst, Handbook of High-Resolution Spectroscopy, Hoboken, NJ, 2011, pp. 1551–1594.
28. P. Bartl, C. Leidlmair, S. Denifl, P. Scheier and O. Echt, J. Phys. Chem. A, 2014, 118, 8050.
29. A. Czáscár, T. Szidarovszy, O. Asvany and S. Schlemmer, Mol. Phys., 2019, 117, 1559–1583.
30. I. Simkó, C. Schran, F. Brieuc, C. Fábri, O. Asvany, S. Schlemmer, D. Marx and A. G. Császár, Angew. Chem., Int. Ed., 2023, 62, e202306744.
31. A. Krapp, G. Frenking and E. Uggerud, Chem. – Eur. J., 2008, 14, 4028–4038.
32. M. Meuwly and J. M. Hutson, J. Chem. Phys., 1999, 110, 3418–3427.
33. W. P. Kraemer, V. Spirko and O. Bludský, Chem. Phys., 2002, 276, 225–242.
34. D. Koner, J. Veliz, A. van der Avoird and M. Meuwly, Phys. Chem. Chem. Phys., 2019, 21, 24976–24983.
35. K. Naskar, S. Ravi, S. Adhikari, M. Baer and N. Sathyamurthy, J. Phys. Chem. A, 2023, 127, 3832–3847.
36. N. Alharzali, H. Berriche, P. Villarreal and R. Prosmiti, J. Phys. Chem. A, 2019, 123, 7814–7821.
37. N. Alharzali, R. Rodríguez-Segundo and R. Prosmiti, Phys. Chem. Chem. Phys., 2021, 23, 7849–7859.
38. M. J. Montes de Oca-Estévez and R. Prosmiti, Front. Chem., 2021, 9, 187.
39. J. Bretón, J. Hernández-Rojas, M. I. Hernández, J. Campos-Martínez and T. González-Lezana, ChemPhysChem, 2023, 24, e202300425.
40. M. J. Montes de Oca-Estévez, A. Valdés and R. Prosmiti, Phys. Chem. Chem. Phys., 2024, 26, 7060–7071.
41. M. J. Montes de Oca-Estévez and R. Prosmiti, Molecules, 2024, 29, 4084.
42. M. J. Montes de Oca-Estévez and R. Prosmiti, Artif. Intell. Chem., 2024, 2, 100059.
43. M. J. Montes de Oca-Estévez, A. Valdés and R. Prosmiti, Molecules, 2025, 30, 2440.
44. R. A. Aziz and M. J. Slaman, J. Chem. Phys., 1991, 94, 8047.
45. A. Valdés, R. Prosmiti, P. Villarreal and G. Delgado-Barrio, J. Chem. Phys., 2005, 122, 044305.
46. A. Valdés, R. Prosmiti, P. Villarreal and G. Delgado-Barrio, J. Chem. Phys., 2006, 125, 014313.
47. Á. Valdés, P. Barragán, R. Pérez de Tudela, J. S. Medina, L. Delgado-Tellez and R. Prosmiti, Chem. Phys., 2012, 399, 39–45.
48. A. Valdés and R. Prosmiti, J. Phys. Chem. A, 2013, 117, 7217–7223.
49. R. Pérez de Tudela, P. Barragán, A. Valdés and R. Prosmiti, J. Phys. Chem. A, 2014, 118, 6492–6500.
50. B. Hammond, W. L. Jr and P. Reynolds, Monte Carlo Methods in ab initio Quantum Chemistry, Word Scientific, Singapore, 1994.
51. P. Niyaz, Z. Bačić, J. W. Moskowitz and K. E. Schmidt, Chem. Phys. Lett., 1996, 23.
52. A. Viel, M. Patel, P. Niyaz and K. Whaley, Comput. Phys. Commun., 2002, 24.
53. V. Buch, J. Chem. Phys., 1992, 726.
54. J. Gregory and D. Clary, Chem. Phys. Lett., 1994, 547.
55. D. Blume, M. Lewerenz, F. Huisken and M. Kaloudis, J. Chem. Phys., 1996, 8666.
56. A. Sarsa, K. E. Schmidt and J. W. Moskowitz, J. Chem. Phys., 2000, 44.
57. F. Paesani, F. A. Gianturco and K. B. Whaley, J. Chem. Phys., 2001, 10225.
58. A. Viel and K. B. Whaley, J. Chem. Phys., 2001, 10186.
59. G. Franke, E. Hilf and L. Polley, Z. Phys. D:At., Mol. Clusters, 1988, 9, 343–349.
60. T. Beck, J. Doll and D. Freeman, J. Chem. Phys., 1989, 90, 5651–5656.
61. P. Barragán, R. Pérez de Tudela, C. Qu, R. Prosmiti and J. Bowman, J. Chem. Phys., 2013, 139, 024308.

---