Ground-state correlation energy of beryllium dimer by the Bethe-Salpeter equation

TL;DR

This paper uses advanced many-body perturbation theory methods, specifically GW and Bethe-Salpeter equations, to accurately compute the ground-state correlation energy of the challenging beryllium dimer, aligning well with experimental data.

Contribution

It presents the first ab initio calculations of Be₂'s potential energy using GW and BSE, demonstrating improved accuracy over previous methods.

Findings

GW corrections improve energy estimates at RPA level

BSE calculations match experimental and high-level theoretical results

Reproduces experimentally observed flattening of the potential

Abstract

Since the '30s the interatomic potential of the beryllium dimer Be$_2$ has been both an experimental and a theoretical challenge. Calculating the ground-state correlation energy of Be$_2$ along its dissociation path is a difficult problem for theory. We present ab initio many-body perturbation theory calculations of the Be$_2$ interatomic potential using the GW approximation and the Bethe-Salpeter equation (BSE). The ground-state correlation energy is calculated by the trace formula with checks against the adiabatic-connection fluctuation-dissipation theorem formula. We show that inclusion of GW corrections already improves the energy even at the level of the random-phase approximation. At the level of the BSE on top of the GW approximation, our calculation is in surprising agreement with the most accurate theories and with experiment. It even reproduces an experimentally observed flattening of the interatomic potential due to a delicate correlations balance from a competition between covalent and van der Waals bonding.

Tables

Table 1. Table 1: Be 2 bond length R e subscript 𝑅 𝑒 R_{e} and energy D e = E ​ ( R e ) subscript 𝐷 𝑒 𝐸 subscript 𝑅 𝑒 D_{e}=E(R_{e}) .

Method	$R_e$ [bohr]	$D_e$ [mHa]
MP2 Klopper and Almlöf (1993); Toulouse et al. (2009)	5.0	-1.99
CCSD Magoulas et al. (2018)	8.37	-0.26
CCSD(T) Mentel and Baerends (2014)	4.64	-3.00
QMC DMC Toulouse and Umrigar (2008)	4.65	-2.82
BSE@$GW$@HF by TF Eq. (1)	4.65	-4.66
Accurate Røeggen and Veseth (2005)	4.63	-4.34
Experiment Merritt et al. (2009)	4.64	-4.24

Equations

$$E_0^{\mathrm{c}}=\frac{1}{2}\left[{\sum }_{i>0}{\Omega }_i-\mathrm{Tr}(A)\right]$$

$$E_0^{\mathrm{c}}=\frac{1}{2\pi }{\int }_0^{i\infty }d\omega \mathrm{Tr}\left[\ln \left(1-w{\chi }^0(\omega )\right)+w{\chi }^0(\omega )\right].$$

$$L=L^0+L^0\Xi L,$$

$$\left(\begin{array}[]{cc}A&B\\ B^{*}&A^{*}\end{array}\right)\left(\begin{array}[]{c}X^{\lambda}\\ Y^{\lambda}\end{array}\right)=\Omega_{\lambda}\left(\begin{array}[]{cc}1&0\\ 0&-1\end{array}\right)\left(\begin{array}[]{c}X^{\lambda}\\ Y^{\lambda}\end{array}\right),$$

$$⟨{\Psi }_0|[{\hat{C}}_t,[\hat{H},{\hat{C}}_{t^{\prime }}^{†}]]|{\Psi }_0⟩$$

$$⟨{\Psi }_0|[{\hat{C}}_t,[\hat{H},{\hat{C}}_{t^{\prime }}]]|{\Psi }_0⟩$$

$$\hat{H}=E_0^0+\sum _{t{t}^{\prime }}{A}_{t{t}^{\prime }}{\hat{C}}_t^{†}{\hat{C}}_{t^{\prime }}+\frac{1}{2}\sum _{t{t}^{\prime }}[B_{t{t}^{\prime }}{\hat{C}}_t^{†}{\hat{C}}_{t^{\prime }}^{†}+B_{t{t}^{\prime }}^{\ast }{\hat{C}}_t{\hat{C}}_{t^{\prime }}],$$

