$H_2$ and $(H_2)_2$ molecules with an ab initio optimization of wave functions in correlated state: Electron-proton couplings and intermolecular microscopic parameters

The hydrogen molecules $H_2$ and $(H_2)_2$ are analyzed with electronic correlations taken into account between the $1s$ electrons exactly. The optimal single-particle Slater orbitals are evaluated in the correlated state of $H_2$ by combining their variational determination with the diagonalization of the full Hamiltonian in the second-quantization language. All electron--ion coupling constants are determined explicitly and their relative importance is discussed. Sizable zero-point motion amplitude and the corresponding energy are then evaluated by taking into account the anharmonic contributions up to the ninth order in the relative displacement of the ions from their static equilibrium value. The applicability of the model to the solid molecular hydrogen is briefly analyzed by calculating intermolecular microscopic parameters for $2 \times H_2$ rectangular configurations.

The hydrogen molecules H2 and (H2) 2 are analyzed with electronic correlations taken into account between the 1s electrons in an exact manner. The optimal single-particle Slater orbitals are evaluated in the correlated state of H2 by combining their variational determination with the diagonalization of the full Hamiltonian in the second-quantization language. All electron-ion coupling constants are determined explicitly and their relative importance is discussed. Sizable zero-point motion amplitude and the corresponding energy are then evaluated by taking into account the anharmonic contributions up to the ninth order in the relative displacement of the ions from their static equilibrium value. The applicability of the model to the solid molecular hydrogen is briefly analyzed by calculating intermolecular microscopic parameters for 2 × H2 rectangular configuration, as well its ground state energy.

I. MOTIVATION
The few-site models of correlated fermions play an important role in singling out, in an exact manner, the role of various local intra-and inter-site interactions against the hopping (i.e., containing both the covalent and the ionic factors) and thus, in establishing the optimal correlated state of fermions [1][2][3][4][5][6][7][8] on local (nanoscopic) scale. The model has also been used to obtain a realistic analytic estimate of the hydrogen-molecule energies of the ground and the excited states in the correlated state [9]. For this purpose, we have developed the so-called ED-ABI method, which combines E xact Diagonalization in the Fock space with a concomitant Ab I nitio determination of the single-particle basis in the Hilbert space. So far, the method has been implemented by taking only 1s Slater orbitals, one per site [10]. The method contains no parameters; the only approximation made is taking a truncated single-particle basis (i.e., one Slater orbital per site) when constructing the field operator, that in turn is used to derive the starting Hamiltonian in the second-quantization representation. This Hamiltonian represents an extended Hubbard Hamiltonian, with all two-site interactions taken into account and the solution comprises not only the exact eigenvalues of the few-site Hamiltonian, but also at the same time evaluation of the adjustable single-particle wave functions in the correlated state. Also, the thermodynamic properties calculated rigorously exemplify [11,12] the low-and high-energy scales, corresponding to spin and local charge fluctuations respectively. The former represents precursory magnetic-ordering effect whereas the latter local effects accompanying with the Mott-Hubbard transition. In general, our approach follows the tradition of accounting for interelectronic correlations via the second-quantization procedure, with the adjustment of the single-particle wave functions, contained in microscopic parameters of the starting model, in the correlated state.
The first aim of this paper is to extend a fully microscopic approach established earlier [9,10] and calculate all six possible electron-ion coupling constants for H 2 as a function of the bond length. As a byproduct, we obtain an accurate estimate of the zero-point-motion amplitude and its energy to a high (ninth) order in the relative displacement of the ions. This evaluation shows explicitly the dominant contributions to the vibronic spectra of the molecule. In effect, the work formulates a complete two-site model of correlated states with all the coupling parameters calculated from an ab initio procedure. It also forms a starting point to a full scale dynamic calculations involving a richer basis in the Hilbert space, at least in the adiabatic limit. So, although the importance of the present results to the discussion of H 2 molecule exact evaluation of the ground-state energy is limited, the approach may be extended to treat the molecular solid hydrogen with inclusion of interelectronic correlations. Explicitly, as a starting point we calculate the intermolecular hopping integrals and the principal electron-electron interaction microscopic parameters as a function of intermolecular distance.
A methodological remark is in place here. As we determine the local ion-electron and electron-electron cou-pling parameters, they can be regarded as a starting estimate of those for the bulk solid molecular hydrogen, as we have studied recently a critical pressure of metallization of the atomic solid (Mott insulating) state [13]. The obtained pressure of atomic-hydrogen metallization is about 100 GP a, the value which squares well with observed for the case of fluid molecular hydrogen (140 GP a) [14], although the recent simulations provide quite different values for fluid hydrogen analyzed at high temperature [15]. Obviously, our previous work is not related directly to the molecular-hydrogen metallization in the solid state [16][17][18][19]. So far, we have discussed rigid-lattice properties. We believe that the present results form the first step in incorporating the vibrational spectrum and correlations to extended systems.
The structure of the paper is as follows. Even though the main purpose of the paper is to calculate the local electron-proton and electron-electron coupling constants, for the sake of the completeness, in Sections II & III we reproduce some of the results of [9] and correct some minor errors (cf. also Appendix A). In Sec. IV we define the method of calculation of both the electron-ion (proton) coupling constants (cf. also Appendix B), as well as estimate the zero-point motion to the ninth order vs. the interionic distance. In Sec. V we extend the single-molecule treatment and provide the intermolecular hopping amplitudes and the electronelectron microscopic parameters which may serve for analysis of the solid molecular hydrogen. Section VI contains physical discussion and a brief outlook, where we also refer to the finite-size Quantum Monte Carlo results. In the series of Appendices we provide some analytical details, as they may form analytical basis for the electronlattice coupling supplementing the classic Slater results for the electronic part of H 2 molecule [20].

