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

We analyze two-particle binding factors of H2, LiH, and HeH+ molecules/ions with the help of our original exact diagonalization ab intio (EDABI) approach. The interelectronic correlations are taken into account rigorously within the second quantization scheme for restricted basis of renormalized single-particle wave functions, i.e., with their size readjusted in the correlated state. This allows us to determine the many-particle covalency and ionicity factors in a natural and intuitive manner in terms of the microscopic single-particle and interaction parameters, also determined within our method. We discuss the limitations of those basic characteristics and introduce the concept of atomicity, corresponding to the Mott and Hubbard criterion concerning localization threshold in many-particle systems. This addition introduces an atomic ingredient into the electron states and thus removes a spurious behavior of covalency with the increasing interatomic distance, as well as provides a more complete physical interpretation of bonding.


Introduction
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,2,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 Intio (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,8,9] and analyzed extensively in the context of correlated states in small clusters and onedimensional 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,11,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 twoatom-molecule situation, to make the discussion analytic and thus provide a degree of clarity.
The structure of the paper is as follows. In Sec. 2 we summarize briefly the EDABI approach. In Sec. 3 we reanalyze the bonding in H 2 , 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 Sec. 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-B.

EDABI method and many-particle bonding
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 secondquantized form iŝ where H.c. denotes the Hermitian conjugation,â iσ (â † iσ ) are fermionic annihilation (creation) operators for state i and spin σ,n iσ ≡â † iσâ iσ , andn i ≡n i↑ +n i↓ ≡n iσ +n iσ . The spin operators are defined asŜ i ≡ 1 2 αβâ † iα σ αβ iâ iβ with σ i representing Pauli matrices. The Hamiltonian contains the atomic and hopping parts (∝ a and t i j , 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 i j , Heisenberg exchange ∝ J H i j , and the two-particle and the correlated hopping terms (∝ J i j and V i j , 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 1 , . . ., r N ) and the corresponding field operatorŝ Ψ(r 1 ), . . .,Ψ(r N ) as 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 1 , . . ., r N ), namely where |λ α is the eigenstate for which the wavefunction Ψ α i explicitly determined. For spin-conserving interaction, α = (σ 1 , σ 2 , . . ., σ 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 2 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 2 molecule the 1s are defined via Slater orbitals {ψ i (r)} in the usual manner where σ ± 1 is spin quantum number i i = 1, 2, and the coefficients β and γ take the form , with α being the inverse size of the orbital, are to be readjusted during the ground-state-energy minimization Figure 1. Flowchart of the EDABI method. The method is initialized by selection of a trial single-particle basis of wave functions, {w i (r)}, and subsequent diagonalization the resulting many-particle Hamiltonian. Optimization of the single-particle states leads to an explicit determination of the trial-wavefunction parameters, microscopic interaction and hopping parameters, ground-state energy, and explicit many-particle wavefunction, all in the correlated interacting state. For detailed discussion see main text.
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, α −1 ) in a recurrent fashion. The whole procedure is schematically illustrated by a flowchart composing Fig. 1.
The selected a single-particle basis for H 2 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 In that situation, the two-particle wave function is defined in accordance with Eq. (3), namely 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Ψ(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 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.

Covalent bonding and ionicity in H 2 case
With the help of analysis presented in Appendix A for the H 2 case, one can write down explicitly the two-electron wave function in the ground state for H 2 molecule (n = 1 case). The lowest-energy spin-singlet state is of the form where the covalent (Ψ c ) and ionic (Ψ i ) components, read Ψ c (r 1 , r 2 ) = [w 1 (r 1 )w 2 (r 2 ) + w 1 (r 2 )w 2 (r 1 )] χ ↑ (r 1 )χ ↓ (r 2 ) − χ ↓ (r 1 )χ ↑ (r 2 ) , Ψ i (r 1 , r 2 ) = [w 1 (r 1 )w 1 (r 2 ) + w 2 (r 1 )w 2 (r 2 )] χ ↑ (r 1 )χ ↓ (r 2 ) − χ ↓ (r 1 )χ ↑ (r 2 ) , with 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 Eq. (10) and (11). Explicitly, the covalency and ionicity coefficients (factors) are defined as Figure 2. The H 2 binding energy versus relative interatomic distance, calculated within EDABI and compared with restricted Hartree-Fock (RHF) and full configuration interaction (full CI) method. a 0 is the Bohr radius. EDABI yields lightly lower energy than full CI calculation at very small interatomic distance, R; this difference does not alter the main point of our qualitative discussion. and 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,17,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 α −1 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 Fig. 2. 3. Two-particle covalency vs. corresponding ionicity for H 2 molecule, calculated within EDABI method and compared with the results of Ref. [19]. Shaded regime marks a gradual evolution towards atomicity, as determined from the Mott-Hubbard criterion (see in the main text). Vertical dotted line marks equilibrium interatomic distance, whereas the horizontal dotted lines illustrate the dominant character of the covalency in that state (with the ratio r = 1.43 ∼ 2 : 1)).
In Fig. 2, we display the H 2 binding energy and have compared our EDABI calculated value with the results of configuration interaction (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 Ry is the energy of 1s state in atomic hydrogen. Next, we define the bonding and ionicity as the corresponding ratios of coefficients in Eq. (9), cf. Fig. 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 towards 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.25Ry then). The regime of strongcorrelations (Mottness) is marked explicitly in Figs. 3 and 4. It specifies a gradual evolution towards 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γ c → 0 and atomicity γ a → 1 as R → ∞ (see the discussion below).

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 i j | 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 2 , 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 Fig. 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 0 , i.e., well above the R bond 1.43a 0 . Note also that even for equilibrium distance R = R bond the hopping/interaction ratio is about ∼ 0.5, i.e., the electrons are moderately correlated.
To complete the picture, we have also plotted in Fig. 6 the antiferromagnetic kineticexchange integral J kex = 4(t + V) 2 /(U − K) versus the direct (Heisenberg) ferromagnetic value J H , both as a function of relative distance R/a 0 . The situation is that J kex > J H for any distance R and this is the reason for the spin-singlet configuration of H 2 in the ground state. In brief, electrons hoping ("resonating") between the sites, possible only in the total spin-singlet state |λ 5 , contribute essentially to the bonding.
To verify the conceptual validity of the introduced Mott threshold for atomicity onset, we have plotted in Fig. 7 the Slater-orbital size α −1 as a function of R. Upon crossing the threshold R Mott , α −1 indeed approaches rapidly with the further increasing R the 1s atomic size value a 0 = 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 Å.
In conclusion, the dominant covalent character of H 2 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 0 . The ground state energy evolves slowly, but steadily towards, the atomic-limit value. Note also that the Hartree-Fock analysis (cf. Fig. 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 2 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.

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 Fig. 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. 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 and its equivalent second-quantized from involving site occupancies 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] ν Explicitly, we propose to decompose single-occupancy probability ν in the following manner where a is the atomicity, and c is called the resonant (true) covalency, and d 2 ≡ n i↑ni↓ 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. Fig.  8a). In the strong-correlation limit (r > 1), it can be defined as c ≡ [|t − V|/(U − K)] 2 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 In the equilibrium state of H 2 , 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 Fig. 9 we provide another characteristic containing atomicity, namely the R dependence d 2 . This double occupancy probability can also characterize the ionicity. The last formula shows that the atomicity is complete when d 2 = 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 2 . 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., true covalency = 1/Motness or c ≡ 1/4r.
In conclusion, based on our analysis of H 2 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 2 . Additionally, in this way the redefined covalency is complementary to the Mottness and vice versa.