$$\hat{H}=E^{0}_{0}-\frac{1}{2}\mathrm{Tr}(A)+\frac{1}{2}\left(\begin{array}[]{cc}\hat{C}^{\dagger}&\hat{C}\end{array}\right)\left(\begin{array}[]{cc}A&B\\ B^{*}&A^{*}\end{array}\right)\left(\begin{array}[]{c}\hat{C}\\ \hat{C}^{\dagger}\end{array}\right),$$

$${\hat{Q}}_{\lambda }^{†}=\sum _t(X_t^{\lambda }{\hat{C}}_t^{†}-Y_t^{\lambda }{\hat{C}}_t),$$

$$\hat{H}=E_0^0-\frac{1}{2}\mathrm{Tr}(A)+\frac{1}{2}\sum _{\lambda }{\Omega }_{\lambda }+\sum _{\lambda }{\Omega }_{\lambda }{\hat{Q}}_{\lambda }^{†}{\hat{Q}}_{\lambda }.$$

$$E_0=E_0^0+\frac{1}{2}\sum _{\lambda }{\Omega }_{\lambda }-\frac{1}{2}\mathrm{Tr}(A).$$

$$E_0^{\mathrm{corr}}=\frac{1}{2}\sum _{\lambda }({\Omega }_{\lambda }-{\Omega }_{\lambda }^{\mathrm{TDA}}),$$

$$E_0^{\mathrm{corr}}=-\sum _{\lambda }{\Omega }_{\lambda }\sum _t|Y_t^{\lambda }{|}^2,$$

Full text

Introduction.

The beryllium dimer Be2 has a long scientific history, with hundreds of experimental and theoretical investigations Merritt et al. (2009); Røeggen and Almlöf (1996). The first synthesis of Be2 was attempted in the ’30s with no success Herzberg (1929, 1933), while Hartree-Fock (HF) or other theoretical modelling Bartlett Jr. and Furry (1931); Bender and Davidson (1967) found a repulsive ground-state, leading to the conclusion that Be2 does not exist. Later studies Blomberg and Siegbahn (1978); Røeggen and Almlöf (1996) pointed to a possible van der Waals binding with a shallow energy minimum at large ($\sim 5$ Å) distance, while other studies yielded a double minimum, at short and long distance separation Wilson (1983); Røeggen and Almlöf (1996). Be2 remained elusive till the ’70s Brom Jr. et al. (1975), and only in the ’80s first rotovibrational spectra were measured Bondybey (1984) and reliable calculations Harrison and Handy (1983) were made, both pointing to a single short-bond minimum at $\sim 2.5$ Å. Today we have very accurate experiments Merritt et al. (2009) and calculations Røeggen and Veseth (2005) for the Be2 interatomic potential. Nevertheless, Be2 remains a severe workbench to check many-body theories, and this is the purpose of this work.

The description of the correlation energy of molecular dimers along their dissociation path is a theoretical challenge Martin et al. (2016); Helgaker et al. (2008). Full configuration interaction (CI) Helgaker et al. (2008) is the only accurate method Røeggen and Veseth (2005); Røeggen and Almlöf (1996) providing results in very good agreement with experiment, but it is limited to very small molecules. Quantum Monte Carlo (QMC) Martin et al. (2016); Casula et al. (2004, 2005); Toulouse and Umrigar (2008); Marchi et al. (2009); Varsano et al. (2014), both variational (VMC) and diffusion (DMC) are also valid alternatives, and DMC is even exact in systems where the ground-state wavefunction does not present nodes, but they can be also very cumbersome. Density-functional theory (DFT) Martin et al. (2016) is in principle exact for calculating ground-state energies, but standard approximations, like the local-density (LDA) and generalized-gradient (GGA) approximations, have shown their limits on dimers binding energies and lengths Fuchs and Gonze (2002), in particular in the dissociation limit. In recent years time-dependent density-functional theory (TDDFT) Runge and Gross (1984); Gross and Kohn (1985) in the adiabatic-connection fluctuation-dissipation theorem (ACFDT) formalism Langreth and Perdew (1975, 1977) has been considered as one promising approach to improve over approximated DFT Furche (2001); Fuchs and Gonze (2002); Heßelmann and Görling (2011); Ren et al. (2012). In any case it relies on TDDFT approximations for the polarizability $\chi$, such as the random-phase approximation (RPA), the adiabatic LDA (ALDA or TDLDA), or beyond. Nevertheless, the interatomic potential of Be2 continues to be a problem also for TDDFT ACFDT, both in the RPA Nguyen and Galli (2010) and also more advanced approximations Toulouse et al. (2009); Ren et al. (2012).