A. Wannier basis
To describe the behavior of an electron in the system of two ions we start from 1s Slater-type orbitals where α is the inverse size of the orbital. To ensure orthogonality we use Wannier functions which in this case reduce the superposition of the atomic states, i.e., with the mixing parameters where S = S(α, R) ≡ Ψ 1 | Ψ 2 is the atomic functions' overlap. Eqs. (3) ensure both the orthogonality and proper behavior in atomic limit i.e., lim The two-site Hamiltonian with one orbital per site has the general form whereâ iσ andâ † iσ are the fermionic operators of annihilation and creation of the electron with spin σ on 1s orbital of hydrogen atom i ∈ {1, 2}.
The microscopic parameters = T 11 , t = T 12 , U = V 1111 , J = V 1122 , K = V 1212 and V = V 1112 correspond to one-and two-particle interactions [9] T ij = w i | T |w j , (5a) where in atomic units T = − 2 −2/|r − R|, and V = 2/|r − r |. In the Appendix A we provide explicitly the form of microscopic parameters as a function of both intersite static distance R and the inverse wavefunction size α. In what follows we diagonalize first (4), and subsequently optimize the wave functions contained in the microscopic parameters of (4). This program will be carried out systematically in what follows.

C. Exact solution
System described by the Hamiltonian (4) has an exact solution previously studied in detailed in [9]. For twoelectron system (n 1 + n 2 = 2), i.e. with one particle per site, the starting basis is representing the intersite spin-triplet states with eigenvalues E 1 = E 2 = E 3 = 2 + K − J, and representing the spin-singlet states, with the corresponding Hamiltonian matrix involving the matrix elements i| H |j ≡ (H ij ): The state (6f) is an eigenvector of (7) with eigenvalue E 6 = 2 + U − J. The diagonalization supplies us with the two other eigenvectors [1] with eigenvalues where D = (U − K) 2 + 16 (t + V ) 2 . The state |− from (8) is the lowest-energy spin-singlet eigenstate. It is this state, for which we determine explicitly the singleparticle wave function and subsequently, determine the microscopic parameters , t, U , J, K, and V explicitly, all as a function of interionic distance R.

D. Optimization of the atomic basis
The ground-state energy is the energy E − of (9), supplemented with the ion-ion repulsion, i.e. by where (2/R) is represented also in atomic units. As all the microscopic parameters are only a function of the distance R and the inverse wave-function size α, we have E G = E G (α, R). For each distance R, we minimize E G with respect to α, thus closing the solution. Finally, we select R = R B as the equilibrium solution, for which still the zero-point motion has to be taken into account.

