In this work, we show that van der Waals molecules X–RG (where RG is the rare gas atom) may be created through direct three-body recombination collisions, i.e., X + RG + RG → X–RG + RG. In particular, the three-body recombination rate at temperatures relevant for buffer gas cell experiments is calculated via a classical trajectory method in hyperspherical coordinates [Pérez-Ríos *et al.*, J. Chem. Phys. **140**, 044307 (2014)]. As a result, it is found that the formation of van der Waals molecules in buffer gas cells (1 K ≲ *T* ≲ 10 K) is dominated by the long-range tail (distances larger than the LeRoy radius) of the X–RG interaction. For higher temperatures, the short-range region of the potential becomes more significant. Moreover, we notice that the rate of formation of van der Walls molecules is of the same order of the magnitude independent of the chemical properties of X. As a consequence, almost any X–RG molecule may be created and observed in a buffer gas cell under proper conditions.

## I. INTRODUCTION

When a three-body process leads to the formation of a molecule as a product state, A + A + A → A_{2} + A, it is labeled a three-body recombination process or a ternary association reaction. These three-body processes are relevant for a wide variety of systems in areas ranging from astrophysics to ultracold physics. In particular, three-body recombination of hydrogen is one of the essential processes to explain H_{2} formation in star-forming regions.^{1,2} In the field of ultracold physics, recent developments in laser technologies and cooling techniques have made it possible to gain a more in-depth insight into the significant role of three-body recombination processes in different phenomena, such as atomic loss processes in ultracold dilute gases^{3–9} and the formation and trapping of cold and ultracold molecules.^{10–14}

