# DynHeMat: A program for zero-point averaged dynamics of pure and doped helium nanodroplets

David A Bonhommeau

► To cite this version:

David A Bonhommeau. DynHeMat: A program for zero-point averaged dynamics of pure and doped helium nanodroplets. Computer Physics Communications, 2026, 321, pp.110014. 10.1016/j.cpc.2025.110014. hal-05472935

HAL Id: hal-05472935

https://hal.science/hal-05472935v1

Submitted on 23 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Computer Programs in Physics

# DynHeMat: A program for zero-point averaged dynamics of pure and doped helium nanodroplets

David A. Bonhommeau

<figure><img src="image_2.png" /></figure>

Université Paris-Saclay, Univ Evry, CY Cergy Paris Université, CNRS, LAMBE, 91025, Evry-Courcouronnes, France

<figure><img src="image_3.png" /></figure>

Check for
updates

## ARTICLE INFO

Editor: Dr. J Ballantyne

**Keywords:**

- Helium nanodroplets
- Rare gas
- Alkali
- Molecular dynamics
- Zero-point energy
- Collision
- Coagulation
- Submersion

## ABSTRACT

DynHeMat is a parallel program aimed at modeling the dynamics of pure and doped helium nanodroplets (HNDs) by means of zero-point averaged dynamics (ZPAD), a method where the quantum nature of helium atoms is taken into account through the use of a He-He pseudopotential which includes zero-point effects of helium clusters on an average manner. Three He-He pseudopotentials, defined for applications in different contexts, are implemented. Large HNDs can be formed by successive coalescences of smaller HNDs keeping in mind that, depending on the HND size and He-He pseudopotential in use, the liquid character of the HND is more or less pronounced. Files containing the positions and velocities of HNDs formed with the three aforementioned He-He pseudopotentials are collected in a local databank, called ZPAD_DB. ZPAD simulations can be carried out at constant energy or temperature, then enabling the user to investigate collision, coagulation or submersion processes in pure or doped HNDs. Impurities can be rare-gas atoms (Ne, Ar, Kr, Xe and Rn), alkali atoms (Li, Na, K, Rb, Cs), or homogeneous clusters composed of such atoms. The program provides information on trajectories, namely positions, velocities, energies, radial distribution functions, and the initial distribution of HND surface atoms. Extension to other impurities or He-He pseudopotentials is made possible by the current structure of the program and keyword system.

## PROGRAM SUMMARY

**Program title:** DynHeMat

**CPC Library link to program files:** https://doi.org/10.17632/3hrfykstvr.1

**Licensing provisions:** GNU General Public License 3 (GPL)

**Programming language:** Fortran 90

**Nature of problem:**

Helium nanodroplets (HNDs) are large quantum systems containing from a few thousands to billion atoms. The more quantum approaches, like quantum Monte Carlo, time-dependent density functional theory and path integral molecular dynamics, are often limited to the treatment of a few hundreds or thousands atoms or to small statistics in terms of projectile velocities or impact parameters, for instance. On the contrary, a classical approach would enable simulations on larger systems provided that the quantum nature of helium atoms is included on an average manner in the calculation in order to ensure that the expected heliophilic or heliophobic nature of impurities can be maintained.

**Solution method:**

The zero-point averaged dynamics (ZPAD) includes zero-point effects on an average manner in classical simulations through the use of an effective He-He potential, and possibly effective He-impurity potentials, which makes the HND liquid and drastically improves the agreement with quantum calculations compared to standard classical simulations. Initially used to tackle the fragmentation of rare-gas clusters embedded in HNDs ionized by electron impact, the ZPAD approach was later applied to study the coagulation of rare-gas atoms in liquid helium and the collision of alkali atoms with small HNDs. In the current version of the program, impurities can be rare gases (from Ne to Rn) or alkali atoms (from Li to Cs). The program is parallelized with MPI in such a way that several independent trajectories can be performed simultaneously, each of them being distributed over a set of CPU cores, which allows simulations on large HNDs (N ≈ 10⁵) to be carried out in a reasonable amount of time.
### Additional comments including Restrictions and Unusual features:

The program was especially tested on x86 and ARM architectures under operating systems linux RedHat 9 and MacOS Ventura, and the number of trajectories run in parallel should not exceed 999. Unlike TDDFT, ZPAD does not take into account superfluidity and quantum vortices cannot be investigated with this approach. Moreover, a significant amount of time may be needed to produce pure HNDs before running ZPAD simulations to investigate collision, coagulation or submersion processes. To compensate for this drawback, a databank which contains XYZ files, including positions and velocities, of He<sub>N</sub> (1000 ≤ N < 90000) is made available to DynHeMat users.

### References:

1. Modeling the fragmentation dynamics of ionic clusters inside helium nanodroplets: The case of He<sub>100</sub>Ne<sub>1</sub><sup>+</sup>, D. Bonhommeau, P. T. Lake, Jr, C. Le Quiniou, M. Lewerenz and N. Halberstadt, J. Chem. Phys. **126**, 051104 (2007). DOI: 10.1063/1.2515225

2. Fragmentation of ionized doped helium nanodroplets: Theoretical evidence for a dopant ejection mechanism, D. Bonhommeau, M. Lewerenz and N. Halberstadt, J. Chem. Phys. **128**, 054302 (2008). DOI: 10.1063/1.2823101

3. Fragmentation dynamics of Ar<sub>4</sub>He<sub>1000</sub> upon electron impact ionization: competition between ion ejection and trapping, N. Halberstadt and D. A. Bonhommeau, J. Chem. Phys. **152**, 234305 (2020). DOI: 10.1063/5.0009363

4. Zero-point averaged dynamics of Ar<sub>n</sub><sup>+</sup>He<sub>1000</sub>: Dopant-size influence on potential-energy surfaces, mass spectra and fragmentation patterns, D. A. Bonhommeau, Chem. Phys. 550 (2021) 111307. DOI: 10.1016/j.chemphys.2021.111307

5. Clustering, collision, and relaxation dynamics in pure and doped helium nanoclusters: Density- vs particle-based approaches, E. García-Alfonso, M. Barranco, D. A. Bonhommeau, N. Halberstadt, M. Pi, et F. Calvo, J. Chem. Phys. **157**, 014106 (2022). DOI: 10.1063/5.0091942

6. Collision of rare-gas atoms on helium nanodroplets: Theoretical evidence for an efficient coagulation of heavy rare-gas atoms, B. Fixot, E. Louaas and D. A. Bonhommeau, J. Chem. Phys. **161**, 044307 (2024). DOI: 10.1063/5.0220027

7. Collision of cesium atoms on helium nanodroplets: Unraveling mechanisms for surface capture at experimental velocities, D. A. Bonhommeau, J. Chem. Phys. **161**, 184302 (2024). DOI: 10.1063/5.0231641

## 1. Introduction

Because of their superfluid nature and temperature of formation (T = 0.37 K) [1,2], helium nanodroplets (HNDs) are often considered as ideal low-temperature nanolaboratories to conduct spectroscopic studies on small molecules [3–5], investigate the active cooling of fragmentation processes [6–8], or produce exotic polymers [9,10], metallic or bimetallic nanoclusters or nanowires by harnessing vortices formed in liquid helium [11,12]. More recently, HNDs have also been used to study ultrafast rotational and vibrational dynamics of molecules (eg, methyliodide [13] or dimers like D<sub>2</sub> [14], D<sub>2</sub><sup>+</sup> [15], Na<sub>2</sub> [16] and In<sub>2</sub> [17]) in the presence of a dissipative environment.

These applications stimulated the development of a number of theoretical approaches although the largest HNDs (N ≈ 10<sup>11</sup>) [18] cannot be modeled yet. Quantum Monte Carlo simulations (QMC) on helium clusters containing at most a few hundreds atoms were undertaken in the 1990's to reproduce the excitation spectra of liquid <sup>4</sup>He [19,20] and qualitatively confirm the experimental number of <sup>4</sup>He atoms required for the appearance of superfluidity in liquid helium, namely N = 60 [2]. Variational and diffusion QMC calculations were later performed to tackle the properties of impurities inside HNDs. In particular, the size-dependent vibrational frequency shift of HF in helium clusters with up to 198 atoms was found in good agreement with experimental data obtained for larger clusters (N > 1000) [21] and small Ne<sup>+</sup>He<sub>n</sub> ions were surprisingly found to have no shell structure unlike Ar<sup>+</sup>He<sub>n</sub> ions [22].

The dynamics of doped helium clusters was sometimes described by hybrid models where helium atoms were modeled by path-integral centroid molecular dynamics (PIMD) while the impurity was described quantum mechanically. For instance, the dynamics of electronically excited Cl<sub>2</sub> embedded in He<sub>N</sub> (N ≤ 200) revealed that the cage effect exerted by helium atoms was drastically reduced, compared with the case where the system was treated classically [23], and the quantum treatment of the photoexcitation of a K atom at the surface of He<sub>300</sub> enabled to pinpoint the possible formation of K<sup>+</sup>He<sub>n</sub> exciplexes (n < 6) [24]. Moreover, the investigation of the submersion of Na<sub>n</sub> clusters in He<sub>N</sub> (N ≤ 1000) by PIMD confirmed that the critical size for

