Quantum-Mechanical Force Balance Between Multipolar Dispersion and Pauli Repulsion in Atomic van der Waals Dimers

The structure and stability of atomic and molecular systems with van der Waals (vdW) bonding are often determined by the interplay between attractive dispersion interactions and repulsive interactions caused by electron confinement. Arising due to different mechanisms -- electron correlation for dispersion and the Pauli exclusion principle for exchange-repulsion -- these interactions do not appear to have a straightforward connection. In this paper, we use a coarse-grained approach for evaluating the exchange energy for two coupled quantum Drude oscillators and investigate the mutual compensation of the attractive and repulsive forces at the equilibrium distance within the multipole expansion of the Coulomb potential. This compensation yields a compact formula relating the vdW radius of an atom to its multipole polarizabilities, $R_{\rm vdW} = A_l^{\,}\, \alpha_l^{{2}/{7(l+1)}}$, where $l$ is the multipole rank and $A_l$ is a conversion factor. Such a relation is compelling because it connects an electronic property of an isolated atom (atomic polarizability) with an equilibrium distance in a dimer composed of two closed-shell atoms. We assess the accuracy of the revealed formula for noble-gas, alkaline-earth, and alkali atoms and show that the $A_l$ can be assumed to be universal constants. Besides a seamless definition of vdW radii, the proposed relation can also be used for the efficient determination of atomic multipole polarizabilities solely based on the corresponding dipole polarizability and the vdW radius. Finally, our work provides a basis for the construction of efficient and minimally-empirical interatomic potentials by combining multipolar interatomic exchange and dispersion forces on an equal footing.