III. STATIONARY STATE FOR H2 MOLECULE
In Fig. 1 we plot the energy of H 2 (dimer) versus the distance R. It is crucial that we obtain a local (and Solutions for the states: the spin-triplets (E1 = E2 = E3) and the spin-singlets (E±, E6) versus the interionic distance R. The spin-singlet state |− is the equilibrium state.
global) minimum at R = R B ≡ 0.757Å. This simple result, obtained in [9] differs with respect to the virtually exact solution by Kołos and Wolniewicz [21,22], In Fig. 2 we plot the sequence of the spin-singlets and the spin-triplet states. Parenthetically, the start from second-quantization language allows for evaluation of the ground-state and the lowest excited states, on an equal footing. This feature provides the difference with purely variational calculations in the first-quantization language. Namely, within this basis the spin-singlet state is stable at arbitrary interionic distance R. In Figs. 3 -5 the inverse wave function size α, as well as all the microscopic parameters, are all displayed as a function of the bond length R. One can see that with the increasing R values tends to the proper free-atom limits. Those quantities form an input for the subsequent evaluation of Microscopic parameters , t, U , and K versus average interionic distance R. Note the convergence of the intersite Coulomb repulsion K to the classical value 2/R (dashed line) at R → ∞. The on-site repulsion U reaches also its atomic limit Uat = 1.25 Ry, whereas the hopping parameter t → 0. electron-proton coupling constants discussed next.

IV. ADIABATIC APPROXIMATION FOR THE ELECTRON-ION COUPLING
Our principal aim here is to extend our previous model [9,10] by allowing the ions oscillate around the equilibrium positions. Thus the interionic distance R is taken now in the form where δR is responsible for the zero-point motion. The electronic part of the ground-state energy is expanded next on δR in terms of a Taylor series, which to the ninth order reads where E and E (1) B = 0, whereas all the remaining terms but for the energy E B describe the oscillations (cf . Table II for numerical values). We have modified the Hamiltonian (4) accordingly by taking into account δR, i.e., where δH is the additional term. Also, H simplifies to the form where Ξ = { , t, U, J, K, V } andÔ i are the corresponding operator parts of Hamiltonian: the two-and four-operator terms of (4) standing next to the respective microscopic parameter (for exampleÔ =n 1 +n 2 ). With the Hamiltonian in this form we now have the energy change due to the change of the microscopic parameters where i are bosonic creation and annihilation operators of the system deformation and the set {ξ i } defines a new set of microscopic parameters -the electron-ion coupling constants. They can be derived by differentiation in a way similar to that of [23,24] (cf. Appendix B for details).
Averages (17) calculated in the ground-state versus distance R. They are of the order of unity.
Coupling constants ξU , ξK , ξJ and ξV versus intersite distance R The shift of the ions changes the system properties in the following manner where the average Ô i 0 = − Ô i − is taken with respect to the ground state. In effect, we obtain The R dependence of the parameters Ô i given by (17) is displayed in Fig. 6. As they are of the order of unity, the principal factor determining the relative strength of the coupling constants are provided by the parameters {ξ i } displayed in Figs. 7 and 8. At the equilibrium bond distance marked by the vertical line, the largest values are (on the absolute scale) those coming from modulation of the hopping parameter (ξ t ) and the change of intersite Coulomb interaction (ξ K ). The first of the two has been included in Su, Schrieffer and Heeger model [25]. The second may play an important role in the high-T C superconductivity [24]. Also, we see that the so-called Holstein coupling [26] is not important if calculated near the hydrogen-molecule equilibrium state.
We determine the value of δR by minimizing the total energy of the system: where where the ionic momentum δP is evaluated via Heisenberg principle and M a.u.
≈ 1836.15267m e is the mass of the proton.