LiH and HeH + cases
We now apply the concepts introduced above for H 2 molecule to the cases of LiH and HeH + . In Figs. 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 2 , depicted in Fig. 2. The quantitative factors are different though and, in particular, the bond length is slightly larger than that in H 2 case. To characterize further those two cases we have plotted in Figs. 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.9 H −0.9 . In HeH + , molecular ion the bonding is largely covalent due to the fact that both two 1s 2 He electrons hop (mix) with the H + state with no electrons in the corresponding 1s state.
One specific feature of those two systems should be noted, which is illustrated in Figs. 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. Figure 10. The HeH + binding energy versus relative interatomic distance, obtained using EDABI method and compared with restricted Hartree-Fock (RHF) and full configuration interaction (full CI) approach. a 0 is the Bohr radius. Figure 11. The LiH binding energy versus relative interatomic distance, obtained using EDABI method and compared with restricted Hartree-Fock (RHF) and full configuration interaction (full CI) approach. a 0 is the Bohr radius. Figure 12. Many-body covalency vs. many-body ionicity for HeH + molecule. The behavior is quite similar to that for H 2 molecule (cf. Fig 3). Figure 13. Many-body covalency vs. many-body ionicity for LiH molecule. The points are taken from [19] for comparison. The covalency shows the same type of the unphysical R-dependence as in the case of H 2 . Figure 14. Atomic orbital size for He 1s and H 1s orbitals in HeH + molecular ion. Note that the Mott-type boundary has been drawn for the 1s states of He as this reached first upon increasing R. Figure 15. Atomic orbital size for H 1s and Li 2s orbitals in LiH molecule.Note the strong renormalization of atomic-orbital sizes, as well as a rapid convergence to the atomic values above R = R bond , the latter being for the ionic bonding.  [28] In Table 1 we display the binding energies of the molecules H 2 , 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 2 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]. Figure 16. Schematic representation of the molecular orbitals for (isosurface probability density cut = 0.02) on the right, expressing relative covalency and ionicity contributions on the left.