I. INTRODUCTION
Noncovalent interatomic and intermolecular interactions represent one of the key factors that determine the physicochemical properties of molecules and materials across chemistry, biology, and materials science [1][2][3][4]. Noncovalent interactions are traditionally classified in a perturbative formalism, from which electrostatics, induction, Pauli (exchange) repulsion, and van der Waals (vdW) dispersion arise as the leading contributions from the first two orders of perturbation theory. From the perspective of computational modeling, the individual terms are usually treated with different effective approaches. Especially the methods used to describe Pauli repulsion and vdW dispersion typically rely on fundamentally different physical models. The vdW dispersion represents a major part of long-range electron correlation forces arising from Coulomb-coupled instantaneous quantum fluctuations of the electronic charge distribution [5][6][7][8][9][10]. Common (semi)local approximations to density-functional theory (DFT), representing one of the main workhorse methods in atomistic modeling, neglect long-range correlation forces and thus do not account for vdW interactions. In recent years, an intense effort has been devoted to develop robust approaches to address this challenge [11][12][13][14][15][16]. Although a unified vdW functional valid for all kinds of systems is still under construction [17], significant progress has been achieved to include dispersion interactions in the form of nonlocal (vdW) density functionals [18][19][20][21]. Furthermore, coarse-grained vdW models have shown great success in describing dispersion interactions at lower computational costs [22][23][24][25][26]. Among them, the quantum Drude oscillator (QDO) model [27][28][29][30][31][32][33] has been firmly established as an efficient and accurate approach for modeling and understanding vdW interactions [32][33][34][35][36][37][38]. Within this approach, each QDO models an atom or a molecule, representing the effective, localized response and polarization fluctuation of its valence electrons. The success of the coupled-oscillator model is exemplified by its excellent description of the electronic response properties of atoms and molecules. In a continuous formalism, with one oscillator at every point in space, coupled oscillators can describe any response allowed by quantum field theory and thus model the response of arbitrary molecules or materials [39,40]. In the common practical coarse-grained formalism, with each oscillator representing one atom, the QDO framework reproduces the leading-order behavior of the electronic polarizability of atoms [41], providing an accurate and reliable description of polarization effects in molecules and materials [42,43]. Moreover, the QDO model allows one to describe excess electrons in matter [28] and to reproduce dispersion-polarized electron densities [38] as well as Coulomb interactions between dipolar quantum fluctuations [44].
Extending the applicability of the QDO framework towards a more complete and systematic description of noncovalent interactions necessitates the incorporation of the exchangeinduced repulsion [45]. Recently, we made a first step in this direction by evaluating the exchange energy between two QDOs within the dipole approximation of the Coulomb potential [46]. Here, we take the next natural step by constructing a common coarse-grained approach for the multipolar dispersion and exchange interactions in vdW-bonded atomic or molecular dimers.
It is important to embed our developments of coarsegrained models into the broader context given by the theory of intermolecular interactions for systems composed of nuclei and electrons [1], which states that the equilibrium geometry of two vdW-bonded atoms or molecules is governed by an interplay of several interactions. The generalized Heitler-London (GHL) theory [47] offers one of the most compact schemes for the interatomic energy decomposition. In the GHL approach, isotropic closed-shell atoms only experience mutual exchange-repulsion and dispersion forces. Another very successful scheme to describe intermolecular interactions and analyze their complex interplay is based on the symmetry-adapted perturbation theory (SAPT) decomposition [48][49][50]. The higher-level SAPT methods, while being computationally expensive, approach a "gold standard" accuracy [51] comparable to the coupled-cluster method with single, double, and perturbative triple excitations [CCSD(T)] for small molecules. Within second-order SAPT, which is the most practical approach, one obtains six contributions [1,45]: (i) electrostatics, E (1) elst ; (ii) exchange, E (1) ex ; (iii) induction, E (2) ind ; (iv) exchange induction, E (2) ex-ind ; (v) dispersion, E (2) disp ; and (vi) exchange dispersion, E (2) ex-disp . Here, the superscripts "(1)" and "(2)" denote the order of the perturbation theory required to derive the corresponding term. In the case of neutral and isotropic fragments, the two induction contributions, E (2) ind and E (2) ex−ind , practically compensate each other [52]. Then, the problem reduces to four remaining terms, which still yield significant contributions to the interaction energy for noble-gas dimers [53]. On the other hand, the Tang-Toennies (TT) model [54], relying just on the exchange-repulsion and dispersion-attraction interactions, is known to reproduce the binding energy curves of closed-shell dimers with high accuracy and efficiency [55]. Recently, an extension (TT2) of this model was proposed [56] to accurately describe noble-gas dimers also at relatively short internuclear distances. Based on the concepts of the GHL theory for interatomic interactions [47], the TT model can be considered as one of the most compact yet accurate models for closed-shell vdW dimers. According to the discussion in Ref. [47], the simplicity of the TT potential arises due to the used analytical asymptotic form of the exchange energy obtained by the surface integral method [57,58]. Since this method is known to deliver the same asymptotic result [59] as the approach based on the multipole expansion of the perturbing potential [60,61], the latter can be used as an alternative way to construct compact TT-like potentials. This idea is supported by our recent study [46], which established a quantum-mechanical scaling law, α 1 ∝ R 7 vdW , between the atomic dipole polarizability and the vdW radius from the force balance between exchange repulsion and dispersion attraction at the equilibrium distance. The corresponding analysis in Ref. [46] was based on the consideration of these two forces stemming from the dipolar term in the multipole expansion [2,3] of the interatomic Coulomb potential. Subsequently, we have derived [62] the proportionality coefficient, which finally led to the relation α 1 = (4π 0 /a 4 0 )α 4/3 f R 7 vdW , as expressed in terms of the vacuum permittivity 0 , the Bohr radius a 0 , and the fine-structure constant α f . Such a relation is not trivial because it connects an electronic polarizability of an atom with an equilibrium distance in a dimer composed of two closed-shell atoms.
In this paper, we build on our previous study by going beyond the dipole approximation and considering further terms in the multipole expansion of the interatomic Coulomb potential. This is performed for both exchange and dispersion interactions between closed-shell systems described within the QDO model. To this end, we investigate the balance between the two types of forces, which yield the dominant contributions in vdW-bonded systems. For atomic dimers at the vdW equilibrium distance, this allows us to study a term-by-term compensation of the attractive (dispersion) and repulsive (exchange) forces for each contribution in the multipole expansion of the full Coulomb interaction between the QDOs. This mutual compensation yields a relation between atomic multipole polarizabilities and the vdW radius as first empirically obtained in Ref. [46]. The presented relation enables a practical and seamless determination of vdW radii as an effective atomic length scale from atomic polarizabilities. From the opposite perspective, the generalized relation also allows one to obtain atomic polarizabilities across the periodic table and at arbitrary multipole rank based on (pretabulated) vdW radii and without the need to resort to the otherwise challenging direct computation of electronic response properties. Altogether, our results deliver deeper insights into the connection between Pauli repulsion and dispersion attraction-two forces which appear at different orders of SAPT. The existence of a quantum-mechanical relation between the two main contributions to the vdW interaction energy at the equilibrium distance reveals a strong connection between exchange and correlation effects and should have implications for achieving an improved understanding of the stability of vdW-bonded matter.

II. METHOD: QUANTUM DRUDE OSCILLATORS
Let us consider two vdW-bonded atoms, A and B, separated by a distance R, and describe them within the QDO model, as illustrated in Fig. 1. Each of the two QDOs representing atoms has three effective parameters-mass μ, charge q, and characteristic frequency ω-which are parametrized to reproduce three atomic observables {α 1 , C 6 , C 8 } [33]: where the Drude (quasi)particle and the related nucleus have charges (−q) and q, respectively. The conditions of Eq. (1) use the dipole polarizability α 1 and the dominant dispersion coefficients C 6 (induced-dipole-induced-dipole interaction) and Schematic representation of two QDOs separated by the distance R = |R|. The black and white spheres represent the two Drude particles and fixed nuclei, respectively. The Coulomb interactions between the two QDOs (gray arrows) and the harmonic intra-QDO potentials (blue arrows) are explicitly highlighted with their connection to Eqs. (2) and (4). C 8 (induced-dipole-induced-quadrupole interactions) in order to parametrize this powerful model, able to efficiently reproduce long-range forces and electronic response properties of atoms and molecules.
The Hamiltonian of the interacting QDOs is given by The corresponding wave functions are given by The full Coulomb interaction between the two QDOs is where r 1 and r 2 are the coordinates of the Drude particles measured from the corresponding fixed nuclei. The Coulomb interaction can be written as a multipole expansion as [63] V = n=1,2,...
where l A and l B refer to the rank of the multipole moments of the two interacting QDOs. Here, we restrict our consideration to the first five terms in the multipole expansion of Eq. (5). The first term, V 1 ∝ R −3 , corresponds to the dipole approximation of the Coulomb potential, ), describing the dipoledipole (d-d) electrostatic interaction (l A = l B = 1). The higher terms arise from the dipole-quadrupole (d-q) interaction for The formulas for V n , with n = 2, 3, 4, and 5, are given in Appendix A. Within the next section, we consider the multipolar contributions to the dispersion and exchange interaction between two QDOs. The analytical formulas are derived in the most general form valid in any system of units, whereas we employ atomic units (a.u.; with 4π 0 =h = 1) to present our numerical results in Sec. III D.