Na<sub>n</sub> submersion should be n<sub>c</sub> ≈ 20 [25], in agreement with experimental data [26] and theoretical predictions based on energy considerations [27]. More recently, the collision and coagulation of argon atoms in Ar<sub>n</sub>He<sub>1000</sub> (n ≤ 4) has been studied by PIMD and compared with results from other theoretical methods, namely time-dependent density functional theory (TDDFT) and zero-point averaged dynamics (ZPAD), a mixed quantum-classical method where the quantum nature of He atoms is taken into account through the use of a He-He pseudopotential [28]. At experimental projectile velocities, the formation of argon clusters within the HND was more favored in TDDFT calculations than in other methods and coagulation was always prevented in ZPAD simulations.

Unlike particle-based methods, TDDFT includes helium superfluidity in its original implementation which makes it a valuable tool to investigate hydrodynamic effects in HNDs, like quantum vortices. As an example, ring vortices can appear upon piercing of the HND by a Cs atom, which occurs at Cs speeds above ~ 200 m s<sup>-1</sup> [29]. Vortex arrays can also be defined in static He-DFT calculations to assess the ability of vortices to favor the coagulation of heliophilic species. When six Ar atoms are projected on a He<sub>5000</sub> nanodroplet containing a six vortex array, a loosely bound argon cluster is formed at the droplet center more slowly than in a vortex-free HND, and the argon cluster structure depends on initial conditions like the time delay between argon atom captures [30]. Finally, TDDFT has been recently employed to investigate the solvation of alkali ions inside helium in close collaboration with experimentalists [31–34].

The principle of ZPAD is simpler than that of PIMD or TDDFT since helium atoms and impurities evolve classically on an effective potential-energy surface which includes zero-point effects due to helium atoms on an average manner. The idea behind ZPAD was originally proposed by Sterling et al. who investigated the dynamics of pure and Li-doped solid H<sub>2</sub> where Li and H<sub>2</sub> molecules were modeled by gaussians whose center evolved classically [35]. In practice, the potential-energy surface (PES) of the system was defined as a sum of pairwise Li-H<sub>2</sub> and H<sub>2</sub>-H<sub>2</sub> pseudopotentials obtained by convolving the corresponding pair potential with a suitable gaussian function and the particles evolved classically
on this effective PES. This method was later generalized to model the photodissociation of HBr and HCl in $Ne_n$ ($n \le 300$) by converging Ne-Ne pseudopotentials iteratively [36], and the fragmentation dynamics of neon and argon clusters embedded in small HNDs ($N \le 1000$) upon electron impact ionization [37–40], then giving birth to the ZPAD approach to tackle the dynamics of pure and doped HNDs. However, the structures of $He_N$ ($N \ge 1000$) based on He-He pseudopotentials converged in these calculations were too solid and alternative pseudopotentials were proposed to correct this shortcoming. Panzou and Lewerenz adjusted a He-He pseudopotential to a Morse potential in order to make $He_N$ ($2500 \le N \le 9000$) liquid and investigate the photodissociation of molecular iodine in HNDs [41]. The same analytical form was later considered by Fixot et al. to define a He-He pseudopotential able to reproduce roughly the total energy of $He_{1000}$ obtained in He-DFT calculations [42]. ZPAD simulations performed with this pseudopotential yielded results in qualitative agreement with PIMD and TDDFT calculations regarding the coagulation of argon clusters in $He_{1000}$. Cs atoms colliding with $He_{1000}$ at experimental velocities were also found to be possibly captured by the HND in ZPAD simulations based on the same He-He pseudopotential, while TDDFT calculations suggested that the HND should be transparent for Cs atoms at such speeds [43].

The DynHeMat program synthesizes recent numerical developments on ZPAD for the study of neutral species (ie, rare-gas and alkali atoms and clusters) located within or at the surface of HNDs [28,42,43] and it is written in such a way that extensions to other impurities or He-He pseudopotentials should require little effort (see Appendix A for more information). The main features of ZPAD are recalled in Section 2. The general expression of the PES used in ZPAD simulations is first given in Section 2.1 together with references to all potential-energy curves implemented in DynHeMat. Information on initial conditions, equations of motion, and MPI parallelization are then reported in Sections 2.2 to 2.4, respectively. Instructions for compiling and running DynHeMat and test cases, information on the framework of input and output files, and details on the content of the ZPAD DB databank are provided in Section 3. Recent applications of DynHeMat to model the coagulation of rare-gas atoms in $RG_N$ $He_{1000}$ ($RG = Ne, Ar, Kr, Xe, or Rn$) [42] and the collision of Cs atoms with $He_{1000}$ [43] are summarized in Section 4. Coalescence processes of HNDs enabling the formation of large droplets commonly found in experiments, like those devoted to the study of alkali cluster submersion in liquid helium [26,44,45], are finally mentioned in the same section.

## 2. Method

### 2.1. Potential-energy surface

The potential-energy surface (PES) used in ZPAD simulations is a sum of He-He and He-impurity interaction potentials combined with a PES for impurities

$$U(\{\mathbf{R}\}) = V_{\text{impurities}}(\{\mathbf{R}_i\}) + \sum_{a=1}^{N} \sum_{i=1}^{n} V_{\text{He-impurity}}(R_{ai}) + \sum_{a=1}^{N} \sum_{\substack{b=1 \\ b>a}}^{N} V_{\text{He-He}}^{\text{eff}}(R_{ab}) \quad (1)$$

where $N$ and $n$ are the numbers of He atoms and impurities, respectively, $\{\mathbf{R}_i\} = \{\vec{R}_1, \dots, \vec{R}_n\}$ represent impurity atomic positions, and $\{\mathbf{R}\}$ represents both helium and impurity atomic positions. Four choices are allowed for He-He potential-energy curves: (i) an improved Tang-Toennies potential [46] which can be used to carry out reference classical simulations with no inclusion of zero-point effects (labelled “TT2”), (ii) an effective potential (or pseudopotential) converged iteratively by truncating the He wavefunction at a distance of $1.6$ Å at each iteration [40] (labelled “IPTW”), (iii) an effective potential adjusted on a Morse potential by Panzou and Lewerenz [41] (labelled “PL”), (iv) a modified version of the PL pseudopotential proposed by Fixot et al. [42] to better reproduce the total energy of $He_{1000}$ found in He-DFT calculations (labelled “mPL”). He-impurity potential-energy curves depend on the nature of the impurity. He-RG interactions ($RG = Ne, Ar, Kr, Xe, and Rn$) can be either described by a Tang-Toennies potential [46] or by a pseudopotential fitted to a Morse potential after convolving the Tang-Toennies potential with the square of an approximate He wave function and a Dirac function for the rare-gas atom according to [42]

<table><caption>**Table 1**

Chemical symbols of elements which are currently available in DynHeMat.</caption><thead><tr><th>Element</th><th>DynHeMat symbols</th></tr></thead><tbody><tr><td>Helium</td><td>He, 4He</td></tr><tr><td>Neon</td><td>Ne, 20Ne, 21Ne, 22Ne</td></tr><tr><td>Argon</td><td>Ar, 36Ar, 38Ar, 40Ar</td></tr><tr><td>Krypton</td><td>Kr, 78Kr, 80Kr, 82Kr, 83Kr, 84Kr, 86Kr</td></tr><tr><td>Xenon</td><td>Xe, 124Xe, 126Xe, 128Xe, 129Xe, 130Xe, 131Xe, 132Xe, 134Xe, 136Xe</td></tr><tr><td>Radon</td><td>Rn</td></tr><tr><td>Lithium</td><td>Li, 6Li, 7Li</td></tr><tr><td>Sodium</td><td>Na, 23Na</td></tr><tr><td>Potassium</td><td>K, 39K, 40K, 41K</td></tr><tr><td>Rubidium</td><td>Rb, 85Rb, 87Rb</td></tr><tr><td>Cesium</td><td>Cs, 133Cs</td></tr></tbody></table>

$$V_{He-RG}^{eff}(R) = \int V_{He-RG}(\|\vec{R} - \vec{r}\|) \psi^2(\vec{r}) d^3r. \quad (2)$$

At the beginning of ZPAD simulations, the DynHeMat program checks the compatibility between the selected He-He and He-RG pseudopotentials because the values of the He wave function used in Eq. (2) depends on the He-He potential-energy curve. For alkali metals (Alk = Li, Na, K, Rb, and Cs), the He-Alk analytical expressions proposed by Patil are considered [47]. No effective potential was created for alkali atoms because these atoms are heliophobic and Eq. (2) originally assumes that the impurity is surrounded by helium atoms.

The PES for impurities is non-zero when $n \ge 2$. For rare gases, it is computed as a sum of Tang-Toennies RG-RG potentials [46]

$$V_{RG_n}(\{\mathbf{R}_i\}) = \sum_{i=1}^{n} \sum_{\substack{j=1 \\ j>i}}^{n} V_{RG-RG}(R_{ij}). \quad (3)$$

On the contrary, alkali clusters are described by second moment approximation (SMA) of the density of states to the tight-binding model (TB) [48]

$$V_{Alk_n}(\{\mathbf{R}_i\}) = \epsilon_0 \sum_i \sum_{\substack{j=1 \\ j \neq i}}^{n} \exp \left[ p \left( 1 - \frac{R_{ij}}{r_0} \right) \right] - \zeta_0 \sum_i \left\{ \sum_{\substack{j=1 \\ j \neq i}}^{n} \exp \left[ 2q \left( 1 - \frac{R_{ij}}{r_0} \right) \right] \right\}^{1/2} \quad (4)$$