In this work, we calculate the Be2 interatomic potential in the framework of ab initio many-body perturbation theory using the $GW$ approximation Hedin (1965) and the Bethe-Salpeter equation (BSE) Salpeter and Bethe (1951); Hanke and Sham (1979); Martin et al. (2016). The ground-state correlation energy is calculated by the trace formula (TF) Ring and Schuck (1980); Rowe (1968a)

[TABLE]

where $\Omega_{i}$ are the positive eigenvalues of the full, i.e. beyond the Tamm-Dancoff approximation, BSE equation, and $A$ is the only-resonant part of the BSE excitonic Hamiltonian. This formula was introduced by Sawada to calculate the correlation energy of the electron gas by summing over the contributions from zero-point plasma oscillations (plasmons) Sawada (1957) as an alternative to the Gell-Mann & Brueckner formula which integrates along the adiabatic connection path over the interaction switch-on $\lambda$ parameter Gell-Mann and Brueckner (1957). Later these two formulae were shown to be equivalent within RPA, both on the electron gas Sawada et al. (1957); Fukuta et al. (1964) and for inhomogeneous systems in using the DFT Kohn-Sham (KS) as reference Furche (2008). This allows us to validate our results from Eq. (1), at least at the RPA level, by also evaluating the correlation energy via the RPA ACFDT formula starting from KS eigenstates Langreth and Perdew (1975, 1977); Fuchs and Gonze (2002)

[TABLE]

Here $w(r,r^{\prime})=1/|r-r^{\prime}|$ is the Coulomb interaction, $\chi^{0}(r,r^{\prime},\omega)$ is the Kohn-Sham polarizability, and the integral over $\omega$ is along the positive imaginary axis. Our results show that $GW$ corrections introduce large improvements already at the level of the RPA. The interatomic potential we obtain at the level of the BSE on top of the $GW$ approximation is in surprising agreement with the most accurate calculations and experiment. $GW$+BSE correlations even seem able to describe the unusual shape of the experimental Merritt et al. (2009) interatomic potential, namely a flattening of the Morse potential in the range between 6 and 9 bohr (3 to 4.5 Å) towards an expanded Morse potential which better fits the experimental vibrational spectrum Merritt et al. (2009).

The solution of the historical problem represented by the paradigmatic Be2 interatomic potential shows that the use of Eq. (1) in the BSE framework can reveal an accurate methodology for calculating the ground-state energy and stability of atoms and molecules, solids, and even nuclei, with an important advance in all these fields.

Theory.

Eq. (1) was derived previously in different ways Sawada (1957); Fukuta et al. (1964); Furche (2008). Here we present a modification of the Thouless derivation Ring and Schuck (1980); Brown and Wong (1967) which extends its validity to starting points different from HF towards DFT or $GW$, and to kernels beyond RPA towards BSE or TDDFT kernels. We start from the Bethe-Salpeter equation,

[TABLE]

