Many-particle covalency, ionicity, and atomicity revisited for a few simple example molecules

Determination of the microscopic nature of chemical bonding has been regarded as a problem of fundamental significance since the advent of quantum chemistry and solid state physics [1–3]. The qualitative classification of the valence-electrons state character as covalent, ionic or atomic helps to rationalize their overall features and select a detailed approach to analyze their detailed electronic properties. In this respect, the role of interactions and associated with them interelectronic correlations is crucial in discussing the evolution of bonding from either atomic or ionic character to predominantly covalent or band states of valence electrons. The many-electron approaches, such as configuration interaction (CI) [4] and others [5, 6], are particularly well suited for this task.

In this work we follow a different route and employ exact diagonalization ab initio (EDABI) method, combining the second-quantization formulation of quantum many-particle Hamiltonian with a concomitant readjustment of the single-particle wave functions in the correlated state of the system. This allows us to reinterpret some of the chemical bonding characteristics using concepts originating from condensed-matter physics, such as Mott–Hubbard localization. EDABI has been formulated in our group some time ago [7–9] and analyzed extensively in the context of correlated states in small clusters and one-dimensional solid-state systems. Apart from providing rigorous description of selected properties, EDABI has supplied us with the evolution from the atomic- to a coherent-metallic state with decreasing interatomic distance. Also, modeling the metallization of molecular hydrogen solid has revealed a series of discontinuous first-order Mott-type transitions as a function of applied pressure [10–12]. The explicit question we would like to address here is to what extent the concepts essential extended to lattice quantum systems, such a Mottness [13, 14], may also be qualitatively applicable to finite molecular systems. Answering this question forced us to reanalyze the meaning of the two-particle covalency and related to it ionicity factors by starting from an analytic form of many-particle wave function. We suggest that such analysis may be useful in practical treatment of bonding, here carried out in two-atom–molecule situation, to make the discussion analytic and thus provide a degree of clarity.

The structure of the paper is as follows. In section 2 we summarize briefly the EDABI approach. In section 3 we reanalyze the bonding in H₂, HeH⁺, and LiH systems. We also discuss there validity of the concept of atomicity—Mottness, with the help of which we single out the resonant covalency and atomicity factors. This discussion offers a resolution of the longstanding paradox of the increasing covalency with the increasing interatomic distance. Finally, in sections 4 and 5 we overview our approach. Formal details and tabulated values of the calculated microscopic parameters as a function of interatomic distance are provided in appendices A and B.

The EDABI method has been proposed by us and formulated in detail earlier [7, 9]. Below, we provide a brief summary of its main features, as this should by helpful in grasping the essence of our approach which will be needed in a subsequent interpretation of the results regarding many-particle covalency and ionicity, as well as the concept introduced by us of atomicity.

The starting point is the Hamiltonian containing all pairwise interactions in the second-quantized form is

$$\begin{align}\hfill \hat{\mathcal{H}}& ={{\epsilon}}_{a}\sum\limits _{i}{\hat{n}}_{i\sigma }+{\sum\limits _{ij\sigma }}^{\prime }{t}_{ij}\,{\hat{a}}_{i\sigma }^{{\dagger}}\,{\hat{a}}_{j\sigma }+U\sum\limits _{i}{\hat{n}}_{i{\uparrow}}\,{\hat{n}}_{i{\downarrow}}\hfill \\ \hfill & \enspace +\frac{1}{2}{\sum\limits _{ij}}^{\prime }{K}_{ij}{\hat{n}}_{i\sigma }\,{\hat{n}}_{j{\sigma }^{\prime }}-\frac{1}{2}{\sum\limits _{ij}}^{\prime }{J}_{ij}^{\text{H}}\left({\hat{\mathbf{S}}}_{i}\cdot {\hat{\mathbf{S}}}_{j}-\frac{1}{4}{\hat{n}}_{i}{\hat{n}}_{j}\right)\hfill \\ \hfill & \enspace +\frac{1}{2}{\sum\limits _{ij}}^{\prime }{J}_{ij}^{\prime }({\hat{a}}_{i{\uparrow}}^{{\dagger}}{\hat{a}}_{i{\downarrow}}^{{\dagger}}{\hat{a}}_{j{\downarrow}}{\hat{a}}_{j{\uparrow}}+\mathrm{H}.\mathrm{c}.)\hfill \\ \hfill & \enspace +\frac{1}{2}{\sum\limits _{ij}}^{\prime }{V}_{ij}({\hat{n}}_{i\sigma }+{\hat{n}}_{j\sigma })({\hat{a}}_{i\bar{\sigma }}^{{\dagger}}{\hat{a}}_{j\bar{\sigma }}+\mathrm{H}.\mathrm{c}.)+{\mathcal{H}}_{\text{ion\mbox{--}ion}},\hfill \end{align} \tag{ 1 }$$

where H.c. denotes the Hermitian conjugation, ${\hat{a}}_{i\sigma }$ $({\hat{a}}_{i\sigma }^{{\dagger}})$ are fermionic annihilation (creation) operators for state i and spin σ, ${\hat{n}}_{i\sigma }\equiv {\hat{a}}_{i\sigma }^{{\dagger}}{\hat{a}}_{i\sigma }$, and ${\hat{n}}_{i}\equiv {\hat{n}}_{i{\uparrow}}+{\hat{n}}_{i{\downarrow}}\equiv {\hat{n}}_{i\sigma }+{\hat{n}}_{i\bar{\sigma }}$. The spin operators are defined as ${\hat{S}}_{i}\equiv \frac{1}{2}{\sum }_{\alpha \beta }{\hat{a}}_{i\alpha }^{{\dagger}}{\sigma }_{i}^{\alpha \beta }{\hat{a}}_{i\beta }$ with σ_i representing Pauli matrices. The Hamiltonian contains the atomic and hopping parts (∝_a and t_ij , respectively), the so-called Hubbard term ∝U; representing the intra-atomic interaction between the particles on the same atomic site i with opposite spins, the direct intersite Coulomb interaction ∝K_ij , Heisenberg exchange $\propto {J}_{ij}^{\text{H}}$, and the two-particle and the correlated hopping terms ($\propto {J}_{ij}^{\prime }$ and V_ij , respectively). The last term describes the ion–ion Coulomb interaction which is adopted here in its classical form.