A. Dispersion interaction
The multipole expansion has been the starting point for quantum-mechanical perturbation calculations of the vdW dispersion interactions of Coulomb-coupled Drude oscillators [15,28,64]. Owing to this approach, the vdW dispersion energies can be expressed in terms of the atomic multipole polarizabilities (with l = 1, 2, . . .) (6) by using the series expansion [33] where T AB l A , l B represents the multipole-multipole interaction tensor. We remark that T AB l A , l B above has been obtained using a spherical harmonic expansion of the Coulomb potential instead of the Cartesian multipolar potential described in Appendix A. Both expansions yield equivalent results [65,66]. In the Supplemental Material of Ref. [33], the following spherical components of this tensor were given: For our derivations here, we further introduce the higher-order coupling components via a generalized expression of the multipolar interaction tensor, This tensor is derived based on the approach of Ref. [67] used by some of us in Ref. [68] as well. Popelier et al. [67] employed the relation where the expression in the large parentheses is a Wigner's 3j symbol [69] and the irregular normalized spherical harmonics are 033181-3 The spherical harmonics are defined as [69] Y l,m (θ, φ) = 2l + 1 4π where P m l (cos θ ) are the associated Legendre polynomials. If we assume now that the distance between two atoms is along the z axis, R = (0, 0, R), then cos θ = 1. Due to P m l (1) = δ m,0 , one can then easily obtain Consequently, we have Now we use the following property of Wigner's 3j symbols [69]: which gives us Obviously, |T AB l A , l B | 2 = |T AB l B , l A | 2 . Therefore it is enough to derive the components with l A l B . This means that where the factorials have been rewritten in terms of the binomial coefficients ( n k ) = n! (n−k)! k! . Then, we obtain in addition to the results of Eq. (8). With the above expressions, the first few multipolar contributions to the dispersion energy between two oscillators become where k e = (4π 0 ) −1 is the Coulomb constant. Based on the above formulas, we can now rewrite the dispersion energy in its conventional expansion [1,33] where C AB (2n+4) are the dispersion coefficients and all the contributions to E AB,disp n , with n up to 5, are given by Eq. (19). Equation (20) arises from second-order perturbation theory with the multipole expansion of the Coulomb potential, as an interaction potential between spherically symmetric atoms. The leading term is the dipole-dipole (d-d) interaction, E AB,disp 1 ∝ R −6 , stemming from the dipolar potential, V dip ∝ R −3 . The higher-order terms in the multipole expansion of the Coulomb interaction yield the dispersion energies E AB,disp and E AB,disp 5 ∝ R −14 coming, respectively, from the instantaneous dipole-quadrupole (d-q); dipole-octupole (d-o) and quadrupole-quadrupole (q-q); dipole-hexadecapole (d-h) and quadrupole-octupole (q-o); and dipole-triakontadipole (d-t), quadrupole-hexadecapole (q-h), and octupole-octupole (o-o) interactions. For noncentrosymmetric molecules, Eq. (20) would have terms with odd powers in R starting with ∝ R −7 [1]. However, here we restrict our consideration to vdWbonded atoms assumed to possess closed valence-electron shells with a spherically symmetric charge density, for which the dispersion terms proportional to R −(2i+1) , with i ∈ N, vanish.