where $L$ is the two-particle correlation function, $L^{0}=GG$ with $G$ the one-particle Green function, and $\Xi$ the two-particle interaction. This is an exact equation for calculating the excitation spectrum. By knowing the exact $\Xi$ and $L^{0}$ (via $G$), the BSE can be solved for $L$ whose poles $\Omega_{\lambda}=E_{\lambda}-E_{0}$ and their associated residuals provide the neutral (optical) excitation energies and oscillator strengths. In practice, approximations are unavoidable. We consider a quite general case of an approximated electronic structure, e.g. HF, $GW$ or DFT, with real energies $\epsilon_{i}$ and orthonormal wavefunctions $\phi_{i}(r)$ used to build $G$ and $L^{0}$. And we consider an approximated static kernel, e.g. the kernel TDH $i\Xi_{ijkl}=w_{iljk}$ (with $w$ the bare Coulomb interaction), or the TDHF $i\Xi_{ijkl}=w_{iljk}-w_{ilkj}$, or the $GW$+BSE $i\Xi_{ijkl}=w_{iljk}-W_{ilkj}$ (with $W$ the screened Coulomb interaction), or even the TDLDA. Under these conditions the Bethe-Salpeter equation can be reduced to the well known Ring and Schuck (1980); Martin et al. (2016); Rowe (1968b); Casida (1995, 1996) RPA equation

[TABLE]

with the Hermitian matrix $A_{tt^{\prime}}=A_{php^{\prime}h^{\prime}}=(\epsilon_{p}-\epsilon_{h})\delta_{pp^{\prime}}\delta_{hh^{\prime}}+i\Xi_{php^{\prime}h^{\prime}}$ and the symmetric matrix $B_{tt^{\prime}}=B_{php^{\prime}h^{\prime}}=-i\Xi_{pp^{\prime}hh^{\prime}}$ indexed by particle-hole transition indices $t=\{ph\}$ from an occupied state $h$ to an empty state $p$, $\hat{c}^{\dagger}_{p}\hat{c}_{h}$, on the orthonormal basis set $\phi_{i}(r)$. Eq. (4) provides a full spectrum of excitations $|\Psi_{\lambda}\rangle$, both the excitation energies $\Omega_{\lambda}=E_{\lambda}-E_{0}$ (with respect to the ground state energy $E_{0}$) and the eigenvectors $(X^{\lambda}\,Y^{\lambda})$. This spectrum constitutes an approximation to the exact spectrum. We now introduce boson Rowe (1968b) transition operators $\hat{C}_{t}$ replacing the $ph$ operator bilinears, $\hat{c}^{\dagger}_{p}\hat{c}_{h}\to\hat{C}^{\dagger}_{t}$, i.e. fulfilling exact boson canonical commutation relations, $[\hat{C}_{t},\hat{C}_{t^{\prime}}]=0,[\hat{C^{\dagger}}_{t},\hat{C}^{\dagger}_{t^{\prime}}]=0,[\hat{C}_{t},\hat{C}^{\dagger}_{t^{\prime}}]=\delta_{tt^{\prime}}$. This is not the case for the fermion bilinears, $[\hat{c}_{p}\hat{c}^{\dagger}_{h},\hat{c}^{\dagger}_{p}\hat{c}_{h}]\neq\delta_{pp^{\prime}}\delta_{hh^{\prime}}$. So the boson operators can be seen as an approximation to the fermion bilinears. We now express the full many-body Hamiltonian $\hat{H}$ in terms of the boson operators imposing the condition that we get, by construction, the same excitation spectrum of Eq. (4). This is achieved if Rowe (1968b)

[TABLE]

Up to quadratic terms only Jancovici and Schiff (1964) the Hamiltonian is written

[TABLE]