We now proceed to definition of two-particle bonding in general situation and within the second-quantization representation. The N-particle state, |Ψ_N ⟩, may be expressed in terms of the N-particle wave function Ψ(r₁, ..., r_N ) and the corresponding field operators $\hat{{\Psi}}({\mathbf{r}}_{1}),\dots ,\hat{{\Psi}}({\mathbf{r}}_{N})$ as

$$\begin{align}\hfill \vert {{\Psi}}_{N}\rangle & =\frac{1}{\sqrt{N!}}\int {\mathrm{d}}^{3}{r}_{1}\dots {\mathrm{d}}^{3}{r}_{N}\,{{\Psi}}_{N}({r}_{1},\dots ,{r}_{N})\hfill \\ \hfill & \quad \times {\hat{{\Psi}}}_{1}^{{\dagger}}({r}_{1})\dots {\hat{{\Psi}}}_{N}^{{\dagger}}({r}_{N})\vert 0\rangle ,\hfill \end{align} \tag{ 2 }$$

with |0⟩ being the universal vacuum state in the Fock space (for pedagogical exposition see, e.g., [15]). Here we employ a short-hand notation r_i ≡ (r_i , σ_i ), where σ_i = ±1 is the spin quantum number. We can revert this relation to determine the wave function Ψ_N (r₁, ..., r_N ), namely

$$\begin{equation}{{\Psi}}_{\alpha }({r}_{1},\dots ,{r}_{N})=\frac{1}{\sqrt{N!}}\langle 0\vert {\hat{{\Psi}}}_{1}({r}_{1})\dots {\hat{{\Psi}}}_{N}({r}_{N})\vert {\lambda }_{\alpha }\rangle ,\end{equation} \tag{ 3 }$$

where |λ_α ⟩ is the eigenstate for which the wavefunction Ψ_α i explicitly determined. For spin-conserving interaction, α = (σ₁, σ₂, ..., σ_N ) is fixed N-spin-configuration. In effect, we determine |λ_α ⟩ states around ground eigenstate |λ_α ⟩ ≡ |λ_min⟩. Hamiltonian (1) is used to obtain the eigenstates which for two-electron H₂ system are discussed analytically in appendix A.