where $\epsilon_0$, $\zeta_0$, $r_0$, $p$ and $q$ are parameters whose values depend on the alkali atom considered. A distance-dependent tight binding model (DDTB) devised by Calvo and Spiegelmann for neutral and ionic sodium clusters [49] is also implemented in DynHeMat. When keeping only energy terms specific to $Na_n$, this model corresponds to a reparameterization of the SMATB model where parameters $\epsilon_0$ and $\zeta_0$ are replaced by $\epsilon_0^+ = \gamma\epsilon_0$ and $\zeta_0^+(n) = \gamma\zeta_0(1 + \alpha/n^\beta)$ with $\alpha = 6.6$, $\beta = 1.27$, and $\gamma = 1.0055$.

### 2.2. Initial conditions

#### 2.2.1. Initial positions and isotopic masses

Atomic symbols and initial atomic positions of the system under interest must be provided in a configuration file before starting a simulation (see Section 3.2 for details on the name and content of this input file). Atomic symbols define the masses of the isotopes to be considered in the calculation, standard chemical symbols (eg, “Ar” for argon) being reserved for the selection of atomic masses based on natural isotopic abundances. An exception is made for radon for which the longest-lived isotope, namely $^{222}$Rn ($t_{1/2} \approx 3.8$ days), is only included in DynHeMat. Table 1 collects all atomic symbols which are recognized in DynHeMat.
### 2.2.2. Atomic velocities

Atomic velocities can be supplied by the user in the input configuration file after atomic positions but this is not mandatory because these velocities can be randomly generated. The Maxwell-Boltzmann distribution at temperature T for the ith velocity component of a system containing N<sub>tot</sub> particles (atoms or rigid molecules) with mass m is given by

$$p(v_i) = \sqrt{\frac{m}{2\pi k_B T}} \exp\left(-\frac{mv_i^2}{2k_B T}\right) \quad (5)$$

which leads to $\langle v_i^2 \rangle = \frac{k_B T}{m}$. In DynHeMat, random velocity components corresponding to a given temperature T can be defined according to the formula

$$v_i = \xi_i \sqrt{\frac{k_B T}{m_i}} \quad (6)$$

where $\xi_i$ is a real random number describing a normal distribution obtained by applying the Box-Muller method to numbers from a uniform (0,1) distribution, as explained in the DLARNV LAPACK subroutine [50], and $m_i$ can be the mass of a He atom or that of an impurity. In practice, the random number generator is initialized by four integral seeds $S_i$ with $i \in [1, 4]$, among which $S_4$ should be odd. Since $N_{traj}$ trajectories can be run in parallel with DynHeMat, $N_{traj}$ sets of integral seeds $S'_{ij}$ ($i \in [1, 4]$ and $j \in [1, N_{traj}]$) can be generated by fixing four scaling factors $\sigma_i$ according to

$$S'_{ij} = S_i + j \times \sigma_i \quad \text{if } i \in [1, 3] \quad (7a)$$

$$S'_{4j} = S_4 + 2j \times \sigma_4 \quad (7b)$$

because $S'_{4j}$ should be odd just like $S_4$.

### 2.2.3. Initial total angular momentum

After producing random velocities, the angular momentum of the system may be large and it may be practical to fix it to a given value. As an example, ZPAD simulations for the collisions of rare-gas and Cs atoms with pure or doped He<sub>1000</sub> discussed in Section 4 were carried out after fixing the total angular momentum $\vec{L}$ of the target HND to zero. The approach followed to fix $\vec{L}$ is based on the work of Jellinek and Li to separate the overall rotation in a N-body system [51]. It was successfully applied to isolate the vibrational component of nonadibatic coupling vectors used to adjust velocities at hopping events during surface-hopping simulations devoted to the fragmentation dynamics of excited Ne<sub>n</sub><sup>+</sup> in the gas phase [52]. In practice, the new velocity vector $\vec{v}_i$ ($i \in [1, N_{tot}]$) is given by

$$\vec{v}_i = \vec{v}_{i0} + (\vec{\omega} - \vec{\omega}_0) \times \vec{r}_i \quad (8)$$

where $\vec{r}_i$ and $\vec{v}_{i0}$ are the position and initial velocity vectors of atom i while $\vec{\omega}_0$ and $\vec{\omega}$ are the initial and target angular velocities defined by their relation with the inertia matrix $\vec{I}: \vec{I} \cdot \vec{\omega}_0 = \vec{L}_0$ and $\vec{I} \cdot \vec{\omega} = \vec{L}$.

The generation of initial velocities and the adjustment of the total angular momentum of the system are optional calculations in DynHeMat. On the contrary, positions and velocities of the system are systematically placed in the center-of-mass frame at the end of the preparatory phase of the simulation. This was not found to drastically alter the value of the total angular momentum in ZPAD simulations.

### 2.2.4. Initial HND radius

After preparing an initial HND configuration, it is possible to assess its radius by averaging the distance between surface particles and the HND center of mass, where a particle is assumed to lie on the HND surface if it is accessible by a fictitious probe particle approaching the HND from above its surface [53]. This method requires to define the diameter of a typical particle (ie, a He atom here) and that of a smaller probe particle. The values of these diameters depend on the selected He-He pseudopotential and are hardcoded in DynHeMat. Triplets of particles in contact with the probe particle are then pinpointed and assumed to

lie on the HND surface if the probe particle does not overlap with any other He atom belonging to the HND. Since this algorithm is based on the definition of a single type of cluster and probe particles, it is better suited to estimate the radius of pure HNDs or HNDs doped with helophilic species, whose surface is only composed of He atoms. Note that an average HND radius can also be computed at the end of the equilibration phase (see Section 3.2.2).

## 2.3. Molecular dynamics and thermodynamic ensembles

### 2.3.1. Integration of classical equations of motion

In ZPAD simulations, helium and impurity nuclei evolve classically. The time evolution of a system composed of N He atoms and n impurities is determined by the classical equations of motion expressed in atomic units:

$$\frac{\partial r_i}{\partial t} = v_i \quad (9a)$$

$$\frac{\partial p_i}{\partial t} = -\nabla_i U(\{\mathbf{R}\}) \quad (9b)$$

where $i \in [1, 3N_{tot}]$ with $N_{tot} = N + n$ and $\nabla_i U(\{\mathbf{R}\})$ is the analytical derivative of the effective potential-energy surface $U(\{\mathbf{R}\})$ with respect to coordinate $r_i$. The classical dynamics is performed by using a velocity Verlet integrator

$$r_i(t + dt) = r_i(t) + v_i(t)dt - \frac{dt^2}{2m_i}\nabla_i U(\{\mathbf{R}(t)\}) \quad (10a)$$

$$v_i(t + dt/2) = v_i(t) - \frac{dt}{2m_i}\nabla_i U(\{\mathbf{R}(t)\}) \quad (10b)$$

$$v_i(t + dt) = v_i(t + dt/2) - \frac{dt}{2m_i}\nabla_i U(\{\mathbf{R}(t + dt)\}) \quad (10c)$$

These equations impose that the total energy, linear and angular momenta are maintained during ZPAD simulations.

### 2.3.2. Thermalization

HNDs are known to have a temperature $T \approx 0.37$ K [1] and equilibration in the canonical ensemble is required to prepare the HND at a realistic temperature. Two methods are proposed to perform an equilibration at a constant temperature T in DynHeMat:

* Standard velocity rescaling where velocities are adjusted on the fly: The ith component of the new velocity vector is given by $v_{i,new} = \sqrt{\frac{E_T}{E_{kin}}} v_i$ where $E_{kin}$ is the current kinetic energy of the system and $E_T$ is the target thermal energy at temperature T.

* Velocity rescaling with canonical sampling [54]: This approach uses the same formula as the standard velocity rescaling method with a kinetic energy $E'_{kin}$ which evolves stochastically in such a way that $E'_{kin} = E_{kin} + dE_{kin}$ with

$$dE_{kin} = \frac{(E_T - E_{kin})dt}{\tau} + 2\sqrt{\frac{E_{kin}E_T}{N_f}} \frac{dW}{\sqrt{\tau}} \quad (11)$$

where $dt$ is the integrator time step, $\tau$ is a time constant for velocity rescaling, $N_f = 3(N_{tot} - 6)$ is the number of degrees of freedom and $dW$ is a Wiener noise represented by a random number normally distributed.

After proper thermalization, the temperature should barely fluctuate about the expected reference temperature during constant-energy ZPAD simulations.

## 2.4. MPI parallelization and benchmarks

A typical ZPAD simulation gathers $N_{traj}$ trajectories whose calculation is distributed over $N_{core}$ CPU cores. Indeed, ZPAD simulations become highly time consuming for large HNDs and the calculation of most theoretical expressions involving He atoms have been parallelized
with MPI (Message Passing Interface) to speed up the calculations. This includes the calculation of He-He interactions and their derivatives, He-impurity interactions and their derivatives, equations of motion of He atoms, and radial distribution functions between He atoms (see Section 3.3), for each trajectory. Short single-trajectory simulations have been achieved to estimate the efficiency $E_j$ of DynHeMat on parallel architectures:

$$E_j = \frac{i \times T_i}{j \times T_j} \qquad (12)$$

where $T_i$ (resp. $T_j$) is the elapsed time of a calculation performed on i (resp. j) CPU cores. Four ZPAD simulations have been performed as benchmarks:

* He$_{1000}$ modeled by the mPL He-He pseudopotential with a time step of 1 fs and a simulation time of 10 ps.

* He$_{1000}$ modeled by the IPTW He-He pseudopotential with a time step of 1 fs and a simulation time of 10 ps.

* He$_{79989}$ modeled by the mPL He-He pseudopotential with a time step of 1 fs and a simulation time of 1 ps.

* He$_{79989}$ modeled by the IPTW He-He pseudopotential with a time step of 1 fs and a simulation time of 1 ps.

DynHeMat efficiencies computed after distributing these simulations over 1 to 100 CPU cores are illustrated in Fig. 1. Simulations on He$_{1000}$ do not strongly benefit from MPI parallelization and the efficiency drops below 0.5 as soon as the number of CPU cores exceeds about 25 and 40 when the mPL and IPTW pseudopotentials are used, respectively. ZPAD simulations based on the IPTW pseudopotential are systematically more efficient because the analytical form of this pseudopotential is more complex and corresponding simulations are therefore more time consuming. Since simulations on He$_{79989}$ were carried out for shorter times, the influence of the reading and treatment of input data by MPI nodes is no longer negligible, especially for large numbers of CPU cores, and two curves are supplied to characterize the efficiency of ZPAD simulations on He$_{79989}$. The first one (solid line with open symbols) depicts the efficiency of the full ZPAD simulation and the second one (dashed line) is obtained by subtracting the time needed to perform a simulation with no dynamics. The efficiency curve of a realistic simulation (eg, $t_{max}$ = 100 ps) should lie in between these two curves, though presumably closer to the dashed one. An efficiency above 0.5 is ensured for high numbers of CPU cores (≥ 80) which confirms that MPI parallelization unleashes its full power for large HNDs.

## 3. Code instructions

### 3.1. Compilation

The DynHeMat program is written in fortran 90 and was parallelized with MPI. It was especially tested on x86 linux supercomputers (Intel Cascade Lake processors) and ARM architectures (Apple M1 Pro), compiled with the mpif90 command linked to standard f90 compilers (ifort or gfortran) and MPI libraries (IntelMPI or OpenMPI). It is worth noting that the option "-fallow-argument-mismatch" should be added if a compilation with gfortran is planned. This option will convert some improper MPI-related errors into warnings and enable the compilation: this is a known flaw in some gfortran versions. The user is provided with a makefile to compile DynHeMat. After modifying compiler options and the path to the mpif90 compiler, the compilation is simply achieved by running the command "make" in the source directory. By default, the executable will be called "DynHeMat".

### 3.2. Input files

A ZPAD simulation only requires two input files: a configuration file, which contains initial atomic positions and velocities, and a setup file based on a keyword system freely adapted from the last version of the MCMC² program distributed under GPL version 3 [55].

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 1. Efficiency of short ZPAD simulations of He$_{1000}$ (a) and He$_{79989}$ (b) achieved on the CSL partition of the Jean Zay supercomputer from the IDRIS French computer center with the mPL and IPTW He-He pseudopotentials. Dashed curves in figure b are obtained by subtracting the elapsed time of a simulation with no dynamics, whose duration is therefore mainly due to the reading and treatment of input data.</figcaption></figure>

#### 3.2.1. Input configuration file

Input configuration files use a XYZ format, namely:

Line 1: Number of atoms, $N_{tot}$.

Line 2: Any comment on the system (may be left blank).

Lines 3 to $N_{tot}$ + 2: atomic symbols (column 1), atomic cartesian coordinates (x: column 2, y: column 3, z: column 4), atomic cartesian velocity components ($v_x$: column 5, $v_y$: column 6, $v_z$: column 7).

Examples of input configuration files are given in the ZPAD_DB databank presented in Section 3.4 and in the test cases discussed in Section 3.5. Atomic symbols should follow the nomenclature detailed in Table 1 because the values of atomic masses are assigned on the basis of atomic symbols appearing in the configuration file. Initial velocities are optional because they can be generated randomly in DynHeMat (see Section 2.2.2). The position and velocity vector components of impurities (ie, rare-gas or alkali atoms) must be placed on top of the configuration file before the position and velocity vector components of He atoms. There is no internal restriction on the number of impurities or He atoms but mixing of unlike impurities is not currently allowed in DynHeMat.
<table><caption>**Table 2**
Setup variables with predefined set of values in DynHeMat.</caption><thead><tr><th>Keyword</th><th>Subkeyword</th><th>Variable</th><th>Values & definitions</th></tr></thead><tbody><tr><td>CONFIG</td><td>UNIT</td><td>input_unit</td><td>0: Atomic units<br>1: Positions in Å and velocities in Å/fs<br>2: SI units</td></tr><tr><td>ENERGY</td><td>UNIT</td><td>output_unit</td><td>0: Atomic unit (Hartree, Ha)<br>1: Kelvin (K)<br>2: Wavenumber (cm<sup>-1</sup>)<br>3: Electron-volt (eV)</td></tr><tr><td></td><td>HE2</td><td>ipotHe2</td><td>0: TT2 potential [46]<br>1: IPTW potential [40]<br>2: PL potential [41]<br>3: mPL potential [42]</td></tr><tr><td></td><td>IMPURITY</td><td>ipot</td><td>0: No impurity<br>1: TT2 He-Ne and Ne-Ne potentials<br>2: TT2 He-Ar and Ar-Ar potentials<br>3: TT2 He-Kr and Kr-Kr potentials<br>4: TT2 He-Xe and Xe-Xe potentials<br>5: TT2 He-Rn and Rn-Rn potentials<br>6: He-Ne pseudopotential from IPTW potential and TT2 Ne-Ne potential<br>7: He-Ar pseudopotential from IPTW potential and TT2 Ar-Ar potential<br>8: He-Kr pseudopotential from IPTW potential and TT2 Kr-Kr potential<br>9: He-Xe pseudopotential from IPTW potential and TT2 Xe-Xe potential<br>10: He-Rn pseudopotential from IPTW potential and TT2 Rn-Rn potential<br>11: He-Ne pseudopotential from PL potential and TT2 Ne-Ne potential<br>12: He-Ar pseudopotential from PL potential and TT2 Ar-Ar potential<br>13: He-Kr pseudopotential from PL potential and TT2 Kr-Kr potential<br>14: He-Xe pseudopotential from PL potential and TT2 Xe-Xe potential<br>15: He-Rn pseudopotential from PL potential and TT2 Rn-Rn potential<br>16: He-Ne pseudopotential from mPL potential and TT2 Ne-Ne potential<br>17: He-Ar pseudopotential from mPL potential and TT2 Ar-Ar potential<br>18: He-Kr pseudopotential from mPL potential and TT2 Kr-Kr potential<br>19: He-Xe pseudopotential from mPL potential and TT2 Xe-Xe potential<br>20: He-Rn pseudopotential from mPL potential and TT2 Rn-Rn potential<br>21: Li-He potential [47] and SMATB model for Li<sub>n</sub> [48]<br>22: Na-He potential [47] and SMATB model for Na<sub>n</sub> [48]<br>23: K-He potential [47] and SMATB model for K<sub>n</sub> [48]<br>24: Rb-He potential [47] and SMATB model for Rb<sub>n</sub> [48]<br>25: Cs-He potential [47] and SMATB model for Cs<sub>n</sub> [48]<br>26: Na-He potential [47] and DDTB model for Na<sub>n</sub> [49]</td></tr><tr><td>THERMALIZATION</td><td>METHOD</td><td>itherm</td><td>0: No thermalization<br>1: Velocity rescaling<br>2: Velocity rescaling with canonical sampling [54]</td></tr></tbody></table>

### 3.2.2. Setup file

The setup file enables the user to define a number of parameters to activate or deactivate functionalities of ZPAD simulations. A list of the ten keywords defined in DynHeMat is provided below in alphabetical order, together with the hardcoded default values for each parameter. Keywords are written in italic capital letters, subkeywords in capital letters, and variable names in italic. Some variables can only take values among a predefined set of values which is reported in Table 2.

* **CONFIG** `FILE_ROOT root_iniconf UNIT input_unit V0_SETTING lvgen Tref L0_SETTING llgen Lref(1) Lref(2) Lref(3):` Setting parameters for initial conditions. `root_iniconf` is the first part of input configuration filenames, which is identical for all trajectories. The end of the input configuration filename is always `"_AAA.xyz"` where AAA is the trajectory number. `input_unit` is an integer to specify units used in the input configuration files (see Table 2). `lvgen` is a boolean to accept or decline the random generation of initial velocities. `Tref` is the reference temperature (in K) used to generate initial velocities. `llgen` is a boolean to accept or decline the alteration of velocities to fix the angular momentum. `Lref(1) Lref(2) Lref(3)` are the values of the target angular momentum components, $L_x$, $L_y$ and $L_z$. Default values: `root_iniconf = iniconf`, `input_unit = 0`, `lvgen = false`, `Tref = 0.4`, `llgen = false`, `Lref(1) Lref(2) Lref(3) = 0.0 0.0 0.0`.

* **CONSERVATION** `ETOT letot detotlim JTOT ljtot djtotlim COM lrpcom comlim`: Setting parameters to check the conservation of