V. EVALUATION OF THE MICROSCOPIC PARAMETERS FOR THE TWO-MOLECULE SYSTEM
We extend our approach by considering a system of two H 2 dimers at the relative distance a from each other FIG. 9. The system of two H2 molecules at the relative distance a. The hopping integrals ti are marked next to the respective dashed lines. Note that the orthogonalization procedure for four sites produces a different basis than that obtained in (3). (cf. Fig. 9). We calculate the respective hopping integrals, where t 12 should approach t defined in (A1b) for large a, as well as the new single-particle energy should again converge to previously obtained value (A1a). We determine all the two-particle interaction integrals, thus going beyond the Hubbard model solved in [5]. Additionally, in Table I  Explicitly, in Fig. 10 we display the intermolecular dependence of single-particle parameters. For the distances a 2 a 0 the hoppings t 13 and t 14 can be regarded as small on the scale t 12 = t. Hence, the system in solid will preserve its molecular character, with no magnetism involved even though we have nominally one electron per FIG. 12. Two-particle microscopic parameters: intramolecular spin-exchange J12 and correlated hopping V12, as well as the intermolecular parameters J13, J14, V13 and V14 for two H2 molecules system vs. intermolecular distance a. Note that all the intermolecular parameters converge to zero quickly.
atom. In other words, the lowest band will be full and no simple-minded Hubbard subband (HOMO-LUMO) picture in the ground state appears. In Fig. 11 we compare the relative values of intramolecular (U , K 12 = K) versus intermolecular (K 13 , K 14 ) Coulomb interactions. Again, intramolecular interactions dominate for a 2 a 0 . From Figs. 10 and 11 it follows then, that in the insulating (molecular-crystal) state virtual hopping processes in a similar manner to the kinetic exchange, will contribute and renormalize the gap between the full-band (valence) and the conduction-band (excited single electron) states. This gap will have the form of the Hubbard gap, as the value of U , corresponding to the transition 2H 2 → H − 2 + H + 2 will have the value U − K 12 ≈ 0.6 Ry, by far the largest energy in the insulating state. For the TABLE I. Numerical values of single-particle energy ( ), the hopping integrals (t αβ ), the on-site Coulomb repulsion U and the intramolecular Coulomb interaction K12 for the two-molecule system, values refer to the points marked in the Fig. 13 FIG. 13. The one-particle microscopic parameters (in Ry) for two H2 molecules system vs. intermolecular distance a and interionic distance R. Note the symmetry of and t14.
As expected for relatively small distances values of t12 and t13 are negative, whereas t14 is positive. When approaching point (0, 0) all parameters diverge to minus (t12 and t13) or plus ( and t14) infinity. The explicit values of the marked points are given in Tab. I sake of completeness, we have plotted in Fig. 12 the remaining interaction parameters: the exchange integrals, intra-(J 12 ) and inter-molecular (J 13 and J 14 ), as well as the correlated hopping amplitudes: V 12 and (V 13 and V 14 ), respectively.
In Fig. 14 we show the difference between the energies of the (H 2 ) 2 system and that of the free molecules (per where E (H2) 2 is the energy of the (H 2 ) 2 system and E B = −2.29587 Ry is the energy of single molecule. The equilibrium parameters are ∆E H2 = −0.01129 Ry and a = 4.5 a 0 . Those results are in agreement with the earlier estimates [28]. The stability of hydrogen molecular clusters were also studied in [29,30]; our approach coherently incorporates electronic correlations (a necessity in describing the non-polar systems) into the molecular picture, that plays an important role in view of the existence of the minimas of ∆E H2 (a) curve [28,29].

VI. DISCUSSION OF RESULTS AND OUTLOOK
The evaluation of the global parameters of the system can be summarized as follows: 1. The H 2 binding energy is E B = −2.29587 Ry. This value can be compared to the Kołos-Wolniewicz value [21]: E K-W = −2.349 Ry, which is about 2.26% lower than the value obtained here.
2. The increase of the binding energy here is due to the zero-point motion, which is and is of the order 1.0485% of the binding energy.
3. Whereas the bond length is here R B = 1.43042 a 0 (as compared to R K-W = 1.3984 a 0 , which is 2.29% lower), the zero-point motion amplitude is |δR| = 0.189028 a 0 , a rather large value. Note that the optimal size of the inverse orbit of the 1s component hydrogen orbital is α B = 1.19378 a −1 0 , so that the effective Bohr orbit is a ≡ α −1 = 0.8377 a 0 . The Bohr orbit decrease is due to the increased binding of electron in the molecule ∼ 0.2932 Ry with respect to that in the hydrogen atom. The size a is substantially smaller than that of 1s orbit (1.06 Å) in H atom [9].  Table III. We also provide the second-order coupling constants values at the hydrogen-molecule equilibrium in the Table IV. Our first-principle calculations allow us to claim that the coupling constant appearing in the Holstein model [26] (ξ ) is decisively smaller than those of Su, Schrieffer and Heeger model [25] (ξ t ) as well as of those coming from the intersite direct Coulomb interaction (ξ K ). This should not be surprising, as the dominant coupling parameters represent interatomic-vibration contributions. A separate branch is represented by phonon excitations, but their analysis requires a construction of a spatially extended system of the molecules.
The question is to what extent the calculated local characteristics will represent their local counterparts in the molecular-solid phase. Certainly, the phonons and the molecule-rotational degrees of freedom will represent low-energy excitations. But the zero-point motion energy of the whole molecule should have to be added. Will this provide a reliable description of the molecular or atomic hydrogen even though our accuracy in determining the individual-molecule energy is about 2% higher than the virtually exact value of Kołos and Wolniewicz [21]? One has to check and this task is under consideration in our group. Such consideration must include the inter-molecular hopping integrals t 13 and t 14 . B all of the derivatives are calculated analytically from equation (10). Orders seventh-ninth (marked by an asterisk) were calculated numerically due to complicated analytical expression for groundstate energy. One should also note that the proper method of treating the few-site H 2 -molecule system is the quantum-Monte-Carlo-method [31][32][33]. Nonetheless, our method evaluates both the system energetics and the wave-function renormalization at the same time in the correlated state.
One can also extend the present model of the molecular binding by including also 2s and 2p adjustable hydrogen orbitals in the Hilbert space of the single-particle states via the corresponding Gaussian representation. The first estimate of the 2s-orbital contribution to selected microscopic parameters is briefly discussed in D. Their numerical values are provided in VI. One can see that the basis extension leads to the numerically relevant corrections. This is an additional route to follow, but only after the first-principle calculations of the solid phase along the lines discussed here is undertaken and tested.
Very recently [34], the dynamical mean field theory (DMFT) has been applied to H 2 molecule and its accuracy tested. Our approach in this respect is much simpler, but still provides a comparable accuracy. Also, we have calculated the vibronic coupling constants, which have been determined accurately recently [35]. Those results compare well with our estimates. This circumstance shows again that our method forms a proper starting point for treatment of solid molecular hydrogen, as a correlated state, at least in the insulating phase.

VII. ACKNOWLEDGMENTS
The authors (APK & JS) are grateful the Foundation for Polish Science (FNP) for financial support within the TEAM Project. M.A. is grateful to the Jagiellonian University for hospitality and FNP for a partial support for his visit. J.S. is grateful to the National Science Center (NCN) for the support within the MAE-   For the sake of completeness we express the microscopic parameters defined in (5) in terms of single-particle parameters via (2) where Ξ parameters are with T = − 2 −2/|r − R|, and V = 2/|r − r |. The single-particle parameters read where the overlaps are given by and Ei(x) is the Exponential Integral.
Appendix B: Adiabatic-approximation details For the sake of completeness, we provide also the explicit form of the coupling constants ξ i ≡ dΞ i /dR. For editorial purposes we calculate first the single-particle coupling constants from (A3). They are The corresponding derivatives of the two-particle parameters are with basis mixing-parameters β and γ (cf. (3)) changes where S is the overlap (A4) and δS reads Our final results are These parameters are displayed vs. R in Figures 7 and  8. We ask the question how important is to include the quantum nature of the electronic interaction in our evaluation of zero-point motion energy. Let us consider, following [9], the energy of the ions as where δP and δR are the momentum and position uncertainties, M is ion mass and e its charge. By expressing δP by δR via uncertainty relation δP 2 δR 2 3 4 2 we obtain that we can minimize with respect to δR. We calculate δR 0 a.u.
where a a.u.
We take the mass of the ion M a.u.
≈ 1836.15267m e and the interionic distance R = R B = 1.43017a 0 . The results are presented in the Table V.

Appendix D: Inclusion of 2s orbitals
We would like to estimate the role of higher orbitals both for more realistic description of H 2 systems and future consideration of other elements. We start by taking the 2s Slater-type orbital where α 2s is the inverse wave function size. It is obviously non-orthogonal with the 1s orbital (1)  We perform the orthogonalization by introducing realistic orbital functions [36] χ 1s (r) = Ψ 1s (α, r) , (D3) χ 2s (r) = AΨ 1s (α 2s , r) + BΨ 2s (α 2s , r) , where A and B are mixing parameters obtained via orthonormality conditions χ 2s χ 2s = 1.
One can obtain the exact values for the optimal interionic distance R = R B = 1.43042 and α = α B = 1.19378. The results, together with comparison to the one-orbital case, are listed in Table VI. Note that the new estimates are carried out for the optimal bond length and the inverse wave-functions size for the case of 1s functions only (R = R B = 1.43042 and α = α B = 1.19378, respectively). The values of single-particle microscopic parameters without optimization of 2s orbitals for 1s and 2s band (R = RB, α = αB, α2s = αB/2). Note that last column describes atomic limit, where α → 1, α2s → 0.5. We can observe that the model fulfills requirements, as the Wannier function w 2s i approaches the exact solution of Hydrogen atom.