Since we focus explicitly on the two-site systems, the set of microscopic parameters (_a , t, U, K, J, J', and V) is defined through integrals of orthogonalized single particle basis functions, {w_i (r)}, used next to define the field operators in turn needed to construct the Hamiltonian (1). They are defined briefly first [7, 9], whereas the values of the microscopic parameters are defined in appendix B, starting from the nonorthogonal basis set of adjustable Slater functions. Namely, the orthogonalized atomic (Wannier) orbitals for H₂ molecule the 1s are defined via Slater orbitals {ψ_i (r)} in the usual manner

$$\begin{equation}{w}_{i\sigma }(\mathbf{r})=\beta [{\psi }_{i\sigma }(\mathbf{r})-\gamma {\psi }_{j\sigma }(\mathbf{r})],\end{equation} \tag{ 4 }$$

where σ ± 1 is spin quantum number i ≠ i = 1, 2, and the coefficients β and γ take the form

$$\begin{equation}\begin{cases}\beta =\frac{1}{\sqrt{2}}\sqrt{\frac{1+\sqrt{1-{S}^{2}}}{1-{S}^{2}}},\quad \hfill \\ \gamma =\frac{S}{1+\sqrt{1-{S}^{2}}},\quad \hfill \end{cases}\end{equation} \tag{ 5 }$$

so that β² + γ² = 1. The Slater orbitals ${\psi }_{i}(r)\equiv \sqrt{\frac{{\alpha }^{3}}{\pi }}\,\mathrm{exp}\left(-\alpha (\mathbf{r}-{\mathbf{R}}_{i})\right)$, with α being the inverse size of the orbital, are to be readjusted during the ground-state-energy minimization in the correlated state. Optimization over parameter α is motivated by the circumstance that the selected single-particle basis {ψ_i (r)} is never complete in the quantum mechanical sense and thus such a procedure allows for a better estimate of the ground-state energy. This allows for the orbital-size adjustment in the interacting environment of remaining particles. In brief, the standard procedure of determining the quantum-mechanical state of the system is inverted in the sense that we first diagonalize many-particle Hamiltonian for fixed values of the microscopic parameters and, subsequently, readjust the wave function (inverse size, α⁻¹) in a recurrent fashion. The whole procedure is schematically illustrated by a flowchart composing figure 1.

Zoom In Zoom Out Reset image size

The selected a single-particle basis for H₂ is composed of four orthogonalized wave functions with indices i = 1, 2 enumerating hydrogen atoms and σ = ±1 for each i. Thus, the truncated field operator takes the form

$$\begin{equation}\hat{{\Psi}}(\mathbf{r})=\sum\limits _{i=1}^{2}{w}_{i}(\mathbf{r}){\chi }_{\sigma }(\mathbf{r}){\hat{a}}_{i\sigma }.\end{equation} \tag{ 6 }$$

In that situation, the two-particle wave function is defined in accordance with equation (3), namely

$$\begin{equation}{{\Psi}}_{\alpha }({\mathbf{r}}_{1},{\mathbf{r}}_{2})=\frac{1}{\sqrt{2}}\langle 0\vert {\hat{{\Psi}}}_{1}({\mathbf{r}}_{1}){\hat{{\Psi}}}_{2}({\mathbf{r}}_{2})\vert {\lambda }_{\alpha }\rangle ,\end{equation} \tag{ 7 }$$

where |λ_α ⟩ is the eigenstate (expressed in second quantization representation). Note that this method of approach allows to determine both the ground-state and the lowest excited states for a single optimal value of α.

Note that here the subscripts 1 and 2 of $\hat{{\Psi}}(\mathbf{r})$ contain both site and spin indices for brevity of notation. Parenthetically, one may generalize the above definition to the case of multiple (n) bonds (with n ⩾ 1) as

$$\begin{equation}{{\Psi}}_{\text{a}}({\mathbf{r}}_{1},\dots ,{\mathbf{r}}_{2n})=\frac{1}{\sqrt{2n!}}\langle 0\vert \prod\limits _{i=1}^{n}{\hat{{\Psi}}}_{\text{i}}({\mathbf{r}}_{i})\prod\limits _{j=n+1}^{2n}{\hat{{\Psi}}}_{j}({\mathbf{r}}_{j})\vert {\lambda }_{\alpha }^{(n)}\rangle .\end{equation} \tag{ 8 }$$

In this manner the double (n = 2) and triple (n = 3) bonds can be defined, albeit numerically only and in more complex situations, e.g., in the case of carbon–carbon bonds. This scenario is not addressed here. Instead, we focus on the covalency and ionicity, as well as introduce atomicity + covalency factor, all for selected two-electron systems. However, we discuss first the inherent paradox of the increasing covalency with the increasing interatomic distance.

3.1. Covalent bonding and ionicity in H₂ case

With the help of analysis presented in appendix A for the H₂ case, one can write down explicitly the two-electron wave function in the ground state for H₂ molecule (n = 1 case). The lowest-energy spin-singlet state is of the form

$$\begin{align}{{\Psi}}_{0}({\mathbf{r}}_{1},{\mathbf{r}}_{2})& =\frac{2(t+V)}{\sqrt{2D(D-U+K)}}{{\Psi}}_{\text{c}}({\mathbf{r}}_{1},{\mathbf{r}}_{2})\\ & \quad -\frac{1}{2}\sqrt{\frac{D-U+K}{2D}}{{\Psi}}_{\text{i}}({\mathbf{r}}_{1},{\mathbf{r}}_{2}),\end{align} \tag{ 9 }$$

where the covalent (Ψ_c) and ionic (Ψ_i) components, read

$$\begin{align}\hfill {{\Psi}}_{\text{c}}({\mathbf{r}}_{1},{\mathbf{r}}_{2})& =\left[{w}_{1}({\mathbf{r}}_{1}){w}_{2}({\mathbf{r}}_{2})+{w}_{1}({\mathbf{r}}_{2}){w}_{2}({\mathbf{r}}_{1})\right]\hfill \\ \hfill & \quad \times \left[{\chi }_{{\uparrow}}({\mathbf{r}}_{1}){\chi }_{{\downarrow}}({\mathbf{r}}_{2})-{\chi }_{{\downarrow}}({\mathbf{r}}_{1}){\chi }_{{\uparrow}}({\mathbf{r}}_{2})\right],\hfill \end{align} \tag{ 10 }$$

$$\begin{align}\hfill {{\Psi}}_{\text{i}}({\mathbf{r}}_{1},{\mathbf{r}}_{2})& =\left[{w}_{1}({\mathbf{r}}_{1}){w}_{1}({\mathbf{r}}_{2})+{w}_{2}({\mathbf{r}}_{1}){w}_{2}({\mathbf{r}}_{2})\right]\hfill \\ \hfill & \quad \times \left[{\chi }_{{\uparrow}}({\mathbf{r}}_{1}){\chi }_{{\downarrow}}({\mathbf{r}}_{2})-{\chi }_{{\downarrow}}({\mathbf{r}}_{1}){\chi }_{{\uparrow}}({\mathbf{r}}_{2})\right],\hfill \end{align} \tag{ 11 }$$

with

$$\begin{equation}D\equiv \sqrt{{(U-K)}^{2}+16{(t+V)}^{2}}.\end{equation} \tag{ 12 }$$

The ratio of the coefficients in (9) provides us with the relative ratio of covalency to ionicity in the ground-state spin-singlet configuration. The spin-singlet part is the same in both equations (10) and (11). Explicitly, the covalency and ionicity coefficients (factors) are defined as

$$\begin{equation}{\gamma }_{\text{c}}=\frac{16{(t+V)}^{2}}{16{(t+V)}^{2}+{(D-U+K)}^{2}},\end{equation} \tag{ 13 }$$

and

$$\begin{equation}{\gamma }_{\text{i}}=\frac{{(D-U+K)}^{2}}{16{(t+V)}^{2}+{(D-U+K)}^{2}},\end{equation} \tag{ 14 }$$

respectively, so that the condition γ_i + γ_c = 1 holds. The factor γ_i asymptotically approaches zero in the limit of large interatomic distance (R → ∞), as expected. However, γ_c → 1 for R → ∞ which represents an unphysical behavior [16–18]. This last feature will be discussed in detail below.

The obtained formulas are interpreted as follows. First, the coefficients γ_c and γ_i in the wave function (9) depend on all the interactions which are present in (1), i.e., they contain the effects electronic correlations. Second, the wave functions (10) and (11) take formally the Heitler–London form, but they are self-consistently optimized in the correlated state (their size α⁻¹ is adjustable). Thus, the present formulation in its simplest from contains a semiquantitatively correct behavior in the large-R limit, as is demonstrated explicitly in figure 2.

Zoom In Zoom Out Reset image size

In figure 2, we display the H₂ binding energy and have compared our EDABI calculated value with the results of CI and restricted Hartree–Fock (RHF) analysis. Note difference between the results for small R < R_bond (at minimum), as EDABI method provides slightly lower energies compared to those of full CI. This behavior should not influence the subsequent discussion in the large-R limit, which concerns us mainly here. Nonetheless, it is worth noting both CI and EDABI are variational approaches and perhaps our optimization of the wave-function size at small R is as important as the inclusion of the higher excited states. The binding energy is defined as E_bind = λ₅ − 2E_H (λ₅ is defined in appendix A), where E_H = −1 Ry is the energy of 1s state in atomic hydrogen. Next, we define the bonding and ionicity as the corresponding ratios of coefficients in equation (9), cf figure 3. We note that the covalency increases with the increasing interatomic distance at the expense of ionicity. However, this apparent inconsistency ignores the possibility of incipient atomicity of the Mott–Hubbard type, i.e., the tendency toward localization of electrons on parent atoms with increasing R (called briefly the Mottness). The Mott-type criterion for the localization of electron on H⁺ ion (i.e., formation of renormalized atomic states) takes the form 2|t + V|/(U − K) = 1. This condition expresses the fact that the of bare kinetic energy is then equal to the effective repulsive Coulomb interaction (U − K). In the strong correlation limit, the ratio is below unity, meaning this repulsive interaction becomes predominant (note that in the strict atomic limit, t + V ≡ 0 whereas U − K = 1.25 Ry then). The regime of strong-correlations (Mottness) is marked explicitly in figures 3 and 4. It specifies a gradual evolution toward the atomic state. Namely, the shaded area should be regarded as the regime with steadily increasing atomicity of the electronic states with increasing R. Thus, the question of unphysically increased covalency for R > R_Mott is resolved in a natural manner as within the shaded area the covalency, γ_c, is composed of a sum of true (resonant) covalency ${\bar{\gamma }}_{\text{c}}\to 0$ and atomicity γ_a → 1 as R → ∞ (see the discussion below).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Correlation effects and incipient Mottness

The general meaning of the Mott (or Mott–Hubbard) effects is as follows. In condensed-matter physics the criterion takes the form of U ≃ W [20], where W ≡ |∑_j(i) t_ij | is the bare bandwidth. The transition takes the form of often discontinuous metal–insulator transition for odd integer number of relevant valence electrons per atom. In molecular system, such as H₂, the HOMO–LUMO splitting ≃U − K must overcome the effective interatomic hopping amplitude W = 2|t|. For Hamiltonian (1), the Mott–Hubbard criterion takes then the form r ≡ (U − K)/(2|t + V|) ≃ 1 so that both the correlated hopping and intersite Coulomb interaction contribute, in addition to t and U. In the present situation, the criterion separates only qualitatively the regime of strong correlations (r ⩾ 1) from that with moderate to weak correlations (r < 1). Various versions of the criterion have been shown in figure 5, depending on the theoretical model selected. Namely, the uppermost curve (in green) provides the criterion for the Hubbard model, which does not yield any Mottness point in the present situation. On the other hand, both the model with V = 0 and the full model (represented by the starting Hamiltonian (1)), are almost identical and yield the critical interatomic distance for localization R = R_Mott ≃ 2.3a₀, i.e., well above the R_bond ≃ 1.43a₀. Note also that even for equilibrium distance R = R_bond the hopping/interaction ratio is about $\sim 0.5$, i.e., the electrons are moderately correlated.

Zoom In Zoom Out Reset image size

To complete the picture, we have also plotted in figure 6 the antiferromagnetic kinetic-exchange integral J_kex = 4(t + V)²/(U − K) versus the direct (Heisenberg) ferromagnetic value J^H, both as a function of relative distance R/a₀. The situation is that J_kex > J^H for any distance R and this is the reason for the spin-singlet configuration of H₂ in the ground state. In brief, electrons hoping ('resonating') between the sites, possible only in the total spin-singlet state |λ₅⟩, contribute essentially to the bonding.

Zoom In Zoom Out Reset image size

To verify the conceptual validity of the introduced Mott threshold for atomicity onset, we have plotted in figure 7 the Slater-orbital size α⁻¹ as a function of R. Upon crossing the threshold R_Mott, α⁻¹ indeed approaches rapidly with the further increasing R the 1s atomic size value a₀ = 0.53 Å. Instead, the main physical process contributing to the bonding are the virtual process between the sites. In effect, the ionicity and covalency factors lose their principal meaning for R ≫ 2.3 Å.

Zoom In Zoom Out Reset image size

In conclusion, the dominant covalent character of H₂ molecule has a well defined meaning for R ≃ R_bond, as it is twice as large as the corresponding ionicity factor. However, this decomposition loses gradually its principal meaning as R increases and crosses beyond R_Mott = 2.3a₀. The ground state energy evolves slowly, but steadily toward, the atomic-limit value. Note also that the Hartree–Fock analysis (cf figure 2) provides unphysical results as this critical value of R is crossed. This means that, in the regime of large interatomic distance, the role of correlation becomes essential. In effect, our analysis is applicable then and can be systematically extended numerically by, e.g., enriching the single-particle basis. It would be also of general interest to ask if those concepts could be tested quantitatively by putting H₂ molecules on surfaces of other systems which would stretch the hydrogen-molecule size beyond the Mott–Hubbard threshold. Obviously, the analysis should then incorporate also the presence of the external surface potential of the substrate. However, this type of analysis goes beyond our goals here.

3.3. Physical reinterpretation of atomicity, covalency, and ionicity: resonant covalency

In order to provide a purely physical reinterpretation of covalency and ionicity we note that the form (3) of the covalent part contains sum of static products of the single-particle wave functions located on the sites 1 and 2 and their reverse; this is due to their indistinguishability in the quantum mechanical sense. On the contrary the coefficients γ_c and γ_i contain also virtual intersite processes depicted schematically in figure 8. In other words, the former factor contains a degree of atomicity in its static form, whereas the latter encompasses true dynamic virtual (hopping) processes of quantum-mechanical mixing. The question is how to separate those two factors into atomicity and resonance covalency parts in an analytic way.

Zoom In Zoom Out Reset image size

To answer this dilemma we propose its following resolution. The allowed local (site) states are |0, i⟩, |↑, i⟩, |↓, i⟩, and |↓, ↑, i⟩, i.e., the empty, single occupied with spin σ =↑ or ↓, or the double atomic occupancies. Therefore, using the following identities

$$\begin{equation}\vert 0,i\rangle \langle 0,i\vert +\sum\limits _{\sigma }\vert \sigma ,i\rangle \langle \sigma ,i\vert +\vert {\uparrow},{\downarrow},i\rangle \langle {\uparrow},{\downarrow},i\vert =\mathbb{I},\end{equation} \tag{ 15 }$$

and its equivalent second-quantized from involving site occupancies

$$\begin{equation}\langle (1-{\hat{n}}_{i{\uparrow}})(1-{\hat{n}}_{i{\downarrow}})\rangle +\langle {\hat{n}}_{i{\uparrow}}(1-{\hat{n}}_{i{\downarrow}})\rangle +\langle {\hat{n}}_{i{\downarrow}}(1-{\hat{n}}_{i{\uparrow}})\rangle +\langle {\hat{n}}_{i{\uparrow}}{\hat{n}}_{i{\downarrow}}\rangle =1.\end{equation} \tag{ 16 }$$

Noting that the probability of empty atomic configuration is equal to that doubly occupied, i.e., physically corresponding to the electron–hole symmetry in condensed-matter systems, we obtain the formula for single-electron occupancy in the final form [22]

$$\begin{equation}\nu \equiv \sum\limits _{\sigma }\langle {\hat{n}}_{i\sigma }(1-{\hat{n}}_{i\bar{\sigma }})\rangle =1-2{d}^{2}.\end{equation} \tag{ 17 }$$

Explicitly, we propose to decompose single-occupancy probability ν in the following manner

$$\begin{equation}\nu \equiv a+c=1-2{d}^{2},\end{equation} \tag{ 18 }$$

where a is the atomicity, and c is called the resonant (true) covalency, and ${d}^{2}\equiv \langle {\hat{n}}_{i{\uparrow}}{\hat{n}}_{i{\downarrow}}\rangle$ denotes atom double occupancy probability. Now, the resonant covalency describes the degree of mixing due to the virtual hopping admixture to the frozen (atomic) configuration (cf figure 8(a)). In the strong-correlation limit (r > 1), it can be defined as c ≡ [|t + V|/(U − K)]² and expresses the contribution of the processes (a) to the two-particle wave function in the second order [23, 24] as expressed by ratio of virtual (double hopping, forth and back) process to the Coulomb interaction change in the intermediate step. Therefore, the atomicity is evaluated as

$$\begin{equation}a=\nu -c=1-2{d}^{2}-c.\end{equation} \tag{ 19 }$$

In the equilibrium state of H₂, the resonant covalency reads c ≃ 0.8, whereas atomicity a ≃ 0.1 is practically negligible. Conversely, with increasing R, c decreases quite rapidly and approaches zero, whereas a → 1, as anticipated.

Finally, in figure 9 we provide another characteristic containing atomicity, namely the R dependence d². This double occupancy probability can also characterize the ionicity. The last formula shows that the atomicity is complete when d² = 0 and then ν = a = 1 (i.e., for R ≫ R_Mott). In the other words, the customarily, defined by (13) covalency, associated with the wave function (3), contains both atomicity and true covalency. For R ≫ R_Mott it involves mainly atomicity with a small admixture of c and d². In this manner, the unphysical increase of γ_c with increasing R is resolved. In brief, fundamentally, we define the resonant (true) covalency c as proportional to the inverse Mottness, i.e.,

$$\begin{equation}\text{true}\enspace \text{covalency}=1/\text{Motness}\quad \mathrm{o}\mathrm{r}\quad c\equiv 1/4r.\end{equation} \tag{ 20 }$$

Zoom In Zoom Out Reset image size

In conclusion, based on our analysis of H₂ molecule we suggest that the covalency definition through the values of γ_c is not conceptually precise, whereas the ionicity is properly accounted for either by γ_i or 2d². Additionally, in this way the redefined covalency is complementary to the Mottness and vice versa.

3.4. LiH and HeH⁺ cases

We now apply the concepts introduced above for H₂ molecule to the cases of LiH and HeH⁺. In figures 10 and 11 we display the binding energies versus interatomic distance for HeH⁺ molecular ion and LiH molecule, respectively. In the former case, the two 1s electrons are regarded as core electrons. Effectively, LiH is regarded as a molecule composed of one 2s electron due to Li and 1s electron due to H, with their orbitals adjustable when the interactions are included. Qualitatively, the character of these curves is similar to those of H₂, depicted in figure 2. The quantitative factors are different though and, in particular, the bond length is slightly larger than that in H₂ case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To characterize further those two cases we have plotted in figures 12 and 13 the covalency and ionicity factors for those two systems, respectively. Note that for LiH the ionicity is predominant in a wide range of R, whereas the opposite is true for HeH⁺. The difference arises from the circumstance that in LiH case the orbital size of 2s electron is decisively larger and has a higher energy leading to predominantly ionic configuration ∼Li^+0.9H^−0.9. In HeH⁺, molecular ion the bonding is largely covalent due to the fact that both two 1s² He electrons hop (mix) with the H⁺ state with no electrons in the corresponding 1s state.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One specific feature of those two systems should be noted, which is illustrated in figures 14 and 15, where the optimized sizes of the relevant orbitals has been shown. Namely, the size of 1s orbital of the He and 2s orbital of Li are strongly renormalized, the former largely expanded, whereas the latter contracted. The principal cause of this effect is the electronic correlation induced by the strong intraatomic (Hubbard) interaction ∼U. As this interaction in He is reduced by the flow of electron to the H⁺ site, it is not so in the case of Li, where presence of the hydrogen electron strongly enhances the role of the interaction. In spite of those differences, both systems exhibit similar span of covalency regime. On the contrary, the incipient Mottness appears for larger distance R_Mott and this is presumably due to a larger renormalized-orbital size for Li. As can be seen from literature [25] and from our results here, HeH⁺ is largely covalent and the whole analysis of a and c factors can be repeated here without any qualitative difference.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In table 1 we display the binding energies of the molecules H₂, HeH⁺, LiH, regarded here as testing ground of our approach. For that reason we compare the obtained results with those deducted from other methods and with use of a richer single-particle basis. Even though our results are quantitatively not too accurate, they are obtained with simplest nontrivial basis, i.e., 1s states only for H₂ and HeH⁺ cases, and with addition of 2s states on Li the LiH case. The EDABI results can be improved in a straightforward manner at the expense of computational resources. However in such a situation our following next discussion of bonding would be purely numerical. In other words, we accept the lower accurateness of E_G value within our method to allow for the analytic character of the subsequent discussion. One can find more accurate value of ground state energy for HeH⁺ [29].

Table 1. Binding energy for H₂, HeH⁺, and LiH molecules (in eV).

Method	H2	HeH+	LiH
EDABI	−4.0749	−1.5803	−1.6537
Full CI	−4.3824	−1.6849	−1.8846
RHF	−3.5963	−1.4839	−1.3616
Reference values	−4.3821 [26]	−2.0542 [27]	−1.3606 [28]

We now compare results for those three model systems qualitatively. First, in figure 16 we plot relative contributions of the covalency and ionicity factors for the two-particle ground state (left), as well as a schematic size of the molecular orbitals relative to their original (atomic) size (right). The atomicity factor is not quantified at this stage. In the first two of them, the dynamics is solely due to 1s electrons, whereas in the LiH case the 1s² configuration of electrons is frozen on Li and the whole dynamics is due to 1s–2s H–Li mixing and the corresponding interactions. This is the reason why LiH is largely ionic, whereas the remaining two are predominantly covalent, as illustrated in figure 13. One sees that the covalency in HeH⁺ is larger than that for H₂ molecule, a rather unexpected intuitively result.

Zoom In Zoom Out Reset image size

A separate discussion should be concerned with other overall properties of the systems studied. In table 2 the Mott (or Mott–Hubbard) critical distance R_Mott (in the units of a₀) is provided. This distance should be compared with the bond length R_bond calculated (cf table 3) according to three independent methods: EDABI, full CI, and RHF methods, respectively. We see that in each case R_bond is decisively lower than R_Mott. This means that the Mott-type boundary can be crossed only in the situation when the molecules are further apart, i.e., obtained artificially when, e.g., they are placed on surfaces with an external force elongating them. Particularly favorable situation occurs when molecules are placed in the environment with a large dielectric constant, as then interaction weakens and the bond length increases. Clearly, then the whole analysis must be revised and a realistic configuration with inclusion of appropriate external (surface) potential. We believe that the essential features of our analysis should survive when the molecule is placed in such environment, i.e., in a potential stretching equally both atoms.

In table 1 the binding energies are listed and compared with those from other methods. These numerical results present probably the weakest point of our EDABI method, since the corresponding values obtained are not very accurate. Nevertheless, we do not consider our method as a practical computing tool. Instead our main aim here was to extend, albeit at best in a semiquantitative manner, the basis for multi-electron covalency and ionicity, enriched by the concept of atomicity, all induced by the electronic correlations. Obviously, the approach can be extended in a straightforward manner by enlarging the single-particle basis and applied for the systems with larger atoms. Both of these factors have been considered by us before [7, 9] for model systems, with one limitation, that we have not analyzed there the bonding properties. This analysis should be explored further along the lines discussed here.

Finally, in table 4 we list the most important microscopic parameters in the equilibrium state. A more detailed analysis of those is presented in appendix B. The values of Coulomb-interaction parameters will be reduced by the dielectric constant factor if system under consideration is placed on surface of an insulating material. This should rescale all the parameter values accordingly.

Table 4. Calculated EDABI values of the equilibrium parameters for LiH, HeH⁺, and H₂ molecules.

Parameter	LiH	H2	HeH+
EG (eV)	−217.45	−31.198	−79.675
U1sLi (eV)	49.062	N/A	N/A
U1sH (eV)	18.559	22.490	19.592
U1sHe (eV)	N/A	N/A	14.925
U2sLi (eV)	12.784	N/A	N/A
t (eV)	−21.150	−9.9049	−15.674
K (eV)	14.368	13.007	11.374
α1sH $({a}_{0}^{-1})$	1.035	1.194	1.240
α1sHe $({a}_{0}^{-1})$	N/A	N/A	2.095
α2sLi $({a}_{0}^{-1})$	1.329	N/A	N/A

The reason for selecting the three systems analyzed here is caused by the circumstance that HeH⁺ is strongly covalent, LiH strongly ionic, and H₂ can be placed in between them. On example of the last of them our novel concept of atomicity and resonant covalency have been proposed.

The introduced here atomicity for the case of molecular system (corresponding to Mott–Hubbard localization effects in periodic systems) amounts to specifying a gradual transformation from molecular to atomic language in describing their electronic states, as a function of interatomic distance. This changeover is the basic feature and is associated with the essential change in regarding those particles as evolving within indistinguishable (molecular) character and acquiring eventually the form of distinguishable (atomic) states.

One must also underline that the concept of atomicity here is quantitative in nature. This is because the Mott–Hubbard localization concept in condensed-matter systems [13, 20] appears usually as a first-order transition, requiring the energy equality of the two macro configuration (delocalized, localized) at this phase transition. Here the evolution may be regarded as a supercritical behavior at best [13, 30, 31]. However, the antiferromagnetic kinetic exchange survives even when the states are becoming orbitally distinguishable [32].

Certainly, a further insight is required to quantify the present discussion for more complex systems. The present concepts are proposed to clarify the obviously unphysical behavior of the increasing covalency with the increasing interatomic distance. As far as we are aware of, this inconsistency, although intuitively understandable, has not been discussed explicitly in the quantum-chemical literature. Also, the emerging atomicity here squares well with the Mott's original argument [20] that the metallic (covalent) state of electrons in a periodic system is ruled out at (semi)macroscopic interatomic distances.

This work was supported by Grants OPUS Nos. UMO-2018/29/B/ST3/02646 and UMO-2021/41/B/ST3/04070 from Narodowe Centrum Nauki. We would like to thank Prof. Ewa Brocławik and Dr Mariusz Radoń for discussions and criticism. We are also grateful to Andrzej Biborski and Andrzej P Kadzielawa for making available to us their QMT library.

All data that support the findings of this study are included within the article (and any supplementary files).

Starting from the orthogonalized restricted basis, we define the field operators as

$$\begin{equation}{\hat{\psi }}_{\sigma }(\mathbf{r})={w}_{1}(\mathbf{r}){\chi }_{\sigma }(1){\hat{a}}_{1\sigma }+{w}_{2}(\mathbf{r}){\chi }_{\sigma }(2){\hat{a}}_{2\sigma },\end{equation} \tag{ A.1 }$$

$$\begin{equation}{\hat{\psi }}_{\sigma }^{{\dagger}}(\mathbf{r})={w}_{1}^{\ast }(\mathbf{r}){\chi }_{\sigma }(1){\hat{a}}_{1\sigma }^{{\dagger}}+{w}_{2}^{\ast }(\mathbf{r}){\chi }_{\sigma }(2){\hat{a}}_{2\sigma }^{{\dagger}},\end{equation} \tag{ A.2 }$$

or, in compact notation, as

$$\begin{equation}\hat{\psi }(\mathbf{r})\equiv \left(\begin{matrix}\hfill {\hat{\psi }}_{{\uparrow}}(\mathbf{r})\hfill \\ \hfill {\hat{\psi }}_{{\downarrow}}(\mathbf{r})\hfill \end{matrix}\right).\end{equation} \tag{ A.3 }$$

In the above, ${\hat{a}}_{i\sigma }$ and ${\hat{a}}_{i\sigma }^{{\dagger}}$ are electron annihilation and creation operators in the state w_iσ (r) ≡ w₁(r)χ_σ (1). Also, as we restrict here to s-orbital systems, the molecular (Wannier) functions can be taken as real if the condition J' = J^H holds. Using the representation (A.3), we obtain Hamiltonian (1) with the microscopic parameters expressed through the Slater orbitals and coefficients β and γ (cf equation (5)), or explicitly through inverse orbital size α and interatomic distance R (see e.g. [7, 9]). The relevant physical quantities may be thus obtained as a function of R, with the orbital parameter α optimized in each case.

To obtain the ground state energy E_G for fixed R, the Hamiltonian (1) is diagonalized. This is carried out by making use of the global symmetry respecting two-particle states, leading to block-diagonal many-body Hilbert space, with specified values of the total spin, S, and its z-component, S^z , as well with transposition antisymmetry preserved. In effect, one can start the basis of $\left(\genfrac{}{}{0.0pt}{}{4}{2}\right)=6$ following states

$$\begin{equation}\left\{ \begin{array}{ll}\quad & \vert 1\rangle ={\hat{a}}_{1{\uparrow}}^{{\dagger}}{\hat{a}}_{2{\uparrow}}^{{\dagger}}\vert 0\rangle ,\\ \quad & \vert 2\rangle ={\hat{a}}_{1{\downarrow}}^{{\dagger}}{\hat{a}}_{2{\downarrow}}^{{\dagger}}\vert 0\rangle , \\ \quad & \vert 3\rangle =\frac{1}{\sqrt{2}}({\hat{a}}_{1{\uparrow}}^{{\dagger}}{\hat{a}}_{2{\downarrow}}^{{\dagger}}+{\hat{a}}_{1{\downarrow}}^{{\dagger}}{\hat{a}}_{2{\uparrow}}^{{\dagger}})\vert 0\rangle , \\ \quad & \vert 4\rangle =\frac{1}{\sqrt{2}}({\hat{a}}_{1{\uparrow}}^{{\dagger}}{\hat{a}}_{2{\downarrow}}^{{\dagger}}-{\hat{a}}_{1{\downarrow}}^{{\dagger}}{\hat{a}}_{2{\uparrow}}^{{\dagger}})\vert 0\rangle , \\ \quad & \vert 5\rangle =\frac{1}{\sqrt{2}}({\hat{a}}_{1{\uparrow}}^{{\dagger}}{\hat{a}}_{1{\downarrow}}^{{\dagger}}+{\hat{a}}_{2{\downarrow}}^{{\dagger}}{\hat{a}}_{2{\uparrow}}^{{\dagger}})\vert 0\rangle , \\ \quad & \vert 6\rangle =\frac{1}{\sqrt{2}}({\hat{a}}_{1{\uparrow}}^{{\dagger}}{\hat{a}}_{1{\downarrow}}^{{\dagger}}-{\hat{a}}_{2{\downarrow}}^{{\dagger}}{\hat{a}}_{2{\uparrow}}^{{\dagger}})\vert 0\rangle . \end{array}\right.\end{equation} \tag{ A.4 }$$

The first three are the spin-triplet states with S^z = +1, −1, 0, whereas the next three are inter- and intra-site singlets, respectively. The triplet state does not hybridize with other states and provides three 1 × 1 irreducible blocks with eigenvalues λ₁ = λ₂ = λ₃ = ₁ + ₂ + K − J^H. The remaining three singlet states compose the 3 × 3 block, so the Hamiltonian in that Fock subspace takes the form

$$\begin{equation}\hat{\mathcal{H}}=\left(\begin{matrix}\hfill {\epsilon}+K+{J}^{\text{H}}\hfill & \hfill 2(t+V)\hfill & \hfill 0\hfill \\ \hfill 2(t+V)\hfill & \hfill 2{\epsilon}+J+U\hfill & \hfill \frac{1}{2}({U}_{1}-{U}_{2})\hfill \\ \hfill 0\hfill & \hfill \frac{1}{2}({U}_{1}-{U}_{2})\hfill & \hfill 2{\epsilon}+U-{J}^{\text{H}}\hfill \end{matrix}\right),\end{equation} \tag{ A.5 }$$

where ≡ (₁ + ₂)/2 and U ≡ (U₁ + U₂)/2. This formulation allows to apply this formalism to both H₂ (where U₁ = U₂ = U and ₁ + ₂ = ), and to HeH⁺ and LiH, where those simplifications are not met due to inequivalent atoms involved.

In the case of H₂, the eigenvalues of equation (A.5) take the form

$$\begin{equation}{\lambda }_{4,5}=2{\epsilon}+\frac{1}{2}(K+U)+J\pm \frac{1}{2}D,\end{equation} \tag{ A.6 }$$

$$\begin{equation}{\lambda }_{6}=2{\epsilon}+U-J,\end{equation} \tag{ A.7 }$$

with $D\equiv {[{(U-K)}^{2}+16{(t+V)}^{2}]}^{\frac{1}{2}}$. The corresponding eigenstates are

$$\begin{equation}\vert {\lambda }_{4,5}\rangle \equiv \vert {\lambda }_{\pm }\rangle =\frac{1}{{[D(D\pm U\mp K)]}^{\frac{1}{2}}}[(4(t+V))\vert 4\rangle \pm (D\pm U\mp K)\vert 5\rangle ],\end{equation} \tag{ A.8 }$$

where, for simplicity, we have defined the atomic-limit energy as the reference point, = 0. We note that the eigenstates |λ_4,5⟩ are superposed of the symmetric ionic state |5⟩ and covalent part |4⟩. The state |λ₅⟩ is the ground state as the λ₅ eigenvalue is the lowest one. In the limit U ≫ |t + V|, the λ₅ eigenvalue reads

$$\begin{equation}{\lambda }_{5}\simeq 2{\epsilon}+{J}^{\text{H}}+K-\frac{4{(t+V)}^{2}}{U-K}.\end{equation} \tag{ A.9 }$$

The last term on the right-hand side of equation (A.9) is the so-called kinetic-exchange contribution. It competes with ferromagnetic Heisenberg exchange ∼J^H. In similar manner, the two-particle states for HeH⁺ and LiH are obtained, except that in those two cases, the diagonalization of the Hamiltonian matrix (A.5) cannot be carried out analytically, since the ₁ ≠ ₂ and U₁ ≠ U₂. The singlet state |λ₅⟩ is elaborated further throughout the main text.

In tables B1–B3 we provide relevant quantities and microscopic parameters versus R, obtained within EDABI scheme for the three systems discussed in main text, i.e., H₂, HeH⁺, and LiH.

Table B1. Ground state energy and microscopic parameters for H₂ molecule (in eV).

R (a0)	EG/N		t	U	K	J	V
0.5	−657	−14.29	−31.75	30.06	19.43	0.43	−0.28
1	−1486	−22.51	−15.95	25.29	15.43	0.36	−0.19
1.5	−1558	−23.84	−9.21	22.08	12.67	0.29	−0.16
2	−1516	−23.41	−5.79	19.96	10.75	0.23	−0.16
2.5	−1461	−22.56	−3.84	18.61	9.34	0.18	−0.16
3	−1418	−21.67	−2.62	17.81	8.24	0.13	−0.16
3.5	−1389	−20.86	−1.82	17.38	7.35	0.09	−0.16
4	−1373	−20.15	−1.26	17.18	6.60	0.06	−0.15
4.5	−1365	−19.52	−0.86	17.09	5.95	0.04	−0.14
5	−1362	−18.98	−0.59	17.05	5.40	0.02	−0.12

Table B2. Ground state energy and microscopic parameters for HeH⁺ molecular ion (in eV).

R (a0)	EG/N	H	He	t	UH	UHe	K	V
0.5	−27.84	−22.42	−14.32	−24.07	22.35	36.70	11.17	−0.99
1	−38.54	−32.81	−32.43	−18.96	20.57	21.85	8.48	−0.75
1.5	−39.84	−34.59	−33.13	−15.49	19.54	14.62	6.72	−0.60
2	−39.63	−34.10	−29.57	−13.15	18.95	11.10	5.66	−0.50
2.5	−39.38	−32.20	−25.93	−11.55	18.61	9.39	4.95	−0.44
3	−39.21	−29.67	−22.86	−10.38	18.41	8.56	4.54	−0.39
3.5	−37.68	−27.20	−20.45	−9.45	18.30	8.15	4.26	−0.36
4	−39.09	−25.08	−18.62	−8.64	18.23	7.96	4.05	−0.33
4.5	−38.97	−23.37	−17.23	−7.90	18.19	7.86	3.86	−0.31
5	−38.92	−22.01	−16.16	−7.20	18.17	7.82	3.69	−0.31

Table B3. Ground state energy and microscopic parameters for LiH molecule (in eV).

R (a0)	EG/N	H	2sLi	t	UH	U2sLi	K	V
1	−98.21	−43.25	−39.32	−55.89	38.29	33.08	21.19	−1.70
1.5	−105.89	−44.47	−40.24	−48.02	35.13	26.34	19.79	−1.52
2	−108.10	−45.402	−40.97	−37.97	32.02	22.01	17.10	−1.02
2.5	−108.88	−45.73	−41.97	−30.93	28.89	17.98	16.77	−0.89
3	−109.29	−46.12	−42.53	−24.43	24.80	14.04	15.12	−0.78
3.5	−109.35	−45.93	−41.55	−19.89	22.08	11.35	13.29	−0.69
4	−109.29	−45.81	−41.53	−14.05	20.30	9.44	11.98	−0.51
4.5	−109.04	−45.50	−41.44	−11.23	19.19	8.41	10.13	−0.38
5	−108.92	−45.75	−41.39	−8.09	18.09	7.75	9.48	−0.32
5.5	−108.71	−44.96	−40.91	−5.48	18.25	7.48	6.32	−0.28
6	−108.32	−44.49	−40.76	−4.01	17.98	7.25	4.01	−0.26

Note that the value of |t| is comparable to U in the limit R < R_bond and diminishes spectacularly when R > R_bond (i.e. in the strong-correlation regime).

The RHF and CI computations were carried out using the GAMESS code and the 6–31 G basis set to represent the Slater functions. Numerical accuracy for the EDABI calculations is 10⁻⁴ a₀ for R and 10⁻⁵ eV for energy, respectively.