total energy, total angular momentum and center-of-mass position and linear momentum. `letot`, `ljtot`, and `lrpcom` are booleans to accept or decline tests on conservation laws. `detotlim`, `djtotlim` and `comlim` are the corresponding thresholds (in atomic units) above which classical quantities are assumed not to be maintained.

Default values: `letot = true`, `detotlim = 10⁻⁵`, `ljtot = false`, `djtotlim = 10⁻³`, `lrpcom = false`, `comlim = 10⁻¹⁰`.

* **ENERGY** `UNIT output_unit HE2 ipotHe2 IMPURITY ipot CUTOFF_HE rcut_HeHe rcut_HeX`: Setting parameters for output energies and potential-energy surfaces. `output_unit` is an integer to define the unit of energies in output files (see Table 2). `ipotHe2` is an integer to define the He-He pseudopotential to be used (see Table 2). `ipot` is an integer to define He-impurity and impurity-impurity models to be used (see Table 2). `rcut_HeHe` and `rcut_HeX` are cutoff distances (in Å) to calculate He-He and He-impurity interactions, respectively. Default values: `output_unit = 1`, `ipotHe2 = 1`, `ipot = 0`, `rcut_HeHe = 40.0`, `rcut_HeX = 10³`.

* **EQUILIBRATION** `T0 VVI_eq(1) TASTEP VVI_eq(2) TMAX VVI_eq(3) PRINT print_every_eq`: Setting parameters for equilibration runs. `VVI_eq(1)`, `VVI_eq(2)` and `VVI_eq(3)` are the time origin, the time step, and the final time of the dynamics (in ps), respectively. `print_every_eq` is the frequency (in number of time steps) between data printings. Default values: `VVI_eq(1) = 0.0`, `VVI_eq(2) = 10⁻³`, `VVI_eq(3) = 0.0`, `print_every_eq = 10`.
<table><caption>**Table 3**
Configuration files stored in the three directories of the ZPAD_DB databank.</caption><thead><tr><th>Directory IPTW/</th><th>Directory PL/</th><th>Directory mPL/</th></tr></thead><tbody><tr><td>He1000_IPTW.xyz</td><td>He1000_PL.xyz</td><td>He1000_mPL.xyz</td></tr><tr><td>He70000_IPTW.xyz</td><td>He5000_PL.xyz</td><td>He5000_mPL.xyz</td></tr><tr><td>He79989_IPTW.xyz</td><td>He70000_PL.xyz</td><td>He10000_mPL.xyz</td></tr><tr><td rowspan="8">He89986_IPTW.xyz</td><td>He79989_PL.xyz</td><td>He19998_mPL.xyz</td></tr><tr><td rowspan="7">He89986_PL.xyz</td><td>He29997_mPL.xyz</td><td>He39996_mPL.xyz</td></tr><tr><td>He59992_mPL.xyz</td><td>He79996_mPL.xyz</td></tr><tr><td>He70000_mPL.xyz</td><td>He79989_mPL.xyz</td></tr><tr><td>He89986_mPL.xyz</td><td>He10000_mPL.xyz</td></tr><tr><td></td><td></td></tr><tr><td></td><td></td></tr><tr><td></td><td></td></tr></tbody></table>

doped HNDs, by replacing inner helium atoms by helophilic species or depositing heliophobic species at the HND surface. Any user can also use the DynHeMat program to produce large HNDs by merging two smaller HNDs whose positions and velocities are provided in configuration files available in the ZPAD_DB databank. In such a case, it is recommended to merge HNDs based on the mPL He-He pseudopotential because it is more fluid and the merging process will be faster. HNDs based on deeper He-He pseudopotentials could then be obtained by equilibrating the subsequent HND with the desired pseudopotential.

### 3.5. Test cases

Five test cases, which cover a sample of current program applications, are made available to check the proper functioning of DynHeMat. They encompass three collision processes (test cases 1 to 3) and two simulations related to the formation of He<sub>10000</sub> (test cases 4 and 5). In test cases 1 to 3, the equilibration of the target HND was achieved before simulating the collision process in a production run. On the contrary, test cases 4 and 5 are typical equilibration runs dedicated to the formation of large HNDs. Benchmark results for these five test cases were computed on an iMac computer (MacOS Ventura version 13.5.1) with a gfortran compiler and an OpenMPI library taken from the MacPorts project [57]. Data on test cases are gathered in a directory called TEST Cases/ which contains five subdirectories, one by test case, called TEST CASE A where A = 1 to 5 is the test case number. Each of these subdirectories contains the XYZ input configuration files required to run the test case, the appropriate DynHeMat setup file, a "README.txt" file which details the content of the setup file and how to run the test case, and a directory called REF/ which gathers benchmark results. It is worth noting that the maintenance of total energy, linear and angular momenta as well as the elapsed time of the simulation can be checked in the standard output files called DynHeMat_AAA.out (see Section 3.3). A brief description of test cases is provided below.

* Test case 1: 4 independent trajectories distributed over 4 CPU cores (one CPU core per trajectory) to simulate the collision between one Ar atom ($v_0 = 500 \text{ m s}^{-1}$) and ArHe<sub>10000</sub> with zero impact parameter.

* Test case 2: A single trajectory distributed over 4 CPU cores to simulate the collision between one Xe atom ($v_0 = 276 \text{ m s}^{-1}$) and Xe<sub>4</sub>He<sub>10000</sub> with zero impact parameter.

* Test case 3: 2 independent trajectories distributed over 4 CPU cores (2 CPU cores per trajectory) to simulate the collision between one Cs atom ($v_0 = 400 \text{ m s}^{-1}$) and He<sub>10000</sub> with zero impact parameter.

* Test case 4: A single trajectory distributed over 8 CPU cores to simulate the merging of two He<sub>5000</sub> clusters at $T = 0.47 \text{ K}$.

* Test case 5: A single trajectory distributed over 8 CPU cores to simulate the thermalization of He<sub>10000</sub> at $T = 0.37 \text{ K}$.

Most of these calculations use the mPL He-He pseudopotential because it was the most commonly used He-He pseudopotential in recent ZPAD simulations on collision processes [42,43] and it is more practical

for coalescence processes because HNDs formed with this pseudopotential are more fluid.

### 4. Applications

The DynHeMat program enables the user to carry out ZPAD simulations to investigate the dynamics of pure and doped HNDs where the impurity can be a rare-gas or alkali atom or cluster. Recent publications [28,42,43] focused on small HNDs ($N = 1000$) to conduct systematic studies on different collisions as a function of the projectile speed and impact parameter for a significant number of trajectories (up to 100) and to compare the subsequent results with those from more sophisticated theoretical approaches, namely TDDFT and ring polymer molecular dynamics (RPMD). Some results of these publications are synthesized in Sections 4.1 and 4.2 before briefly discussing the formation of large HNDs and its possible application to the study of cluster coagulation or submergence in liquid helium.

#### 4.1. Coagulation of rare-gas atoms within HNDs

In experiments on doped HNDs, the chamber after the helium source is a pickup cell where the impurity of interest is expected to be captured by HNDs, then raising the question of the propensity of these impurities to remain at the HND surface or to submerge in helium. Rare-gas atoms were one of the first systems whose capture and coagulation in liquid helium has been investigated [58,59] but dynamical simulations remain sparse due to the large size of HNDs and the quantum nature of helium atoms. In this context, the coagulation of rare-gas atoms ($RG = Ne, Ar, Kr, Xe$ and $Rn$) in pure and doped He<sub>10000</sub> has been tackled by ZPAD simulations using the mPL, PL, and IPTW He-He pseudopotentials [42]. These calculations were called ZPAD-mPL, ZPAD-PL, and ZPAD-IPTW (or simply ZPAD), respectively. $RG + RG_nHe_{10000}$ ($n = 1, 2$ or $4$) collisions have been considered for various projectile speeds and two impact parameters, $b = 0$ and $10 \AA$, but meaningful sample results are only discussed here. It is worth noting that it is good practice to prepare such simulations by equilibrating $RG_nHe_{10000}$ independently with DynHeMat before adding the projectile position and velocity vector components on top of the input configuration file.

The main goal of this work was to assess the ability of the rare-gas projectile to penetrate inside the droplet to form bound structures (ie, structures close to the gas-phase structure), saddle structures (ie, compact structures different from the gas-phase structure), or unbound structures. The analysis of $Ar(v_0 = 500 \text{ m s}^{-1}) + Ar_4He_{10000}$ collisions proved that ZPAD-mPL simulations yield results in qualitative agreement with TDDFT and RPMD simulations, as illustrated in Fig. 2a, while other ZPAD-like methods underestimate the number of compact structures. When the speed of the projectile increases up to $800 \text{ m s}^{-1}$, the amount of bound and saddle structures rises and that of unbound structures drastically decreases to reach no more than ~20% (see Fig. 2b). At higher velocities, the abundance of unbound structures starts to stabilize while that of bound structures increases more sharply at the expense of saddle structures. Interestingly, ZPAD-PL and ZPAD-mPL simulations lead to the same trends in Fig. 2b.