where the constant $E^{0}_{0}$ is approximated by the expectation value of the Hamiltonian, $E^{0}_{0}=\langle\Phi_{0}|\hat{H}|\Phi_{0}\rangle$, over the zero-order ground-state Slater determinant $\Phi_{0}$ constructed with the occupied wavefunctions $\phi_{h}(r)$ 111More exactly $E^{0}_{0}=\langle\Psi^{0}_{0}|\hat{H}|\Phi^{0}_{0}\rangle$ where $\Psi^{0}_{0}$ is the ground-state killed by the boson operators, $\hat{C}_{t}|\Psi^{0}_{0}\rangle=0$. Since $\hat{C}_{t}\simeq\hat{c}^{\dagger}_{p}\hat{c}_{h}$, the 0-order ground state $\Phi_{0}$, killed by the $\hat{c}^{\dagger}_{p}\hat{c}_{h}|\Phi_{0}\rangle=0$, is an approximation for $\Psi_{0}^{0}$, and $\hat{C}_{t}|\Phi_{0}\rangle\simeq 0$ is approximated to zero.. For HF orbitals $\phi^{\mathrm{HF}}(r)$, the constant is the HF ground-state energy $E^{0}_{0}=E^{\mathrm{HF}}_{0}$, whereas in the case of DFT the constant is the sum of the kinetic, external, Hartree and exchange operators evaluated on the Kohn-Sham wavefunctions, $E^{0}_{0}=E^{\mathrm{DFT}}-E_{\mathrm{xc}}+E_{\mathrm{EXX}}$. The Hamiltonian can be recast into a matrix form,

[TABLE]

and the solutions $(X^{\lambda}\,Y^{\lambda})$ to Eq. (4), subject to the orthonormality condition $\sum_{t}(X^{\lambda*}_{t}X^{\lambda^{\prime}}_{t}-Y^{\lambda*}_{t}Y^{\lambda^{\prime}}_{t})=\delta_{\lambda\lambda^{\prime}}$, allows us to define new boson operators $\hat{Q}_{\lambda}$ by a Bogoliubov transformation of the $\hat{C}_{t}$,

[TABLE]

which diagonalizes the Hamiltonian

[TABLE]

The operators $\hat{Q}_{\lambda}$ act on the (approximated) ground $|\Psi_{0}\rangle$ and excited $|\Psi_{\lambda}\rangle$ states, $\hat{Q}^{\dagger}_{\lambda}|\Psi_{0}\rangle=|\Psi_{\lambda}\rangle$, $\hat{Q}_{\lambda}|\Psi_{0}\rangle=0$. And the expectation value of the Hamiltonian Eq. (5) over the excited states $|\Psi_{\lambda}\rangle$ provides, by construction, the excitation energies $\Omega_{\lambda}$ with respect to the ground state. Finally, the expectation value of the Hamiltonian Eq. (5) over $|\Psi_{0}\rangle$ provides the total energy of the ground state including correlation (within the stated approximations),

[TABLE]

Thus, for a starting HF electronic structure, the two terms after the constant $E^{0}_{0}=E^{\mathrm{HF}}_{0}$ in Eq. (6) provide the correlation energy Eq. (1), alternatively recast in the form

[TABLE]

or in the form

[TABLE]

where $\Omega_{\lambda}^{\mathrm{TDA}}$ are the eigenvalues of Eq. (4) in the Tamm-Dancoff approximation (TDA), i.e. taking the matrix $B=0$. Eqs. (7) and (8) show more clearly that the ground state correlation energy arises only beyond the TDA approximation and from terms in $B$ and $Y$. The TDA is only able to introduce correlations in the excited states.

To summarize, Eq. (6) is an expression for the ground-state total energy $E_{0}$ including correlation at the same level of approximation (e.g. RPA or TDHF, $GW$+BSE, etc.) as taken for the $A$ and $B$ matrices, and so at the same level of approximation of the associated excitation spectrum $\Omega_{\lambda}$. In contrast to the ACFDT formalism, which provides an in principle exact formula for correlations within TDDFT, Eq. (6) is only an approximate expression since the beginning that relies on the validity of the boson approximation between operators, $\hat{C}_{t}\simeq\hat{c}^{\dagger}_{p}\hat{c}_{h}$, and the validity of the killing condition on the zero-order Slater determinant ground-state, $\hat{C}_{t}|\Phi_{0}\rangle\simeq 0$. However for both formulae most critical are the approximations on the kernel (RPA and beyond) and on the starting electronic structure (LDA or else). Comparing the merits of these various approximations is still in its infancy for atoms and molecules Olsen and Thygesen (2014); Holzer et al. (2018).