B. Exchange-repulsion interaction
The above derivation of the dispersion energy was performed for the general case of two QDOs with arbitrary parameters; however, the description of the exchange repulsion between two QDOs is more delicate. The exchange interaction should obviously be present for two different  QDOs, as caused by the Pauli repulsion between electrons constituting the two Drude particles. Nonetheless, in order to construct the exchange interaction, one needs to deal with indistinguishable particles, a concept that requires generalization for two Drude particles possessing different parameters. Our starting assumption is that the exchange energy should be proportional to the overlap integral S between the wave functions of two different QDOs, similar to the case of two identical Drude particles [46]. This assumption was recently employed [70] for a simplified generalization of the coarsegrained dipole-dipole exchange energy of a homonuclear dimer, E ex (d−d) ≈ k e q 2 S/2R, derived in Ref. [46]. The authors of Ref. [70] have simply replaced the overlap integral S of two identical QDOs by its counterpart obtained for different QDOs and shown that already such a simplified treatment improves their computational scheme for vdW dispersion interactions. However, due to the coarse-grained treatment of valence electrons within the QDO model, care needs to be taken for the most general definition of the exchange energy between QDOs. This is a subject of our ongoing studies. Here, we follow the approach of Ref. [46] and derive multipole contributions to the exchange energy of two identical QDOs.
Formally, we consider two indistinguishable Drude par- as bosons assuming that they represent closed valence shells with vanishing total spin. Therefore the total wave function of a dimer should be written as a permanent By employing the Heitler-London perturbation theory [71,72], the exchange energy for two identical vdW-bonded QDOs at their equilibrium distance becomes well approximated with its exact asymptotic result given by the exchange integral [46] The evaluation of Eq. (22) with the expansion of Eq. (5) results in multipole contributions to the exchange energy, where each of them is directly proportional to the overlap integral defined as For the dipole-dipole contribution, which reproduces the result of Ref. [46]. Now, we evaluate further contributions going beyond the dipole approximation. For the dipole-quadrupole interaction, described by the second term, V 2 ∝ R −4 , in the multipole expansion of the Coulomb potential, we derive Then, the next term, V 3 ∝ R −5 , has two contributions, V 3(d−o) and V 3(q−q) , related to the dipole-octupole and the quadrupole-quadrupole interaction [3,73], respectively. The corresponding exchange integrals are obtained as Further on, we have two contributions from V 4 ∝ R −6 , the dipole-hexadecapole (d-h) and the quadrupole-octupole (q-o) interaction. The related exchange integrals are Finally, for the dipole-triakontadipole (d-t), quadrupolehexadecapole (q-h), and octupole-octupole According to Eqs. (24) and (25) for all interatomic distances. This is in contrast to the dispersion contributions, where E disp 1(d−d) clearly dominates at large distances. Such a nonmonotonic behavior of the multipole contributions, as we obtain here for the exchange energy, was also found in Ref. [74] for the multipole expansion of the exchange-dispersion energy. Now, we will use the derived dominant multipole contributions to the dispersion and exchange energies, in order to study the balance of the corresponding forces at the equilibrium distance in homonuclear dimers.