Argon projectile speeds close to $1000 \text{ m s}^{-1}$ are not realistic and collisions modeled for other rare-gas atoms were performed by providing the projectile with a speed which corresponds to the kinetic energy of an Ar projectile moving at $v_0 = 500 \text{ m s}^{-1}$. These calculations proved that the coagulation of rare-gas atoms is much more efficient for heavy rare-gas atoms, Xe and Rn, than for light ones, Ne and Ar. This conclusion can be drawn, for instance, from Fig. 3a where abundances of bound and unbound structures were estimated after performing ZPAD-mPL simulations to model $RG + RG_nHe_{10000}$ collisions. A closer inspection of kinetic energy distributions of helium atoms ejected during $Ne(v_0 = 707 \text{ m s}^{-1}) + NeHe_{10000}$ and $Rn(v_0 = 212 \text{ m s}^{-1}) + RnHe_{10000}$ collisional processes revealed that the faster Ne projectiles favor the evaporation of high-energy He atoms during the first picoseconds after the impact with the
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Abundances of bound (black), saddle (blue) and unbound (red) structures for collisions Ar($v_0$,b) + Ar$_n$He$_{1000}$. (a) Collisions modeled by TDDFT, RPMD and ZPAD-like methods with $v_0$ = 500 m s$^{-1}$ and b = 0 or 10 Å. (b) Collisions modeled by ZPAD-mPL (solid curves) or ZPAD-PL (dashed curves) as a function of the initial projectile speed at b = 0. Reproduced from "B. Fixot, E. Louaas, D. A. Bonhommeau, J. Chem. Phys. 161, 044307 (2024)", with the permission of AIP Publishing. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

HND but only 7 He atoms are ejected after 200 ps. On the contrary, the collision with the slower Rn projectile yields the ejection of a greater number of lower-energy He atoms and helium evaporation is not ended at t = 200 ps.

These results were obtained for relatively small HNDs ($N$ = 1000) and a direct extension of this work would be to assess the influence of the HND size on the final outcome of RG + RG$_n$He$_N$ collisions. In particular, ZPAD-mPL simulations were found in qualitative agreement with TDDFT simulations for Ar + Ar$_n$He$_{1000}$ collisions whereas superfluidity is not taken into account in ZPAD-like approaches. It is however unclear whether the inclusion of superfluidity in a phenomenological manner would be mandatory for calculations on large HNDs for at least two reasons: (i) coagulation does not systematically occur in superfluid HNDs modeled by TDDFT, and (ii) we have shown [42] that unbound structures in He$_{1000}$ were characterized by projectiles roaming in a small region of the HND at velocities in between 10 and 30 m s$^{-1}$. Such speeds are sufficient to favor coagulation provided that the projectile velocity

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** (a) Abundances of bound (black) and unbound (red) structures for collisions RG + RGHe$_{1000}$ at b = 0. (b) Average kinetic energy and number of ejected He atoms for collisions Ne + NeHe$_{1000}$ (black) and Rn + RnHe$_{1000}$ (red) at b = 0. Reproduced from "B. Fixot, E. Louaas, D. A. Bonhommeau, J. Chem. Phys. 161, 044307 (2024)", with the permission of AIP Publishing. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

vector points in the direction of the embedded dopant, which was not the case in our ZPAD-like simulations, then accounting for the existence of unbound structures.

## 4.2. Dynamics of alkali atoms and clusters in HNDs

Unlike rare-gas atoms, alkali atoms are heliophobic and should remain at the HND surface upon soft landing. TDDFT calculations predict that Cs atoms should be captured by the HND surface if their initial speed remains below 75 m s$^{-1}$. Beyond this value, Cs atoms bounce on the HND surface or pierce the droplet and move away [29]. In contrast, RPMD and former ZPAD simulations predicted soft capture of the Cs atom at initial speeds up to 500 m s$^{-1}$ (higher speeds were not tested with these methods) and impact parameters b = 0 and 10 Å, as depicted in Fig. 4a. ZPAD-mPL simulations yield an intermediate picture: soft Cs capture is a dominant process when $v_0 \lesssim 300$ m s$^{-1}$ but Cs capture remains possible up to 600 m s$^{-1}$ after a transient submersion in helium which extends on the nanosecond time scale (see Fig. 4a and b). At higher velocities, the Cs atom tends to pierce the HND which stimulates the ejection of
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Abundance of Cs captures, bounces, transient submersions and escapes for collisions Cs + He$_{1000}$. (a) Collisions modeled by TDDFT, RPMD and ZPAD-like methods with $v_0 = 500$ m s$^{-1}$ and $b = 0$ or $10$ Å. (b) Collisions modeled by ZPAD-mPL as a function of the initial projectile speed at $b = 0$. Reproduced from "D. A. Bonhommeau, J. Chem. Phys. 161, 184302 (2024)", with the permission of AIP Publishing.</figcaption></figure>

He atoms. Molecular sketches for all these cases are available in our reference publication on this topic [43].

Despite the use of a He-He pseudopotential with a shallow potential-energy well, Cs atoms were found to maintain their heliophobicity. Unpublished ZPAD-mPL simulations performed on other alkali atoms (ie, Li, Na, K and Rb) at low initial projectile speeds also yielded soft capture of the projectile on the HND surface, although the corresponding He-alkali interactions are deeper than the He-Cs interaction, which makes these atoms less heliophobic than cesium [47]. A more challenging problem concerns the modeling of the submersion of alkali clusters in liquid helium. The critical size above which these clusters submerge in helium was estimated experimentally to be $n_c = 20$ for Na [26], $n_c \approx 80$ for K [44], and $n_c \gtrsim 100$ for Rb [45], in agreement with theoretical predictions based on an analysis of the balance between energy terms in doped HNDs [27]. The mPL, PL and IPTW He-He pseudopotentials were originally devised from calculations performed on relatively small HNDs ($N < 10^4$) and their ability to model properly the submersion of alkali clusters in HNDs of experimental size ($N \lesssim 10^5$) still needs to be assessed. DynHemat enables the user to investigate the submersion of alkali clusters in HNDs and such a work is currently in progress in our group. One limitation of this kind of simulations is the need for the formation of large HNDs, representative of experimental sizes, by suc-

<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Coalescence of two He$_{5000}$ clusters modeled with the mPL pseudopotential at a temperature of 0.47 K. Molecular sketches were generated with the VMD software [56].</figcaption></figure>

cessive coalescences of smaller droplets. An example of coalescence between two He$_{5000}$ clusters is presented in Fig. 5. In this particular case, the formation of a nearly spherical HND requires to equilibrate the system for 600 ps but the formation of larger HNDs may require ZPAD simulations extending on the nanosecond time scale, which accounts for the MPI parallelization of DynHeMat.

To avoid DynHeMat users to construct structures of large HNDs from scratch, He$_N$ structures for sizes in between $N = 1000$ and $N = 89986$ are supplied in ZPAD_DB (see Section 3.4). These structures can be used to build structures of doped HNDs by replacing core He atoms by the expected dopant, before equilibrating the new structure. It is also worth noting that the formation of large HNDs by successive coalescence processes may excite some vibrational modes of the HND and long equilibration may be required to damp the oscillations of the droplet radius. In ZPAD simulations based on the IPTW pseudopotential [28] the merging of two He$_{500}$ clusters was very slow and no spherical HND was formed after 200 ps while RPMD simulations yielded an almost spherical HND. In this respect, results from ZPAD simulations based on He-He pseudopotentials with a more pronounced fluid character (eg, mPL pseudopotential) now resemble RPMD outcomes with spherical HNDs produced after a few hundreds picoseconds. However, both methods lead to a decrease
of the HND quadrupole moment due to the steady gain in sphericity of the HND during the merging process whereas quadrupole deformations yielding a slow oscillatory motion with a period of 160 ps are observed in TDDFT simulations which include superfluidity.

Last but not least, the internal structure and keyword system of Dyn-
HeMat were thought to facilitate future extensions of the program, such
as the implementation of He-He pseudopotentials maybe more appro-
priate for submersion studies on large HNDs or that of additional He-
impurity and impurity-impurity potential-energy curves to tackle the
zero-point averaged dynamics of homogeneous or heterogeneous clus-
ters in helium.

CRediT authorship contribution statement

**David A. Bonhommeau:** Conceptualization, Data curation, Fund-
ing acquisition, Investigation, Methodology, Project administration, Re-
sources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.

Data availability

Data will be made available on request.

Declaration of competing interest

The authors declare the following financial interests/personal rela-
tionships which may be considered as potential competing interests:
David A. Bonhommeau reports equipment, drugs, or supplies was pro-
vided by Grand Équipement National De Calcul Intensif. If there are
other authors, they declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

Acknowledgements

This work was granted access to the HPC resources of IDRIS un-
der allocations 2023-AD010814624, 2024-AD010814624R1, and 2025-
AD010814624R2 made by GENCI.

# Appendix A. Extension of DynHeMat to other He-He pseudopotentials and impurities

He-He pseudopotentials used in ZPAD simulations have been con-
structed by following two main approaches. Pseudopotentials like the
IPTW pseudopotential implemented in DynHeMat are constructed it-
eratively by performing MD simulations on small HNDs, extracting
the helium pair correlation function from these simulations, building
the average helium radial potential to solve the corresponding radial
Schrödinger equation, and finally convolving the He-He pair potential
with the squared helium wave function to get a He-He pseudopotential.
This procedure is repeated until convergence by updating the He-He
potential used in MD simulations with the pseudopotential computed at
the end of the last iteration. This method provided a number of He-He
pseudopotentials for small pure and doped HNDs (N ≤ 1000) [37–40]
but the overall procedure is time consuming. Moreover, convergence is
generally ensured by truncating the helium wave function at distances
between 1.5 Å and 1.6 Å and data points are poorly adjusted on a Morse
potential, which imposed to consider a more sophisticated analytical
form for these pseudopotentials. The second approach, which was used
to build the mPL pseudopotential, consists in getting general properties
of a target HND (eg, total energy and radius) from He-DFT calculations
and adjusting the parameters of a Morse potential to reproduce these
properties in ZPAD simulations.

