On the formation of van der Waals complexes through three-body recombination

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 + X + RG $\rightarrow$ X-RG + X. 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 [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 $\lesssim T \lesssim 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 magnitude independently 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 labelled as a three-body recombination process or as 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][4][5][6][7][8][9], and the formation and trapping of cold and ultracold molecules [10][11][12][13][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 1 . The binding mechanism in vdW molecules relies on the compensation between the shortrange repulsive interaction (due to the overlap of closedshell 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 [22][23][24][25]; chemical reactions [24,[26][27][28][29][30]; and physical phenomena like superfluidity of 4 He nanodroplets [31,32]. 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 [33][34][35][36].
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 [23,26,27,[36][37][38][39][40][41][42][43]. Recently, thanks to the development of buffer gas sources [44], it has been possible to investigate the formation of vdW molecules through three-body recombination processes [36, 39-41, 43, 45, 46]. 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 (where RG is the rare gas atom) through direct three-body recombination processes X + X + RG → X-RG + X. Our approach is based on a classical trajectory (CT) methodology in hyperspherical coordinates, which has been already applied to three-body recombination of helium [9,47], and to ion-neutral-neutral threebody recombination processes [12,13,48]. Indeed, a direct three-body approach for forming vdW molecules has never been carried out via CT method to the best of our knowledge. Performing these calculations, we notice a clear distinction between the formation rates at lowenergy 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 independently 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 summarise 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 i , interacting with each other via the potential V ( r 1 , r 2 , r 3 ). Here, we neglect the three-body term of the potential, hence V can be expressed as a summation of pair-wise potentials, i.e., V ( r 1 , r 2 , r 3 ) = U (r 12 ) + U (r 23 ) + U (r 31 ), where r ij = | r j − r i |. The dynamics of these particles is governed by the Hamiltonian with p i 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 [49,50]. For a three-body problem, Jacobi vectors are related to r i vectors as where M = m 1 + m 2 + m 3 is the total mass of the system and ρ CM is the three-body center-of-mass vector. These vectors are illustrated in FIG. 1.
Due to the conservation of total linear momentum (i.e., ρ CM is a cyclic coordinate), the degrees of freedom of the center of mass can be neglected, 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 1 and P 2 indicate the conjugated momenta of ρ 1 and ρ 2 , respectively. V ( ρ 1 , ρ 2 ) is the potential expressed in terms of the Jacobi coordinates.
Noting that the 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 threebody 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 towards a scattering center in a d-dimensional space in which d = 3n−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 ≤ π [51][52][53][54][55][56][57].
Position and momentum vectors in the hyperspherical coordinates can be constructed from Jacobi vectors and their conjugated momenta as and respectively, where µ = m 1 m 2 m 3 /M is the three-body reduced mass (for further details see Refs. [47 and 58]). Consequently, the Hamiltonian H in this 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 to a collision to take place [24]. 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 [47,58]. Using the same analogy, we can define the impact parameter vector b as projection of the position vector in a 5D hyperplane perpendicular to the initial 6D momentum vector P 0 (i.e., b. P 0 = 0). Therefore, the cross section associated with the threebody recombination process, after averaging over different orientations of P 0 , is obtained as follows [47]: where 1 is the solid angle element associated with vector b, and we made use of the relation P 0 = √ 2µE c . 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, and P(E c , 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, P 0 ), which depends on both direction and magnitude of impact parameter and initial momentum vectors, has been averaged out by means of Monte Carlo method [47,59]. Without loss of generality, we choose the z axis in 3D space to be parallel to the Jacobi momentum vector P 2 . The initial hyperangles determining the orientation of vectors P 0 and b in the 6D space are sampled randomly from probability distribution functions associated with the appropriate angular elements in hyperspherical coordinates (see Ref. [47]).
In the next step, the opacity function P(E c , b) for a given collision energy, E c = P 2 0 /(2µ) and magnitude of impact parameter, b, is obtained by dividing the number of classical trajectories that lead to recombination events, n r , by the total number of trajectories simulated n t [47]. Thus, where the second term in Eq. (10) is the statistical error owing the inherent stochastic nature of the Monte Carlo technique.
To solve the Hamilton's equations we made use of the explicit Runge-Kutta (4,5) method, the Dormand-Prince pair [60]. 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 = | ρ 1 × P 1 + ρ 2 × P 2 |, is conserved to at least six significant digits. The initial magnitude of hyperradius, | ρ 0 |, 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 fulfils the condition for three particles to be initially in an uniform rectilinear state of motion. The black dotted curve indicates the He-He interaction based on parameters given in Ref. [65].

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.
FIG. 3 displays the two-body potentials U (r) that have been used in the calculations (Li-He and Na-He from Ref. [61]; Ti-He from Refs. [43 and 66]; As-He, P-He and N-He from Ref. [63]; and He-He from Ref. [65]). 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 spherically symmetric component of the potential given as [66][67][68] where U Σ , U Π , U ∆ , and U Φ are the distinct molecular potentials correlated with Ti-He in the ground electronic state.
In what follows, 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, independently 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 as 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 shows that even though in both regimes, the dependence of k 3 on E c follows a powerlaw, the energy-dependence for the high-energy domain is much steeper than 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 shortrange 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 shortrange 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 [47,58], as it is explained below.

Low-energy regime
Based on our results for the energy-dependent threebody recombination rates of X-He formation, we conclude that the energy-dependence of 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 [69]. In this framework, every trajectory with an impact parameter below some threshold valueb leads with unit prob- 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 the angular momentum quantum number or partial wave. The potential U eff (r) shows a maximum at r 0 = 6µ 0 C 6 ( + 1) Classically, a reaction occurs if and only if E c ≥ U eff (r 0 ). Now, to find the critical impact parameterb which is assigned to E c = U eff (r 0 ), we may use the relation between angular momentum quantum number , the collision energy, and the impact parameter [24,58], i.e., Substituting ( + 1) obtained from E c = U eff (r 0 ) into Eq. (14) yieldsb = 2 3 Applying this model to both X-RG and RG-RG interactions and keeping in mind that 6D impact parameter b is a combination of the 3D impact parameters associated with the Jacobi coordinates ρ 1 and ρ 2 , we expect the same power-law for b max (introduced in Sec. Therefore, in virtue of Eq. (16) and the Langevin assumption that P(E c , b > b max ) = 0 and P(E c , b ≤ b max ) = 1, from Eq. (8) we obtain the low-energy power-law for the three-body recombination cross section as σ rec (E c ) ∝ E −5/6 c , and hence, the three-body recombination rate as which, as expected, it is consistent with the classical threshold law that has been found for low-energy collisions in Ref. [47].
Let us now examine our findings via an example, namely, As + He + He interaction. FIG. 6 shows the opacity function P(E c , 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 ∝ E −1/6 c . 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 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 collision energies D e (top left corner) does not obey the same power-low 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 can not provide the correct power-law for the high energies above dissociation energy D e ≈ 16 K, where b max varies much faster than E −1/6 c . 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 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 ).

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 interac-  The C 6 dispersion coefficient in the LJ and Buckingham potentials is derived from ab initio calculations [43,66]. As expected from our previous discussion, despite the non-negligible differences in the shape of the short-range interaction potential in the short-range region [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 upper panel in FIG. 7). Second, the energy-dependence of the three-body recombination rate do not depend on the equilibrium distance of the X-He molecule as long as the dissociation energy is the same. In other words, the dissociation energy seems to be the most relevant short-range parameter of the two-body potential affecting the recombination rate.
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 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 ≤ T ≤ 20 K. The thermal average of the three-body recombination rate is obtained via integrating the energy-dependent threebody recombination rate Eq. (9) over the appropriate three-body Maxwell-Boltzmann distribution of collision energies, yielding  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 alkalimetals 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 a vdW coefficient C 6 is taken from Ref. [61]. b vdW coefficient C 6 is taken from Ref. [62]. c vdW coefficient C 6 is taken from Ref. [34]. d vdW coefficient C 6 is taken from Ref. [63]. e vdW coefficient C 6 is taken from Ref. [64]. f vdW coefficient C 6 is taken from Refs. [43 and 66].
the same order of magnitude, except for Na. This is due to the rapid decrease of 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 [41], Ag-He [39], and Ti-He [43], 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 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 [58]. For the systems under consideration in this work, the two-body X-RG long-range interaction (i.e., −C 6 /r 6 X−RG ) 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 larger static polarizability and hence a larger C 6 (in general).
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 hyperspherical CT method in threebody collisions see Ref. [47].

V. CONCLUSIONS AND PROSPECTS
In summary, we have shown, via a classical trajectory method introduced in Ref. [47], 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 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 threebody recombination rate for the formation of vdW molecules is of the same order of magnitude, independently 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 X atom. Indeed, we have shown that CT calculations are reliable for relevant temperatures in buffer gas cells (1K 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 in 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 vdW complexes' aggregation.