C. Force balance between multipolar dispersion and exchange contributions
The equilibrium geometry of atomic or molecular systems is dictated by the condition that the net forces acting on each atom vanish. Therefore, for two atoms or molecules separated by a distance R, this condition is determined by F net (R eq ) = −∇ R E tot (R)| R=R eq = 0, where R eq represents the equilibrium distance and E tot is the total interaction energy. The structure, stability, and dynamics of vdW-bonded atomic dimers are governed by the interplay between the dispersion and exchange interactions [47]. This means that at R = R eq the two respective forces have to mutually compensate each other. In what follows, we consider such a compensation by going beyond the dipole approximation for the interaction, in order to obtain higher-order multipole contributions to the attractive and repulsive forces.
At the equilibrium distance, R eq = 2R vdW , in homonuclear dimers composed of two identical Drude particles (μ = μ A = μ B , ω = ω A = ω B , and q = q A = q B ), the exchange force can be well approximated by the expression F ex ≈ −∇ R J ex ,  according to Ref. [46] and our discussion above. In addition, at the internuclear distances comparable to or larger than the equilibrium one, R √h /μω [46]. Then, the corresponding multipole contributions to the exchange force are obtained as From Eq. (19) we calculate the multipole contributions to the dispersion force (for homonuclear dimers) At R = R eq , the attractive and repulsive forces should cancel each other. Within the dipole approximation, from the force balance, This formula not only expresses a relation between α 1 and R vdW = R eq /2 but also contains the QDO parameters μ and ω, which are not uniquely defined for atoms. To obtain a formula connecting atomic parameters R vdW and α 1 , we rewrite Eq. (31) as where σ QDO = √h /2μω is the spatial variance or spread of a QDO. Within the QDO model, σ QDO describes an effective atomic length, which corresponds to the Bohr radius in the case of the hydrogen atom [75]. According to Ref. [46], the ratio R vdW /σ QDO decreases with increasing σ QDO , and the factors σ 4 QDO and e −(R vdW /σ QDO ) 2 in Eq. (32) compensate each other. This compensation allows the QDO model to approximately capture the constant behavior of the ratio α 1 /R 7 vdW confirmed empirically for many atoms [46]. Therefore, within the QDO model, the relation between the vdW radius and the dipole polarizability can be expressed as where the proportionality coefficient, as a function of the product μω and the vdW radius, is given by As was discussed in Ref. [46], this coefficient can be also written in terms of the radial volume occupied by the ground-state charge density of the QDO, n 0 (r) ≡ | 0 (r)| 2 = ( μω πh ) 3 /2 e − μω h r 2 , and its value at the vdW radius, n 0 (R vdW ), as Taking into account that A 1 was found to be essentially a constant for 72 atoms in the periodic table [46], Eq. (36) suggests a relation between an atomic volume and the electron charge density at the vdW radius.
The results of Eqs. (34) and (36) are based on taking into account only the first term, V 1 ∝ R −3 , in the expansion of Eq. (5) for the Coulomb potential. However, it is well known that at least two further terms, V 2 ∝ R −4 and V 3 ∝ R −5 , are important to properly describe the binding curves of vdWbonded atomic dimers [55]. Therefore, here we consider an extension of Eq. (33) by including the higher-order multipole terms from the expansion of the Coulomb potential. To this end, we evaluate the individual multipole contributions to the dispersion and exchange forces given by Eqs. (29) and (30), respectively, at the equilibrium distance of vdW-bonded dimers. Based on the empirical findings of Ref. [46], the mutual compensation of multipolar dispersion and exchange forces can ultimately be represented by the general expression which extends Eq. (33) to the multipole polarizabilities, α l . One can also rewrite Eq. (37) in the following way: where each multipole polarizability is expressed in terms of the vdW radius. This allows one to obtain α l either from R vdW or α 1 , for an arbitrary l. Since first-principles calculations of higher-order polarizabilities are computationally demanding [76], our finding provides an alternative way to approximate multipole polarizabilities. Based on the expressions for the multipolar dispersion and exchange forces, we can now also explicitly calculate the proportionality coefficients A μω l ≡ A l (μω, R vdW ). As a particularly helpful example, we can consider the force balance condition for the dipole-multipole interaction, i.e., the F l[d−z(l )] terms of Eqs. (29) and (30) [85]. For noble-gas atoms and hydrogen, missing in Ref. [85], the reference vdW radii are taken from Refs. [24,86], respectively. The QDO parameters {q, μ, ω} are set according to Eq. (1), to reproduce α ref 1 as well as the homoatomic dispersion coefficients C 6 and C 8 . The three fitted quantities together with the reference quadrupole (α ref 2 ) and octupole (α ref 3 ) polarizabilities are taken from Refs. [46,80] for the noble gases (He, Ne, Ar, Kr, Xe), from Refs. [77,78] for the elements in group I (H, Li, Na, K, Rb, Cs), and from Ref. [79] for the elements in group II (Be, Mg, Ca, Sr, Ba). The QDO multipole polarizabilities are obtained from Eq. (6)  with the l-dependent rational constants D l = l (l+2)(2l+1) 2 l+5 (l+1) . Equation (39) generalizes Eq. (34) and shows an additional factor of R 3(l−1)/7(l+1) vdW arising for l > 1. It is worth noting that the derived A μω l within the QDO model formally still contain R vdW . The values for A μω l calculated from Eq. (39), however, remain almost constant for any choice of realistic parameters (cf. Table I) as was also observed for the corresponding empirical proportionality factors [46]. Moreover, in terms of the quantities related to l = 1, the above expression can be simplified even further The general formula given by Eq. (39) allows us to obtain the proportionality coefficients A μω l for every order in the dipole-multipole interactions, even without deriving further multipolar contributions to the dispersion and exchange forces.
The presented findings based on the dipole-multipole interaction can be generalized via the force balance at each order of the multipole expansion as we highlight for the quadrupole-quadrupole and octupole-octupole interactions in Appendix B. Alternatively, one can use the general expression for the QDO multipole polarizabilities given by Eq. (6), in order to derive A QDO l = R vdW /α 2/7(l+1) l,QDO by means of Eq. (37). A comparison between the two approaches to the proportionality coefficients A l is given in the following section.

D. Assessment of our formalism for atoms
In the previous sections, we have presented a coarsegrained approach to describe dispersion and exchange interactions between two closed-shell atoms within the QDO model. Here, we examine the applicability of the presented formulas and apply them to analyze the ratio between the vdW radius and multipole polarizabilities for atoms, thus demonstrating the validity of the scaling law of Eq. (37) obtained within the QDO model. Our analysis will be focused on hydrogen, noble gases from He to Xe, alkali atoms from Li to Cs, and alkaline-earth elements from Be to Ba. To this end, the atomic multipole polarizabilities, α l , are either taken from high-level ab initio calculations in the literature [77][78][79] While eventually it would be interesting and important to extend our analysis to a broader set of atoms and small molecules, we are not aware of a comprehensive set of accurate data for atomic and molecular multipole polarizabilities. Accurate ab initio reference calculations of α l in general require demanding computational approaches with sophisticated treatment of electron correlation effects and, especially with increasing order l, large and diffuse basis sets [76,81,82].