In order to provide the best comparison with previous literature and with ACFDT results, we use the same large cc-pV5Z Gaussian basis set Dunning (1989) adopted in Ref. Toulouse et al. (2009), together with the auxiliary cc-pV5Z-RI basis set Hattig (2005) in a Coulomb-fitting, resolution-of-identity approach. Input Kohn-Sham or Hartree-Fock eigenstates were calculated by the NWChem Valiev et al. (2010) package, whereas many-body calculations were performed with the Fiesta Jacquemin et al. (2015); Duchemin et al. (2018) code. For improved accuracy Jacquemin et al. (2015); Gui et al. (2018) we performed ev$GW$ calculations, namely partially self-consistent $GW$ on the eigenvalues only. $GW$ corrections were calculated explicitly for all 4 occupied and 14 empty states, while corrections to higher states were extrapolated by a rigid scissor-operator shift from the last calculated level. This provided a result converged up to 0.2 mHa.

Results.

In Fig. 1 we present the interatomic potential of Be2 as derived from recent accurate experimental rotovibrational spectra Merritt et al. (2009), and from the most accurate configuration interaction (CI) calculation Røeggen and Veseth (2005) which is in good agreement with the experiments. On the same figure, we present the theoretical curve that we calculated at the level of the direct RPA approximation on top of DFT in the PBE Perdew et al. (1996) approximation (RPA@PBE) by the TF formula Eq. (6) or equivalently Eq. (8) which, we checked, provide the same result within numerical precision. We show also the same RPA@PBE curve calculated by the ACFDT formula. They are in perfect agreement. We can compare our RPA@PBE curves with those calculated by other authors Nguyen and Galli (2010); Toulouse et al. (2009) with the ACFDT formula. Although we used exactly the same calculation parameters reported in Ref. Toulouse et al. (2009), our result differ from theirs, in particular for the binding energy, due to the basis set superposition error that they mentioned and which they removed by counterpoise. We prefer not to use the counterpoise method because some works Sheng et al. (2011); Mentel and Baerends (2014) showed it is not justified in the case of systems where van der Waals interaction can be important. For this reason our result is closer, surprisingly, to the one of Ref. Nguyen and Galli (2010) obtained by a much different implementation, namely plane-waves and norm-conserving pseudopotentials. Nevertheless, all RPA@PBE calculations present qualitatively the same fake feature, a “bump” (a maximum) at 6–7 Bohr, in contrast to experiment and to the accurate CI calculation. Note that the bump is completely absent in the original DFT PBE interatomic potential which is strongly attractive everywhere with a deep minimum and large binding energy (see Ref. Ren et al. (2012)). Finally, we present curves calculated at the same level of direct RPA but using a $GW$ corrections on top of PBE. The most important point is that there is a large improvement in RPA@$GW$@PBE with respect to RPA@PBE. The “bump” is wiped out. The binding energy is also improved, though still underestimated. The bonding length is now overestimated.