New He-He pseudopotentials and their derivatives can be included
in DynHeMat by adapting subroutines available in lib_Rg-pot.f90, which
gathers interaction potentials and related derivatives for all rare gases

including helium, provided that the pseudopotential is fitted on an ana-
lytical expression implemented in DynHeMat. In particular, subroutines
POT_HE2_MPL and DPOT_HE2_MPL calculate the sum of potential en-
ergies and their derivatives for the mPL pseudopotential, respectively.
The selection of a new potential from the setup file can be implemented
by (i) extending the range of allowed values for the ipotHe2 integer,
which is dedicated to the selection of He-He interaction potentials in
lib_config.f90, and (ii) reflecting this change in subroutines POT_HEN,
DPOT_HEN, DPOT_XHEN and SURFACE, which use ipotHe2 as input in
SELECT CASE statements. The three first subroutines are included in
lib_Rg-pot.f90 and the last one in libsurface.f90.

The implementation of new impurities in DynHeMat implies to ex-
tend the range of values of the *ipot* integer in lib_config.f90 and reflect
this change in the corresponding SELECT CASE statement of subrou-
tine DPOT_XHEN. The definition of interaction potentials and related
derivatives will then strongly depend on the nature of the impurity un-
der scrutiny. Note that the whole program uses MPI blocking routines to
speed up ZPAD simulations and parallelization of fortran subroutines for
the new impurity can follow a similar strategy by adapting subroutines
already present in DynHeMat.

References

[1] J. Harms, M. Hartmann, B. Sartakov, J.P. Toennies, A.F. Vilesov, Rotational structure of the IR spectra of single SF₆ molecules in liquid ⁴He and ³He droplets, J. Mol. Spectrosc. 185 (1997) 204–206. https://doi.org/10.1006/jmsp.1997.7343

[2] S. Grebenev, J.P. Toennies, A.F. Vilesov, Superfluidity within a small helium-4 cluster: the microscopic andronikashvili experiment, Science 279 (1998) 2083-2086. https://doi.org/10.1126/science.279.5359.2083

[3] M. Hartmann, R.E. Miller, J.P. Toennies, A. Vilesov, Rotationally resolved spectroscopy of SF₆ in liquid helium clusters: a molecular probe of cluster temperature, Phys. Rev. Lett. 75 (1995) 1566-1569. https://doi.org/10.1103/PhysRevLett.75.1566

[4] F. Stienkemeier, J. Higgins, W.E. Ernst, G. Scoles, Laser spectroscopy of alkali-doped helium clusters, Phys. Rev. Lett. 74 (1995) 3592-3595.

[5] J.P. Toennies, A.F. Vilesov, Spectroscopy of atoms and molecules in liquid helium, Ann. Rev. Phys. Chem. 49 (1998) 1-41. https://doi.org/10.1146/annurev.physchem.49.1.1

[6] W.K. Lewis, B.E. Applegate, J. Sztáray, B. Sztáray, T. Baer, R.J. Bemish, R.E. Miller, Electron impact ionization in helium nanodroplets: controlling fragmentation by active cooling of molecular ions, J. Am. Chem. Soc. 126 (2004) 11283-11292. https://doi.org/10.1021/ja030653q

[7] W.K. Lewis, R.J. Bemish, R.E. Miller, Fragmentation of HCN in optically selected mass spectrometry: nonthermal ion cooling in helium nanodroplets, J. Chem. Phys. 123 (2005) 141103. https://doi.org/10.1063/1.2046672

[8] A. Braun, M. Drabbels, Imaging the translational dynamics of CF₃ in liquid helium droplets, Phys. Rev. Lett. 93 (2004) 253401. https://doi.org/10.1103/PhysRevLett.93.253401

[9] K. Nauta, R.E. Miller, Nonequilibrium self-assembly of long-chains of polar molecules in superfluid helium, Science 283 (1999) 1895-1897. https://doi.org/10.1126/science.283.5409.18

[10] K. Nauta, R.E. Miller, Formation of cyclic water hexamer in liquid helium: the smallest piece of ice, Science 287 (2000) 293-295. https://doi.org/10.1126/science.287.5451.293

[11] E. Latimer, D. Spence, C. Feng, A. Boatwright, A.M. Ellis, S. Yang, Preparation of ultrathin nanowires using superfluid helium droplets, Nano Lett. 14 (2014) 2902-2906. https://doi.org/10.1021/nl500946u

[12] G. Haberfehler, P. Thaler, D. Knez, A. Volk, F. Hofer, W.E. Ernst, G. Kothleitner, Formation of bimetallic clusters in superfluid helium nanodroplets analysed by atomic resolution electron tomography, Nat. Commun. 6 (2015) 8779. https://doi.org/10.1038/ncomms9779

[13] D. Pentlechner, J.H. Nielsen, A. Slenczka, K. Mølmer, H. Stapelfeldt, Impulsive laser induced alignment of molecules dissolved in helium nanodroplets, Phys. Rev. Lett. 110 (2013) 093002. https://doi.org/10.1103/PhysRevLett.110.093002

[14] J. Qiang, L. Zhou, P. Lu, K. Lin, Y. Ma, S. Pan, C. Lu, W. Jiang, F. Sun, W. Zhang, H. Li, X. Gong, I.S. Averbukh, Y. Prior, C.A. Schouder, H. Stapelfeldt, I.N. Cherepanov, M. Lemeshko, W. Jäger, J. Wu, Femtosecond rotational dynamics of D₂ molecules in superfluid helium nanodroplets, Phys. Rev. Lett. 128 (2022) 243201. https://doi.org/10.1103/PhysRevLett.128.243201

[15] J. Qiang, L. Zhou, Y. Peng, C. Yu, P. Lu, S. Pan, C. Lu, G. Chen, R. Lu, W. Zhang, J. Wu, Femtosecond collisional dissipation of vibrating D₂⁺ in helium nanodroplets, Phys. Rev. Lett. 132 (2024) 103201. https://doi.org/10.1103/PhysRevLett.132.103201

[16] L. Kranabetter, H.H. Kristensen, A. Ghazaryan, C.A. Schouder, A.S. Chatterley, P. Janssen, F. Jensen, R.E. Zillich, M. Lemeshko, H. Stapelfeldt, Nonadiabatic laser-induced alignment dynamics of molecules on a surface, Phys. Rev. Lett. 131 (2023) 053201. https://doi.org/10.1103/PhysRevLett.131.053201

[17] B. Thaler, M. Meyer, P. Heim, M. Koch, Long-lived nuclear coherences inside helium nanodroplets, Phys. Rev. Lett. 124 (2020) 115301. https://doi.org/10.1103/PhysRevLett.124.115301
[18] L.F. Gomez, K.R. Ferguson, J.P. Cryan, C. Bacellar, R.M.P. Tanyag, C. Jones, S. Schorb, D. Anielski, A. Belkacem, C. Bernando, R. Boll, J. Bozek, S. Carron, G. Chen, T. Delmas, L. Englert, S.W. Epp, B. Erk, L. Foucar, R. Hartmann, A. Hexemer, M. Huth, J. Kwok, S.R. Leone, J.H.S. Ma, F.R.N.C. Maia, E. Malmerberg, S. Marchesini, D.M. Neumark, B. Poon, J. Prell, D. Rolles, B. Rudek, A. Rudenko, M. Seifrid, K.R. Siefermann, F.P. Sturm, M. Swiggers, J. Ullrich, F. Weise, P. Zwart, C. Bostedt, O. Gessner, A.F. Vilesov, Shapes and vortivities of superfluid helium nanodroplets, Science 345 (2014) 906–909. https://doi.org/10.1126/science.1252

[19] M.V. Rama Krishna, K.B. Whaley, Collective excitations of helium clusters, Phys. Rev. Lett. 64 (1990) 1126–1129. https://doi.org/10.1103/PhysRevLett.64.1126

[20] M.V. Rama Krishna, K.B. Whaley, Microscopic studies of collective spectra of quantum liquid clusters, J. Chem. Phys. 93 (1990) 746–759. https://doi.org/10.1063/1.459525

[21] D. Blume, M. Lewerenz, F. Huisken, M. Kaloudis, Vibrational frequency shift of HF in helium clusters: quantum simulation and experiment, J. Chem. Phys. 105 (1996) 8666–8683. https://doi.org/10.1063/1.472648

[22] C.A. Brindle, M.R. Prado, K.C. Janda, N. Halberstadt, M. Lewerenz, Structure and stability of (NeHe)n⁺: experiment and diffusion quantum Monte Carlo theory with "on the fly" electronic structure, J. Chem. Phys. 123 (2005) 064312. https://doi.org/10.1063/1.1982790

[23] T. Takayanagi, M. Shiga, Photodissociation of Cl₂ in helium clusters: an application of hybrid method of quantum wavepacket dynamics and path integral centroid molecular dynamics, Chem. Phys. Lett. 372 (2003) 90–96. https://doi.org/10.1016/S0009-2614(03)00368-3