van der Waals (vdW) molecules consist of two atoms held together by the long-range dispersion interaction^{15} presenting binding energies ≲1 meV. Therefore, vdW molecules show the weakest gas-phase molecular bond in nature, except for ultra-long-range Rydberg molecules showing binding energies $\u223c4$ neV.^{65–70} The binding mechanism in vdW molecules relies on the compensation between the short-range repulsive interaction (due to the overlap of closed-shell orbitals) and the attractive −*C*_{6}/*r*^{6} vdW interaction, where the dispersion coefficient *C*_{6} depends on the polarizability of the interacting atoms. Interestingly enough, the study of vdW interactions provides crucial information necessary to investigate the formation and stability of gases, liquids, and materials such as vdW heterostructures and biopolymers,^{16–19} chemical reactions,^{18,20–24} and physical phenomena such as superfluidity of ^{4}He nanodroplets.^{25,26} In particular, investigating the properties of vdW molecules (as the simplest form of vdW complexes) containing a rare gas atom leads to a deeper understanding of the nature of bonding in rare gas crystals and of the dynamics of impurities interacting with dense rare gas vapors.^{27–30}

Despite the significance of vdW molecules in modern chemical physics, the community has been focused on its characterization rather than on revealing how they emerge in different scenarios.^{17,20,21,30–37} Recently, thanks to the development of buffer gas sources,^{38} it has been possible to investigate the formation of vdW molecules through three-body recombination processes.^{30,33–35,37,39,40} However, the field is still lacking a global study on the formation of vdW molecules through three-body collisions.

In the present work, we study the formation of vdW molecules X–RG through direct three-body recombination processes X + RG + RG → X–RG + RG. Here, RG indicates the rare gas atom and atom X has been chosen in a way to cover a broad range of chemical characteristics, i.e., from three different groups of the periodic table: alkali group (Li and Na), transition metals (Ti), and pnictogen group (As, P, and N). Our approach is based on a classical trajectory (CT) methodology in hyperspherical coordinates, which has been already applied to the three-body recombination of helium^{9,41} and to ion-neutral–neutral three-body recombination processes.^{12,13,42} Indeed, a direct three-body approach for the formation of vdW molecules has never been carried out via the CT method to the best of our knowledge. Performing these calculations, we notice a clear distinction between the formation rates at low-energy collisions and high-energy collisions, established by the dissociation energy of the X–RG potential. Moreover, our results show that the three-body recombination rate is of the same order of magnitude independent of the X atom, and hence, most of the vdW molecules X–RG should be observable in buffer gas cells.

This paper is organized as follows: In Sec. II, we summarize the main aspects of the classical trajectory method employed to study direct three-body recombination processes. In Sec. III, we precisely investigate the dependence of three-body recombination rates on the collision energy and temperature, utilizing six different systems. In Sec. IV, we discuss the applicability of the classical treatment at low temperatures. Finally, in Sec. V, we summarize our chief results and discuss their possible applications.

## II. CLASSICAL TRAJECTORY METHOD IN HYPERSPHERICAL COORDINATES

Consider a system of three particles with masses *m*_{i} (*i* = 1, 2, 3) at the respective positions $r\u2192i$, interacting with each other via the potential $V(r\u21921,r\u21922,r\u21923)$. Here, we neglect the three-body term of the potential, which is a good approximation for van der Waals molecules and clusters,^{43} and hence, *V* can be expressed as a summation of pair-wise potentials, i.e., $V(r\u21921,r\u21922,r\u21923)=U(r12)+U(r23)+U(r31)$, where $rij=|r\u2192j\u2212r\u2192i|$. The dynamics of these particles is governed by the Hamiltonian

with $p\u2192i$ being the momentum vector of the *i*th particle.

To solve Hamilton’s equations and find classical trajectories, it is more convenient to employ Jacobi coordinates.^{44,45} For a three-body problem, Jacobi vectors are related to $r\u2192i$ vectors as

where *M* = *m*_{1} + *m*_{2} + *m*_{3} is the total mass of the system and $\rho \u2192CM$ is the three-body center-of-mass vector. These vectors are illustrated in Fig. 1.

Due to the conservation of the total linear momentum (i.e., $\rho \u2192CM$ is a cyclic coordinate), the degrees of freedom of the center of mass can be neglected, and thus, the Hamiltonian (1) transforms to

Here, *μ*_{12} = *m*_{1}*m*_{2}/(*m*_{1} + *m*_{2}) and *μ*_{3,12} = *m*_{3}(*m*_{1} + *m*_{2})/*M*, and $P\u21921$ and $P\u21922$ indicate the conjugated momenta of $\rho \u21921$ and $\rho \u21922$, respectively. $V(\rho \u21921,\rho \u21922)$ is the potential expressed in terms of the Jacobi coordinates.

Noting that Hamilton’s equations of motion are invariant under the canonical transformation (2), it is possible to predict the evolution of the trajectories in terms of Jacobi coordinates from Hamiltonian (3) via

and transform the solutions back to Cartesian coordinates. As an example, Fig. 2 shows the classical trajectories calculated for Li + He + He three-body collisions for different collision energies and the same impact parameter (*b* = 0). The panel (a) of this figure shows an elastic or non-reactive trajectory in which the three-body collision leads to three free particles flying away form each other. On the contrary, in panels (b) and (c), the three-body collision ends up forming a molecule that vibrates rapidly, i.e., a three-body recombination event Li + He + He → Li–He + He.

### A. Classical three-body recombination in hyperspherical coordinates

It is well-known that, classically, *n*-body collisions in a three-dimensional (3D) space can be mapped into a problem involving one particle with a definite momentum moving toward a scattering center in a *d*-dimensional space in which *d* = 3*n* − 3 is equal to the independent relative coordinates of the *n*-body system. Exploiting this point, we define the initial conditions and impact parameter associated with a three-body problem as single entities in a six-dimensional (6D) space. The 6D space is described in hyperspherical coordinates, which consist of a hyperradius *R* and five hyperangles *α*_{j} (*j* = 1, 2, 3, 4, 5), where 0 ≤ *α*_{1} < 2*π* and 0 ≤ *α*_{j>1} ≤ *π*.^{46–52}

The position and momentum vectors in the hyperspherical coordinates can be constructed from Jacobi vectors and their conjugated momenta as

and

respectively, where $\mu =m1m2m3/M$ is the three-body reduced mass (for further details, see Refs. 41 and 53). Consequently, the Hamiltonian *H* in these coordinates reads as

In the 3D space, the collision cross section *σ* is defined as the area drawn in a plane perpendicular to the initial momentum containing the scattering center that the relative motion of the particles (known as trajectory) should cross in order for a collision to take place.^{18} This concept can be extended to a 6D space by visualizing it in a five-dimensional hyperplane (embedded in a 6D space) instead of a plane.^{41,53} Using the same analogy, we can define the impact parameter vector $b\u2192$ as the projection of the position vector in a 5D hyperplane perpendicular to the initial 6D momentum vector $P\u21920$ (i.e., $b\u2192\u22c5P\u21920=0$). Therefore, the cross section associated with the three-body recombination process, after averaging over different orientations of $P\u21920$, is obtained as follows:^{41}

where $d\Omega b=sin3(\alpha 4b)sin2(\alpha 3b)sin(\alpha 2b)d\alpha 4bd\alpha 3bd\alpha 2bd\alpha 1b$ is the solid angle element associated with vector $b\u2192$, and we made use of the relation $P0=2\mu Ec$. The function $P$ is the so-called opacity function, i.e., the probability that a trajectory with particular initial conditions leads to a recombination event. Note that the factor Ω_{b} = 8*π*^{2}/3 is the solid hyperangle associated with $b\u2192$, and $P(Ec,b)=0$ for *b* > *b*_{max}. In other words, *b*_{max} represents the largest impact parameter for which three-body recombination occurs. Finally, the energy-dependent three-body recombination rate is obtained as

### B. Computational details

The angular dependence of the opacity function $P(b\u2192,P\u21920)$, which depends on both the direction and magnitude of the impact parameter and initial momentum vectors, has been averaged out by means of the Monte Carlo method.^{41,54} Without loss of generality, we choose the *z* axis in 3D space to be parallel to the Jacobi momentum vector $P\u21922$. The initial hyperangles determining the orientation of vectors $P\u21920$ and $b\u2192$ in the 6D space are sampled randomly from probability distribution functions associated with the appropriate angular elements in hyperspherical coordinates (see Ref. 41).

In the next step, the opacity function $P(Ec,b)$ for a given collision energy, $Ec=P02/(2\mu )$, and the magnitude of impact parameter, *b*, are obtained by dividing the number of classical trajectories that lead to recombination events, *n*_{r}, by the total number of trajectories simulated *n*_{t}.^{41} Thus,

where the second term in Eq. (10) is the statistical error owing the inherent stochastic nature of the Monte Carlo technique.

To solve Hamilton’s equations, we made use of the explicit Runge–Kutta (4,5) method, the Dormand–Prince pair.^{55} The acceptable error for each time step has been determined by absolute and relative tolerances equal to 10^{−15} and 10^{−13}, respectively. The total energy is conserved during collisions to at least four significant digits and the magnitude of the total angular momentum vector, $J=|\rho \u21921\xd7P\u21921+\rho \u21922\xd7P\u21922|$, is conserved to at least six significant digits. The initial magnitude of hyperradius, $|\rho \u21920|$, is generated randomly from the interval [*R*_{0} − 25, *R*_{0} + 25] *a*_{0} centered around *R*_{0} = 550*a*_{0} (*a*_{0} ≈ 5.29 × 10^{−11} m is the Bohr radius). This value fulfills the condition for three particles to be initially in an uniform rectilinear state of motion.

## III. RESULTS AND DISCUSSION

Throughout this section, we consider the formation of weakly bound He-containing vdW molecules in their electronic ground state, through the three-body recombination process X + ^{4}He + ^{4}He → X–^{4}He + ^{4}He, for six different X atoms from three different groups of the periodic table. We consider ^{7}Li and ^{23}Na from the alkali group, ^{48}Ti from the transition metals, and ^{75}As, ^{31}P, and ^{14}N from the pnictogen group. All these atoms, with the exception of Ti, show an *S* electronic ground state.

Figure 3 displays the two-body potentials *U*(*r*) that have been used in the calculations (Li–He and Na–He from Ref. 56, Ti–He from Refs. 37 and 61, As–He, P–He, and N–He from Ref. 58, and He–He from Ref. 60). Note that all the X–He complexes show a single electronic state correlated with the ground electronic state of the atom and the rare gas atom, which are described as Lennard-Jones (LJ) potentials with the form *U*(*r*) = *C*_{12}/*r*^{12} − *C*_{6}/*r*^{6}. However, since the electronic ground state of Ti presents an *F* symmetry, Ti–He shows four different electronic states correlated with the ground electronic state of Ti and He atoms. In this case, we have taken the LJ potential fitted to the spherically symmetric component of the potential given as^{61–63}

where *U*_{Σ}, *U*_{Π}, *U*_{Δ}, and *U*_{Φ} are the distinct molecular potentials correlated with Ti–He in the ground electronic state. Note that the corresponding well depths range from *D*_{e} ≈ 1.87 K = 1.30 cm^{−1} (for Na–He) to *D*_{e} ≈ 19.74 K = 13.72 cm^{−1} (for N–He).

In the following, we present the three-body recombination rates calculated from the CT method and explore their dependence on the collision energy and on the particular features of the underlying two-body potentials.

### A. Energy-dependent three-body recombination rate for X–He–He systems

The energy-dependent three-body recombination rates, *k*_{3}(*E*_{c}), for the six considered cases are illustrated in Fig. 4. It is quite remarkable that, despite the drastic differences in the properties of X atoms and parameters of X–He interaction potentials, the recombination rates are of the same order of magnitude. Moreover, it is noticed that the energy-dependent three-body recombination rate shows the same trend as a function of the collision energy, independent of the X atom under consideration. In particular, we identify two power-law behaviors (linear in the log–log scale) connected at the dissociation energy, *D*_{e}, represented by the black dashed line in each of the panels of Fig. 4. Indeed, *D*_{e} acts as the threshold energy for two distinct regimes: the low-energy regime, where *E*_{c} < *D*_{e}, and the high-energy regime, where *E*_{c} > *D*_{e}.

The data displayed in Fig. 4 show that even though in both regimes, the dependence of *k*_{3} on *E*_{c} follows a power law, the energy dependence for the high-energy domain is much steeper than that for the low-energy one. In our view, this behavior is related to the interplay between the role of the long-range tail of the X–He potential and its short-range region in the formation of vdW molecules at different energies. In other words, the formation of vdW molecules at low energies is mainly a consequence of the X–He interaction potential’s long-range tail, but this is not the case for high-energy collisions.

To check the validity of this statement, we have computed the energy-dependent three-body recombination rate at different collision energies by varying the short-range part of the interaction potential while keeping *C*_{6} constant, and the results are shown in Fig. 5. In this figure, we observe that only when *E*_{c} > *D*_{e}, the three-body recombination rate shows a variation from its nearly constant value at *E*_{c} < *D*_{e}. Therefore, the precise details of short-range X–He interaction (properties of potential well) only matter at high collision energies, where the three-body recombination rate starts to show a steep behavior.

The power-law behavior at low-energy collisions is expected in the virtue of the classical nature of the collisions,^{41,53} as explained in the following.

#### 1. Low-energy regime

Based on our results for the energy-dependent three-body recombination rates of X–He formation, we conclude that the energy dependence of the recombination rate in the low-energy regime depends chiefly on the dominant long-range 1/*r*^{6} interaction. To study this dependence, we apply a classical capture model following the pioneering ideas of Langevin for ion–neutral reactions.^{64} In this framework, every trajectory with an impact parameter below some threshold value $b\u0303$ leads with unit probability to a reaction event. $b\u0303$ is given by the largest partial wave for which the height of the centrifugal barrier is equal to the collision energy. For neutral–neutral interactions, the effective long-range potential reads as (in atomic units)

The second term in Eq. (12) is the centrifugal barrier, with *μ*_{0} being the two-body reduced mass and *ℓ* being the angular momentum quantum number or partial wave. The potential *U*_{eff}(*r*) shows a maximum at

Classically, a reaction occurs if and only if *E*_{c} ≥ *U*_{eff}(*r*_{0}). Now, to find the critical impact parameter $b\u0303$, which is assigned to *E*_{c} = *U*_{eff}(*r*_{0}), we may use the relation between the angular momentum quantum number *ℓ*, the collision energy, and the impact parameter,^{18,53} i.e.,

Substituting *ℓ*(*ℓ* + 1) obtained from *E*_{c} = *U*_{eff}(*r*_{0}) into Eq. (14) yields

Applying this model to both X–RG and RG–RG interactions and keeping in mind that the 6D impact parameter *b* is a combination of the 3D impact parameters associated with the Jacobi coordinates $\rho 1\u2192$ and $\rho 2\u2192$, we expect the same power law for *b*_{max} (introduced in Sec. II A), i.e.,

Therefore, in virtue of Eq. (16) and the Langevin assumption that $P(Ec,b>bmax)=0$ and $P(Ec,b\u2264bmax)=1$, from Eq. (8), we obtain the low-energy power law for the three-body recombination cross section as $\sigma rec(Ec)\u221dEc\u22125/6$ and, hence, the three-body recombination rate as

which, as expected, is consistent with the classical threshold law that has been found for low-energy collisions in Ref. 41.

Let us now examine our findings via an example, namely, As + He + He interaction. Figure 6 shows the opacity function $P(Ec,b)$ for the formation of As–He due to a three-body recombination in terms of collision energy *E*_{c} and 6D impact parameter *b*. The white dashed line represents the collision energy equal to the dissociation energy, *E*_{c} = *D*_{e}, and the white curve indicates the $b\u221dEc\u22121/6$. As expected, the opacity function has its maximum at *b* = 0 and *E*_{c} = 10^{−3} K (the lowest illustrated energy). By increasing the impact parameter, the opacity function along each line at constant collision energy gradually decreases and eventually vanishes at *b* = *b*_{max}.

The white curve ($b\u221dEc\u22121/6$) in Fig. 6 reasonably resembles the loci of *b*_{max} in the low-energy regime. However, in higher energies, this loci deviates from the white curve, and for *E*_{c} ≳ *D*_{e} (top left corner) does not obey the same power law any more. Consequently, based on Eq. (17), we can now explain the observed trend of the recombination rates displayed in panel (d) of Fig. 4. While the energy dependence of *k*_{3} on the collision energies below 10^{−2} K can be conveniently explained by the adopted classical capture model, this model cannot provide the correct power law for the high energies above the dissociation energy *D*_{e} ≈ 16 K, where *b*_{max} varies much faster than $Ec\u22121/6$. In the intermediate regime connecting these two limits, the dependence of *k*_{3} on the collision energies gradually deviates from the initial relation given by Eq. (17) and is closer to $k3(Ec)\u221dEc\u22121/2$.

In addition, it is worth mentioning that including the three-body interaction term of the X–He–He potential energy surface will lead to a deviation from the derived power-law behavior for *b*_{max} and *k*_{3}(*E*_{c}).

#### 2. High-energy regime

To explore the effect of short-range details of potential on the formation of vdW molecules at higher collision energies, without loss of generality, we focus on the energy-dependent three-body recombination rate for Ti + He + He. The CT calculations are performed by means of three different potentials for the Ti–He interaction, namely, *ab initio* potential from Refs. 37, 61, and 62, Buckingham potential^{27} with the form *U*(*r*) = *C*_{1} exp(−*C*_{2}*r*) − *C*_{6}/*r*^{6}, and the LJ potential used in previous calculations. The results are displayed in Fig. 7. The *C*_{6} dispersion coefficient in the LJ and Buckingham potentials is derived from *ab initio* calculations.^{37,61}

As expected from our previous discussion, despite the non-negligible differences in the shape of the short-range interaction potential [compare red (LJ), blue (*ab initio*), and green (Buckingham) curves in the inset of Fig. 7], the three-body recombination rates *k*_{3}(*E*_{c}) at collision energies below the dissociation energy (*E*_{c} < *D*_{e}) are the same, which confirms our observation that low-energy collisions are dominated by the long-range tail of the X–He potential.

In contrast, proceeding to the high-energy regime, we spot two distinct behaviors. First, the three-body recombination rate related to the shallower X–RG potential (smaller *D*_{e}) shows a slightly steeper energy dependence, and accordingly, the power law is different [see panel (a) in Fig. 7]. Second, the energy dependence of the three-body recombination rate does not depend on the equilibrium distance of the X–He molecule as long as the dissociation energy is the same [see panel (b) in Fig. 7]. In other words, the dissociation energy seems to be the most relevant short-range parameter of the two-body potential affecting the recombination rate. Therefore, the inclusion of non-additive interactions on the X–RG–RG potential energy surface may lead to a slightly different trend.

It is important to note that the long-range tail of the *ab initio* potential contains higher order terms proportional to 1/*r*^{8}, 1/*r*^{10}, … coming from the spherical multipole moment expansion of the involved electronic clouds. However, the three-body recombination rates obtained for both LJ and *ab initio* potentials are identical in the whole energy regime as well as for Buckingham and *ab initio* potentials in the low-energy regime. Therefore, our results strongly suggest that the effect of long-range interaction on the formation of vdW complexes is mainly through the 1/*r*^{6} term of the dispersion potential.

### B. Temperature-dependent three-body recombination rate for X–He–He systems

In the final part of our discussion, we focus on the production of vdW molecules in buffer gas cells. In particular, we investigate the thermal averaged three-body recombination rate as a mechanism for the formation of X–He vdW molecules at temperatures 4 K ≤ *T* ≤ 20 K. The thermal average of the three-body recombination rate is obtained via integrating the energy-dependent three-body recombination rate Eq. (9) over the appropriate three-body Maxwell–Boltzmann distribution of collision energies, yielding

where *k*_{B} is the Boltzmann constant. The results obtained by performing the thermal average (18) for six different X + He + He reactive collisions for 4 K ≤ *T* ≤ 20 K are shown in Fig. 8. Comparing the data in panels (d)–(f) of Fig. 8, it is noticed that the trend and the magnitude of the three-body recombination rate is very similar for all three considered pnictogens, As, P, and N. Similarly, the temperature-dependent three-body recombination rate displays the same tendency for alkali metals Li and Na [panels (a) and (b)].

Finally, we notice that the three-body recombination rate *k*_{3}(*T*) for all the X atoms considered shows nearly the same order of magnitude, except for Na. This is due to the rapid decrease in the Na–He recombination rate *k*_{3}(*E*_{c}) at relatively low collision energies in comparison to the rest of X–He complexes [compare panel (a) with other panels of Fig. 4]. However, it is still fascinating that vdW molecules containing X atoms from totally different groups of periodic table show a similar three-body recombination rate within the range of typical temperatures in buffer gas cells. Indeed, in virtue of the observation of Li–He,^{35} Ag–He,^{33} and Ti–He,^{37} it should be possible to observe any X–He molecule in a buffer gas cell under the proper conditions.

## IV. RELIABILITY OF CLASSICAL TRAJECTORY CALCULATIONS AT LOW TEMPERATURES

The reliability of classical trajectory calculations for scattering observables depends on the collision energy at which the system is studied. In particular, the lower the collision energies, the more the quantum mechanical effects become prominent. In general, the importance of quantum mechanical effects on a system is related to the number of partial waves contributing to the scattering observables. Classically, the largest number of the allowed partial wave is calculated by setting the centrifugal barrier equal to the collision energy *E*_{c}.^{53} For the systems under consideration in this work, the two-body X–RG long-range interaction (i.e., $\u2212C6/rX\u2212RG6$) dominates the potential interaction. Thus, we have

with *μ*_{0} being the two-body reduced mass.

As an example, *ℓ*_{max} values for the X–He pairs are listed in Table I. Note that in the case of three-body collisions, the total angular momentum will be affected by the He–He interaction, and the expected number of partial waves may increase compared with the one obtained from the X–RG long-range interaction potential. Moreover, based on Eq. (19), at a given collision energy, heavier RG atoms show larger *ℓ*_{max}. This can be understood, considering that the heavier the system is, the closer to the classical realm it is. In this case, it is also related to the fact that heavier RG atoms show a larger static polarizability and hence a larger *C*_{6} (in general).

. | E_{c} (K)
. | ||
---|---|---|---|

X–RG . | 100 . | 10 . | 1 . |

Li–He^{a} | 15 | 7 | 3 |

–Ne^{b} | 24 | 11 | 5 |

–Ar^{c} | 32 | 15 | 7 |

–Kr^{c} | 36 | 16 | 7 |

–Xe^{c} | 39 | 18 | 8 |

Na–He^{a} | 17 | 8 | 3 |

–Ne^{b} | 35 | 16 | 7 |

–Ar^{c} | 51 | 24 | 11 |

–Kr^{c} | 61 | 28 | 13 |

–Xe^{c} | 68 | 31 | 14 |

N–He^{d} | 13 | 6 | 2 |

–Ar^{e} | 33 | 15 | 7 |

–Kr^{e} | 39 | 18 | 8 |

P–He^{d} | 16 | 7 | 3 |

As–He^{d} | 17 | 8 | 3 |

Ti–He^{f} | 18 | 8 | 4 |

. | E_{c} (K)
. | ||
---|---|---|---|

X–RG . | 100 . | 10 . | 1 . |

Li–He^{a} | 15 | 7 | 3 |

–Ne^{b} | 24 | 11 | 5 |

–Ar^{c} | 32 | 15 | 7 |

–Kr^{c} | 36 | 16 | 7 |

–Xe^{c} | 39 | 18 | 8 |

Na–He^{a} | 17 | 8 | 3 |

–Ne^{b} | 35 | 16 | 7 |

–Ar^{c} | 51 | 24 | 11 |

–Kr^{c} | 61 | 28 | 13 |

–Xe^{c} | 68 | 31 | 14 |

N–He^{d} | 13 | 6 | 2 |

–Ar^{e} | 33 | 15 | 7 |

–Kr^{e} | 39 | 18 | 8 |

P–He^{d} | 16 | 7 | 3 |

As–He^{d} | 17 | 8 | 3 |

Ti–He^{f} | 18 | 8 | 4 |

Therefore, due to the relatively large number of partial waves, we believe that the CT calculation presented in Sec. II is a reasonable approach to study the formation of vdW molecules even at energies near 4 K. For a more detailed comparison between the quantum and classical results obtained by the hyperspherical CT method in three-body collisions, see Ref. 41.

## V. CONCLUSIONS AND PROSPECTS

In summary, we have shown, via a classical trajectory method introduced in Ref. 41, that van der Waals molecules can be formed through direct three-body recombination. In particular, we have investigated the energy dependence and temperature dependence of the three-body recombination rates for six vdW complexes containing atoms with totally different chemical properties. As a result, we found that the X–RG molecule’s dissociation energy is the determinant parameter to differentiate between the low and high-energy regimes. At low energies, the formation rate of vdW molecules is relatively insensitive to the short-range interaction and is dominated by the long-range tail of the potential, i.e., −1/*r*^{6}. This regime is further explored by a classical capture model, à la Langevin. Conversely, at higher collision energies, the short-range part of the potential plays an important role.

However, the most exciting result is that the three-body recombination rate for the formation of vdW molecules is of the same order of magnitude, independent of the chemical properties of the atom colliding with the remaining two rare gas atoms. In other words, the formation rate of X–RG vdW molecules is almost independent of the chemical properties of the X atom. Indeed, we have shown that CT calculations are reliable for relevant temperatures in buffer gas cells (1 K ≲ *T* ≲ 10 K) based on the number of contributing partial waves to different scattering observables. Therefore, it should be possible to create and study X–RG vdW molecules in buffer gas cells, where some of them have been already observed and studied. Nevertheless, it is necessary to correctly identify the molecular dissociation processes to understand vdW molecules under equilibrium conditions, which can be considered as an extension of this work in the near future. Last but not least, in our view, a complete understanding of the formation process of vdW molecules will become a cornerstone of the physics of aggregation of vdW complexes.

## ACKNOWLEDGMENTS

We thank Dr. Timur V. Tscherbul for sharing with us the *ab initio* potential of Ti–He.

## DATA AVAILABILITY

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