The bump at 6–7 Bohr appears in TF or ACFDT RPA calculations on top of DFT PBE (as well as LDA, see Ref. Nguyen and Galli (2010)) but it is conjured by the exact exchange (EXX) term evaluated on DFT PBE eigenstates (but not HF) and wiped out by a $GW$ calculation with self-consistency on the eigenvalues only. This can be appreciated also in Fig. 2 where we present an RPA calculation on top of HF (RPA@HF), both our TF result and the ACFDT result of Ref. Nguyen and Galli (2010). Both results present a correct attractive behaviour without fake bumps, although the original HF interatomic potential (see Fig. 2) is repulsive. Nevertheless the dissociation curve is too flat (though less in our calculation than in Ref. Toulouse et al. (2009)), and the binding energy is severely underestimated. However, the introduction of $GW$ corrections introduces an important improvement also when starting from HF. The RPA@$GW$@HF curve (Fig. 2) is much closer to the accurate and the experimental result, though it presents a $\sim$10% overestimate of both the binding energy and the bond length. Finally Fig. 2 shows the curve calculated by solving the Bethe-Salpeter equation on top of a $GW$ electronic structure and starting from HF. The result is, beyond any expectation, impressively good. The equilibrium length is within 0.02 bohr from the exact value, and the binding energy is overestimated by only 0.3 mHa (see Table 1). These values are even more accurate than the QMC value Toulouse and Umrigar (2008), outperforming quantum chemistry methods like MP2 and even CCSD as well Lotrich et al. (2005); Toulouse et al. (2009); Klopper and Almlöf (1993); Magoulas et al. (2018). Notice also that any change we have introduced to the standard procedure, like solving the BSE directly on top of HF, or considering an unscreened kernel as in TDHF, immediately destroys the good agreement, even providing results that are much worse in some cases. Our BSE result seems better than the most advanced ACFDT TDDFT approximations of Ref. Ren et al. (2012), like the RPA+ and RPA+SOSEX, which are unbound, and also the RPA+rSE, which provides a good binding energy and bond length but still presents a bump at 6–7 Bohr. The range-separated hybrid (RSH) functional of Ref. Gerber and Ángyán (2005); Toulouse et al. (2009, 2010) is certainly the approximation that is closer to the physics underlying the $GW$ approach in which correlations are addressed by the introduction of screening. RSH+MP2 or RSH+RPAx tries to mimic this physics by introducing two ranges of different screened exchange. In Be2 it obtains a clear improvement toward the good shape of the potential with no bump, but the binding energy and the bond length are not yet sufficiently reproduced. This can only be obtained by an approach presenting continuous variation of the screening at all ranges, like in GW+BSE.

Finally, the BSE on top of the $GW$ approximation result even seems able to reproduce an unusual feature which has been pointed out in the accurate experiment of Ref. Merritt et al. (2009). In the Inset of Fig. 2 we show a zoom on the region 6 to 9 bohr ($\sim$3 to 4.5 Å) where van der Waals (vdW) interaction effects should enter into play and where, in the past, it was conjectured the existence of a vdW secondary (or even the only main) minimum. In this region the experimental curve presents an evident flattening which makes the interatomic potential deviate from the simple Morse potential, towards a more complex expanded Morse oscillator (EMO, see Ref. Merritt et al. (2009) and its Fig. 3). The EMO potential shape seems essential to best fit the experimental vibrational spectrum of Ref. Merritt et al. (2009) and seems a particularity of this dimer in which there is competition between covalent and vdW interactions. The vdW interaction is unable to produce a secondary minimum, but does have an influence on the shape of the interatomic potential, distorting it from the simple Morse oscillator towards a flattening. Our BSE@$GW$@HF curve also presents this flattening, though shifted with respect to experiment. It is nevertheless impressive that we were able to describe such a delicate balance between covalent and vdW bonding which only arises from correlations. This means that the Bethe-Salpeter equation provides an at least qualitatively good description of the high order correlation effects.

Conclusions.

We have calculated the Be2 interatomic potential along its dissociation path by an approach within ab initio many-body perturbation theory relying on the trace formula, as an alternative to the TDDFT ACFDT formalism. Our approach has been validated against TDDFT ACFDT calculations already presented in the literature. The introduction of $GW$ corrections already largely improves the shape of the interatomic potential and also the binding energy and the bond length. At the level of the BSE we have obtained a Be2 interatomic potential in very good agreement with experiment and with accurate CI calculations. We were even able to reproduce a flattening observed in the experimental potential which results from a delicate balance of correlations due to a competition between covalent and van der Waals bonding.

Acknowledgements.

Acknowledgements.

This work is dedicated to the memory of János G. Ángyán, whom we acknowledge for inspiration and useful discussions. We thank Maria Hellgren for useful discussions and Ronald Cox for critical reading of our manuscript.