[24] T. Takayanagi, M. Shiga, Theoretical study on photoexcitation dynamics of the K atom attached to helium clusters and the solvation structures of K⁺Heₐn exciplexes, Phys. Chem. Chem. Phys. 6 (2004) 3241–3247. https://doi.org/10.1039/B400277F

[25] F. Calvo, Submersion of sodium impurities in finite cryogenic droplets: a path-integral molecular dynamics study for ⁴He and para-H₂, Phys. Rev. B 95 (2017) 035429. https://doi.org/10.1103/PhysRevB.95.035429

[26] L. An der Lan, P. Bartl, C. Leidlmair, H. Schöbel, R. Jochum, S. Denifl, T.D. Märk, A.M. Ellis, P. Scheier, The submersion of sodium clusters in helium nanodroplets: identification of the surface → interior transition, J. Chem. Phys. 135 (2011) 044309. https://doi.org/10.1063/1.3610388

[27] C. Stark, V.V. Kresin, Critical sizes for the submersion of alkali clusters into liquid helium, Phys. Rev. B 81 (2010) 085401. https://doi.org/10.1103/PhysRevB.81.085401

[28] E. García-Alfonso, M. Barranco, D.A. Bonhommeau, N. Halberstadt, M. Pi, F. Calvo, Clustering, collision, and relaxation dynamics in pure and doped helium nanoclusters: density- vs particle-based approaches, J. Chem. Phys. 157 (2022) 014106. https://doi.org/10.1063/5.0091942

[29] A. Leal, D. Mateo, A. Hernando, M. Pi, M. Barranco, Capture of heliophobic atoms by ⁴He nanodroplets: the case of cesium, Phys. Chem. Chem. Phys. 16 (2014) 23206–23213. https://doi.org/10.1039/c4cp03297g

[30] F. Coppens, F. Ancilotto, M. Barranco, N. Halberstadt, M. Pi, Dynamics of impurity clustering in superfluid ⁴He nanodroplets, Phys. Chem. Chem. Phys. 21 (2019) 17423–17432. https://doi.org/10.1039/c9cp02789k

[31] S.H. Albrechtsen, C.A. Schouder, A.V.n. Muñoz, J.K. Christensen, C.E. Petersen, M. Pi, M. Barranco, H. Stapelfeldt, Observing the primary steps of ion solvation in helium droplets, Nature 623 (2023) 319–323. https://doi.org/10.1038/s41586-023-06593-5

[32] E. García-Alfonso, M. Barranco, N. Halberstadt, M. Pi, Time-resolved solvation of alkali ions in superfluid helium nanodroplets, J. Chem. Phys. 160 (2024) 164308. https://doi.org/10.1063/5.0205951

[33] S.H. Albrechtsen, J.K. Christensen, C.E. Petersen, C.A. Schouder, P.J. Carchi-Villalta, I. Sánchez-Pérez, M. Bartolomei, T. González-Lezana, F. Pirani, H. Stapelfeldt, Femtosecond-and-atom-resolved solvation dynamics of a Na⁺ ion in a helium nanodroplet, J. Chem. Phys. 162 (2025) 174309. https://doi.org/10.1063/5.0260588

[34] E. García-Alfonso, M. Barranco, M. Pi, N. Halberstadt, Time-resolved solvation of alkali ions in superfluid helium nanodroplets: theoretical simulation of a pump-probe study, J. Chem. Phys. 163 (2025) 144309. https://doi.org/10.1063/5.0291643

[35] M. Sterling, Z. Li, V.A. Apkarian, Simulations of quantum crystals by classical dynamics, J. Chem. Phys. 103 (1995) 5679–5683.

[36] P. Slavíček, P. Jungwirth, M. Lewerenz, N.H. Nahler, M. Fárník, U. Buck, Pickup and photodissociation of hydrogen halides in floppy neon clusters, J. Phys. Chem. A 107 (2003) 7743–7754. https://doi.org/10.1021/jp0357525

[37] D. Bonhommeau, P.T. Lake, Jr., C.L. Quiniou, M. Lewerenz, N. Halberstadt, Modeling the fragmentation dynamics of ionic clusters inside helium nanodroplets: the case of He₁₀₀Ne₄⁺, J. Chem. Phys. 126 (2007) 051104. https://doi.org/10.1063/1.2515225

[38] D. Bonhommeau, M. Lewerenz, N. Halberstadt, Fragmentation of ionized doped helium nanodroplets: theoretical evidence for a dopant ejection mechanism, J. Chem. Phys. 128 (2008) 054302. https://doi.org/10.1063/1.2823101

[39] N. Halberstadt, D.A. Bonhommeau, Fragmentation dynamics of ArₐHe₁₀₀ upon electron impact ionization: competition between ion ejection and trapping, J. Chem. Phys. 152 (2020) 234305. https://doi.org/10.1063/5.0009363

[40] D.A. Bonhommeau, Zero-point averaged dynamics of ArₐHe₁₀₀: dopant-size influence on potential-energy surfaces, mass spectra and fragmentation patterns, Chem. Phys. 550 (2021) 111307. https://doi.org/10.1016/j.chemphys.2021.111307

[41] R. Panzou, M. Lewerenz, Photofragmentation of I₂ molecules inside helium clusters: model calculations using a quantum effective potential, Mol. Phys. 119 (2021) e1977862. https://doi.org/10.1080/00268976.2021.1977862

[42] B. Fixot, E. Louaas, D.A. Bonhommeau, Collision of rare-gas atoms on helium nanodroplets: theoretical evidence for an efficient coagulation of heavy rare-gas atoms, J. Chem. Phys. 161 (2024) 044307. https://doi.org/10.1063/5.0220027

[43] D.A. Bonhommeau, Collision of cesium atoms on helium nanodroplets: unraveling mechanisms for surface capture at experimental velocities, J. Chem. Phys. 161 (2024) 184302. https://doi.org/10.1063/5.0231641

[44] L. An der Lan, P. Bartl, C. Leidlmair, H. Schöbel, S. Denifl, T.D. Märk, A.M. Ellis, P. Scheier, Submersion of potassium clusters in helium nanodroplets, Phys. Rev. B 85 (2012) 115414. https://doi.org/10.1103/PhysRevB.85.115414

[45] A. Schiller, P. Martini, E.E. J.A. Maalouf, P. Scheier, Submersion of rubidium clusters in helium nanodroplets, Eur. Phys. J. D 75 (2021) 124. https://doi.org/10.1140/s10053-021-00112-9

[46] X. Sheng, J.P. Toennies, K.T. Tang, Conformal analytical potential for all the rare gas dimers over the full range of internuclear distances, Phys. Rev. Lett. 125 (2020) 253402. https://doi.org/10.1103/PhysRevLett.125.253402

[47] S.H. Patil, Adiabatic potentials for alkali-inert gas systems in the ground state, J. Chem. Phys. 94 (1991) 8089. https://doi.org/10.1063/1.460091

[48] Y. Li, E. Blaisten-Barojas, D.A. Papaconstantopoulos, Structure and dynamics of alkali-metal clusters and fission of highly charged clusters, Phys. Rev. B 57 (1998) 15519–15532. https://doi.org/10.1103/PhysRevB.57.15519

[49] F. Calvo, F. Spiegelmann, Mechanisms of phase transitions in sodium clusters: from molecular to bulk behavior, J. Chem. Phys. 112 (2000) 2888–2908. https://doi.org/10.1063/1.480862

[50] LAPACK version 3.4.1. http://www.netlib.org/lapack/.

[51] J. Jellinek, D.H. Li, Separation of the energy of overall rotation in any N-body system, Phys. Rev. Lett. 62 (1989) 241–244. https://doi.org/10.1103/PhysRevLett.62.241

[52] D. Bonhommeau, A. Viel, N. Halberstadt, Dissociative ionization of neon clusters Neₙ, n = 3 to 14: a realistic multisurface dynamical study, J. Chem. Phys. 123 (2005) 054316. https://doi.org/10.1063/1.1953530

[53] M.A. Miller, D.A. Bonhommeau, C.J. Heard, Y. Shin, R. Spezia, M.P. Gaigeot, Structure and stability of charged clusters, J. Phys. Condens. Matter 24 (2012) 284130.

[54] G. Bussi, D. Donadio, M. Parrinello, Canonical sampling through velocity-rescaling, J. Chem. Phys. 126 (2007) 014101. https://doi.org/10.1063/1.2408420

[55] D.A. Bonhommeau, MCMC2 (Version 1.1.2): a Monte Carlo code for multiply charged clusters, Comput. Phys. Commun. 207 (2016) 533–535. https://doi.org/10.1016/j.cpc.2016.06.009

[56] W. Humphrey, VMD: visual molecular dynamics, J. Mol. Graph. 14 (1) (1996) 33–38. https://doi.org/10.1016/0263-7855(96)00018-5

[57] The MacPorts project official homepage, The MacPorts Project Official Homepage, (url)https://www.macports.org/.

[58] A. Scheidemann, J.P. Toennies, J.A. Northby, Capture of neon atoms by ⁴He clusters, Phys. Rev. Lett. 64 (1990) 1899–1902. https://doi.org/10.1103/PhysRevLett.64.1899

[59] M. Lewerenz, B. Schilling, J.P. Toennies, Successive capture and coagulation of atoms and molecules to small clusters in large liquid helium clusters, J. Chem. Phys. 102 (1995) 8191–8207. https://doi.org/10.1063/1.469231