Overall properties
We now compare results for those three model systems qualitatively. First, in Fig. 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 2 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 Fig. 13. One sees that the covalency in HeH + is larger than that for H 2 molecule, a rather unexpected intuitively result. 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 0 ) is provided. This distance should be compared with the bond length R bond calculated (cf. Table 3) according to three independent methods: EDABI, full configuration-interaction (CI), and restricted Hartree-Fock (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.  [30] 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 Coulombinteraction parameters will be reduced by the dielectric constant factor if system under

Outlook
The reason for selecting the three systems analyzed here is caused by the circumstance that HeH + is strongly covalent, LiH strongly ionic, and H 2 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,31,32]. However, the antiferromagnetic kinetic exchange survives even when the states are becoming orbitally distinguishable [33].
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.

Acknowledgment
This work was supported by Grants OPUS No. UMO-2018/29/B/ST3/02646 and No. 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.

Appendix A. Eigenvalues and eigenstates for H 2 molecule
Starting from the orthogonalized restricted basis, we define the field operators aŝ or, in compact notation, asψ In the above,â iσ andâ † iσ are electron annihilation and creation operators in the state w iσ (r) ≡ w 1 (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. Eq. (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 4 2 = 6 following states The first three are the spin-triplet states with S z = +1, −1, 0, whereas the next three are interand intra-site singlets, respectively. The triplet state does not hybridize with other states and provides three 1 × 1 irreducible blocks with eigenvalues λ 1 = λ 2 = λ 3 = 1 + 2 + K − J H . The remaining three singlet states compose the 3 × 3 block, so the Hamiltonian in that Fock subspace takes the formĤ where ≡ ( 1 + 2 )/2 and U ≡ (U 1 + U 2 )/2. This formulation allows to apply this formalism to both H 2 (where U 1 = U 2 = U and 1 + 2 = ), and to HeH + and LiH, where those simplifications are not met due to inequivalent atoms involved. In the case of H 2 , the eigenvalues of Eq. (A.5) take the form 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 |λ 5 is the ground state as the λ 5 eigenvalue is the lowest one. In the limit U |t + V|, the λ 5 eigenvalue reads The last term on the right-hand side of Eq. (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 1 2 and U 1 U 2 . The singlet state |λ 5 is elaborated further throughout the main text.

Appendix B. Tables of relevant quantities and parameters for considered systems
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 2 , HeH + , and LiH. 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-31G basis set to represent the Slater functions. Numerical accuracy for the EDABI calculations is