033181-7
As a result, calculating converged multipolar polarizabilities is difficult, a problem which is further enhanced by the numerical aspects associated with finite-field derivative techniques as used in such calculations. Experimental determination, on the other hand, is subject to origin and orientational dependencies as well as a strong influence of thermal effects [76,83,84], which can introduce considerable uncertainties-in particular, with increasing system size or multipole order.
To apply the derived formulas, a set of reference vdW radii, R ref vdW , is required. In the case of alkali as well as alkaline-earth elements, these radii are taken from the recent database of Batsanov [85]. For noble-gas atoms, missing in Ref. [85], we use the database of Bondi [86], which provides often used values of vdW radii for group 18 of the periodic table. In addition, for hydrogen, we use R ref vdW = 3.1 a.u. from Ref. [24], where it was theoretically estimated based on the atomic charge density. This value was shown to work well for the relation between the atomic dipole polarizability and vdW radius in Ref. [46].
Both the vdW radii of Batsanov and those of Bondi are extracted from experimental crystallographic structural data. However, it is important to mention that a straightforward definition of the vdW radius is only possible for noble-gas atoms as inert elements with closed valence shells. For other atoms, R vdW is evaluated by considering a variety of different molecular crystal structures and extracting neighboring atom-atom distances, where each atom belongs to a different closed-shell molecule. This definition is especially subtle for chemical elements with spin-polarized valence shells, such as alkali atoms, which can form bonds with different spin states. Therefore one has to keep in mind that existing vdW radii are just statistical quantities for most chemical elements.
First, we analyze the empirical proportionality constants  Table I, where as a general trend alkali-metal atoms show the largest deviation with an average relative deviation of 4.3% (compared with 2.4% for alkaline-earth metals and 0.5% for noble gases). The alkali atoms possess a relatively weakly bound and, therefore, highly polarizable single valence electron. This feature and the different possible spin states of alkali atoms in molecular solid-state systems arguably allow the vdW radii observed for alkali metals to  change widely, causing the largest deviation for the empirical constants A ref l from their average values, as also illustrated in Fig. 2.
Overall, the observed deviations in A ref l stay within 4% (=0.1 a.u.) for all considered species except for Cs and Rb, which show a slightly higher average relative deviation of 6.3 and 5.5%, respectively. Hence we suggest that the coefficients A ref l can be considered as universal constants for the studied atomic species. Taking the average values A l = A ref l for the noble-gas atoms, we can thus write the unified relations between the vdW radius and the dipole, quadrupole, and octupole polarizabilities, 27 a.u., which are equivalent to the empirical relations reported in Ref. [46]. The relations obtained above can be used for at least three different purposes. First, the vdW radius of atoms can now be calculated given any single multipolar atomic polarizability. This polarizability can correspond to a free atom or an atom in a molecule or material [24]. The vdW radius can then be used for a conceptual understanding of an atom in its environment or for practical calculations of the vdW energy [17,24]. Second, given the dipole polarizability, one can accurately determine multipole polarizabilities (at least up to octupole) from Eq. (42). In fact, this approach is substantially more reliable for atoms than using the QDO model for multipole polarizabilities. A third application would be the possibility to determine atomic multipole polarizabilities from calculated or measured atomic vdW radii. A potential downside of this application is that a small error in the vdW radius would result in a large error for multipole polarizabilities, according to Eq. (42).
The quantity R vdW (α l ) defined in Eq. (42) represents an effective vdW radius expressed in terms of the multipolar polarizability. To demonstrate the validity of this new definition for the vdW radius, Fig. 3(a) compares the results for R vdW (α l ) = A l α 2/7(l+1) l with l = {1, 2, 3} to the reference values, R ref vdW . There is an excellent correlation between R vdW (α 1 ) and its reference value for all considered elements, with a maximum relative error R.E. = |R vdW (α 1 ) − R ref vdW |/R ref vdW of 0.91% for the noble gases, 2.74% for alkaline-earth (group II) elements, and 5.90% for hydrogen and alkali metals (group I). The increasing errors when going from alkaline-earth to alkali metals are related to the increase in the statistical errors of R ref vdW stemming from the increasingly complicated evaluation of the vdW radii based on experimental crystal-structure data. Indeed, this evaluation becomes less accurate for elements with more pronounced metallic properties [85,87]. Comparing groups I and II of the periodic table, the statistical errors in R ref vdW of the alkaline-earth elements are smaller since they have a closed s-electron shell, which makes their behavior closer to that of the noble gases with completely closed valence shells.
Although the dipole polarizability α 1,ref is known with high accuracy for many chemical elements in the periodic table, the accurate determination of higher-order multipolar polarizabilities is more involved. Indeed, Fig. 3(a) shows an increase in the maximum R.E. for R vdW (α 2 ) (within 0.99% for noble gases, 1.48% for group I, and 5.14% for group II) and, subsequently, for R vdW (α 3 ) (within 1.39% for noble gases, 6.60% for group I, and 8.62% for group II). This analysis is validated in Fig. 2 as well, where we compare the average values A 1 = 2.54 a.u., A 2 = 2.45 a.u., and A 3 = 2.27 a.u. with the relative ratios R vdW /α 2/7(l+1) l , for l = {1, 2, 3}. It is also noteworthy that reference values for higher-order polarizabilities are rather limited in the literature, with the exception of hydrogen, the only element in the periodic table for which the multipole polarizability α H l is known analytically [88]. Therefore we employ the multipole polarizabilities obtained within the QDO model by means of Eq. These results are shown in Fig. 3(b), where a good agreement between R QDO vdW (α l ) and R ref vdW is observed for l = {2, 3, 4, 5}, in addition to the case of l = 1, for which we have R QDO vdW (α 1,QDO ) ≡ R vdW (α 1,ref ). We note that the QDO model is constructed, by definition, on the dispersion coefficients, and the QDO polarizabilities (with l > 1) are underestimated for the noble-gas atoms with respect to the reference data [33]. Consequently, the QDO proportionality coefficients A QDO 2 = 2.52 a.u. and A QDO 3 = 2.37 a.u. are overestimated with respect to the determined "universal" values A 2 = 2.45 a.u. and A 3 = 2.27 a.u., as also shown in Fig. 2. Therefore one can  expect that the higher-order QDO coefficients, A QDO 4 = 2.23 a.u. and A QDO 5 = 2.10 a.u., are also overestimated. To further assess the scaling law of Eq. (37), we compare the resulting empirical constants with the proportionality coefficients A μω l , obtained by means of Eq. (39). Table I Table I for the corresponding A μω l . Hence the simplifications in the coarse-grained description of valence electrons within the QDO model and its parametrization lead to a less accurate determination of the proportionality coefficients, based on Eq. (39), than their indirect evaluation based on the QDO multipole polarizabilities calculated by means of Eq. (6). However, even A QDO l obtained for noble gases still show noticeable deviations from a constant behavior. It is unclear yet what aspect of atoms makes A ref l behave as universal constants for different chemical elements. Nevertheless, we expect the clarification of this question to be crucial for the eventual improvement of the QDO model.
Finally, we calculate the ratio R vdW /α 2/7(l+1) l for the hydrogen atom, taking into account the known analytical expression of its multipole polarizabilities given in Ref. [88] as α H l = (4π 0 ) a 2l+1 0 (2l + 1)!(l + 2)/2 2l l, and compare it with the proportionality coefficients A QDO l , which are obtained from α l,QDO provided by Eq. (6). As shown in Fig. 4, the current QDO parametrization yields quite accurate results for dipole, quadrupole, and octupole orders but then exhibits an increasing overestimation of the proportionality constant with increasing multipole rank l. Given that the vdW radius is fixed at its reference value [24], this reflects the fact that the QDO model predicts underestimated multipolar polarizabilities [33]. In order to improve the model, a possible future step is to use highly accurate theoretical or experimental reference data for the quadrupole polarizability instead of the dispersion coefficient C 8 in the parametrization scheme, which would increase the accuracy of the higher-order α l,QDO values. Consequently, such new parametrization would yield smaller values of A QDO l , providing better agreement with the reference data and improving the relations between multipole polarizabilities and the equilibrium distance in vdW-bonded atomic dimers. Figure 4 also shows that the A QDO l curve produces a slight deviation in the general trend of the reference proportionality coefficients at l = 2. This kink stems from the nonmonotonic behavior of the multipole contributions to the exchange energy, where J ex 2(d−q) is larger than J ex 1(d−d) (see our discussion at the end of Sec. III B).
One possible explanation for the difference between the polarizabilities of the QDO model and the hydrogen atom is the contribution of excitations to continuum states in the latter case. The QDO model has no continuum states and can only effectively describe such excitations. Despite the observed deviations in the higher-order multipole polarizabilities and the corresponding proportionality coefficients, the QDO approach allowed us to verify the scaling relation R vdW ∝ α 2/7(l+1) l , which is remarkably valid for atoms.

IV. DISCUSSION AND SUMMARY
We have presented a coarse-grained description of the repulsive force due to the Pauli principle and attractive dispersion forces between two closed-shell atoms or molecules. Our formalism is based on two interacting quantum Drude oscillators, for which the dispersion and exchange-repulsion energies up to an arbitrary order in the multipole expansion of the Coulomb potential were derived. The obtained formulas can be employed for constructing and rationalizing effective interaction potential models, as well as for finding new scaling laws between electronic and geometric properties of atoms and molecules.
As a practical illustration of our theory, we investigated a mutual compensation between the repulsive exchange and attractive dispersion forces for each term in the multipole expansion. The results confirm and extend the recently proposed relations [46] between atomic multipole polarizabilities, α l , and the van der Waals radius, R vdW . The generalized scaling law, R vdW = A l α 2/7(l+1) l , is compelling because it connects an electronic response property of a single atom (atomic multipole polarizability) with the equilibrium distance in a homonuclear dimer.
Let us enumerate some of the potential applications of the formulas presented in this paper and possible future research directions: (i) First and foremost, the relation between atomic vdW radius and atomic polarizabilities, R vdW = A l α 2/7(l+1) l , dispenses with the need to indirectly measure vdW radii. Once the polarizability is calculated for a free atom or an atom in a molecule or material, the vdW radius can be computed from the formula above. Subsequently, the vdW radius can be used as a proxy for an atomic size, as an effective radius in interatomic vdW potentials, or in damping functions for vdW-inclusive electronic-structure calculations. We remark that in quantum mechanics, many possible definitions can be made for an effective atomic size. Our derivations provide a definition of the atomic vdW radius in terms of observable quantities-atomic multipole polarizabilities. Obviously, a more detailed comparison of calculated vdW radii with experiment would be welcome by measuring effective vdW radii in a wide set of systems and comparing the measured radii to first-principles calculations and our formulas.
(ii) Our formula allows a straightforward and accurate calculation of atomic multipole polarizabilities from the dipole polarizability. Given α 1 and the universal values of A ref l obtained in this paper, any multipole polarizability can now be calculated as a function of these two parameters. This is especially important given the high computational cost of calculating multipole polarizabilities from first principles of quantum mechanics. Going further, it would be interesting to assess different recursive relations between α l and α l+1 polarizabilities based on the QDO model and the definition of the vdW radius.
(iii) Our analytical results allow calculating multipole polarizabilities α l for an arbitrary value of l. Such data become increasingly important in coarse-grained models, which describe the molecular response by increasingly larger fragments. For example, one might want to describe the response of a protein, where one QDO models the response of each amino acid. Similar to electrostatics, where higher multipoles become of growing importance when increasing the fragment size, the polarization response follows the same trend. Hence we expect our formulas to play a key role in the development of coarse-grained models for chemical and biological matter.
(iv) While most of the results in this paper were presented for homonuclear dimers, an accurate combination rule is already known for computing equilibrium distances, R AB eq , in heteronuclear dimers [46] This formula allows the calculation of equilibrium distances in heteronuclear closed-shell dimers solely based on the knowledge of atomic polarizabilities of each atom. The derivation of such combination rules from first principles requires generalizing the Pauli principle to QDOs with different parameters, and future work on this subject is needed..
(v) The determination of R AB eq for two atoms A and B from their dipole polarizabilities provides a way to construct generalized Tang-Toennies-type potentials [47,[54][55][56] that require only one adjustable parameter: The equilibrium interaction energy E eq . It remains to be investigated whether the asymptotic dispersion coefficients could be connected to E eq , allowing one to construct parameter-free Tang-Toenniestype interatomic potentials for closed-shell systems.
(vi) Last but not least, the relation between R vdW and the polarizability could be used to develop a more general and more accurate parametrization of the QDO model. Namely, the universality of A ref l coefficients holds for atoms, but it is not such a good approximation within the QDO model itself. One could enforce the obtained relation, R vdW = A l α 2/7(l+1) l , using universal values A l , to hold on average for the QDO model during the parametrization procedure. This is a direction of our current study.
Ultimately, the close connection between vdW attraction and Pauli repulsion unveiled in this paper paves the way for the construction of efficient coarse-grained models for the description of the exchange-repulsion interaction in atomic and molecular systems. Together with the well-established success of the QDO model in describing vdW dispersion, our results also provide the basis for constructing consistent and minimally empirical models for interatomic and intermolecular forces.

ACKNOWLEDGMENTS
The authors acknowledge financial support from the European Research Council via ERC Consolidator Grant "BeStMo (725291)" as well as the Luxembourg National Research Fund via the FNR CORE Jr project PINTA(C17/MS/11686718) and AFR Ph.D. Grant No. CNDTEC(11274975).

APPENDIX A: MULTIPOLE EXPANSION OF COULOMB POTENTIAL
Here, we present the contributions V n to the multipole expansion of the Coulomb potential given by Eq. (5).

APPENDIX B: FORCE BALANCE FOR THE QUADRUPOLE-QUADRUPOLE AND OCTUPOLE-OCTUPOLE INTERACTIONS
In order to demonstrate that the force balance is valid for each term in the multipole expansion, we calculate the proportionality coefficients A 2(q−q) ≡ A 2(q−q) (μω, R vdW ) and A 3(o−o) ≡ A 3(o−o) (μω, R vdW ) for the quadrupole-quadrupole and octupole-octupole interactions, which correspond to the ratios R vdW /α 2/21 2 and R vdW /α 1/14 3 , respectively. According to Eq. (6), one has α 2 = (3/4)(h/μω)α 1 and 033181-12 α 3 = (5/4)(h/μω) 2 α 1 . By employing the aforementioned two relations and considering the quadrupole-quadrupole and octupole-octupole terms in Eqs. (29) and (30) 2(q−q) = 2.08 a.u. and A μω 3(o−o) = 2.00 a.u. with a deviation of 6.3% for the quadrupole-quadrupole interaction and 8.7% for the octupole-octupole interaction. In both cases, the mean values obtained form the high-multipole terms (quadrupole-quadrupole and octupole-octupole) differ most from the universal values A 2 = 2.45 a.u. and A 3 = 2.27 a.u. with respect to the ones obtained from dipole-quadrupole and dipole-octupole interactions. Moreover, the error for A μω is bigger than that for A μω 2(q−q) . This means that, as expected, the lower-order contributions in the multipole expansion are more accurate with respect to the higher-order terms.