Exact Analytical Form of Diatomic Molecular Orbitals.

We provide the exact analytical form of diatomic molecular orbitals, as given by the solutions of a single-electron diatomic molecule with arbitrary nuclear charges, using our recently developed method for solving Schrödinger equations. We claim that the best representation of the wave function is a factorized form including a power prefactor, an exponentially decaying term, a modulator function on the exponential, and additional factors accounting for nodal surfaces and the magnetic quantum number. Applying our method, we have identified unexpected extreme points along the potential energy curves, hence revealing the limitations of the well-known concepts of bonding and antibonding.

Introduction

Since the concept of orbital was proposed in the early days of quantum chemistry, this concept has been widely accepted by chemists as the most powerful theoretical tool to gain deep insight into chemical problems. Although orbitals are based on the independent single-electron picture, which is only an approximation to the correlated many-electron picture described by many-electron wave functions in atoms and molecules, their concise, intuitive, and visualizable features make them the most commonly used terminology in chemists’ routine discussions.

Among various orbitals introduced for different purposes, the most successful ones are undoubtedly the atomic orbitals (AOs) of hydrogenic systems. Rigorously defined as the solutions of a single-electron Schrödinger equation (SE) with spherical symmetry, each of these AOs can be analytically written as a product of a radial function and a spherical harmonic. For many-electron atoms with electron–electron interactions, although one cannot separate variables to deduce single-electron SEs, the clever idea of introducing electronic screening and effective nuclear charge has allowed a reduction to the hydrogenic picture where electrons fill into different energy levels.[1,2] The resulting AOs for many-electron atoms capture the essential physics and have achieved enormous success in explaining the electronic structure of elements in the periodic table.

Yet, chemistry deals with molecules. Apart from understanding AOs, perhaps even more important for chemists is to decipher the molecular orbitals (MOs). In contrast to atoms, molecules display much greater complexity in the presence of multiple electrons and nuclei. The idea of treating electron–electron interactions as an effective screening as used in atoms is not feasible for molecules due to the ambiguity of assigning effective nuclear charges. It turns out that a plausible way is to invoke a fictitious noninteracting system with an effective potential. The resulting MOs can be generated by solving the SE of that particular potential, either local or nonlocal, for example, as has been practiced by the Kohn–Sham density functional theory (KS-DFT)[3] or the Hartree–Fock (HF) theory, giving rise to different MOs associated with different methods or different functional approximations in DFT. This ambiguity of MOs can be once again attributed to the attempt at using approximate orbitals to describe a correlated many-electron wave function.

Nevertheless, for single-electron molecules, MOs can be unambiguously defined. For diatomic molecules, in particular, these MOs are named σ, π, δ, etc., in analogy with the s, p, and d types for AOs, and have given tremendous inspiration to chemists. In fact, they have become the most essential part of modern chemistry textbooks.[1,2,4] In contrast to the s, p, and d orbitals whose analytical forms have been clarified, however, the exact analytical forms of MOs remain a challenge. For example, the MOs for the simplest molecular ion, H2+, have been studied intensively, but their compact analytical forms are still elusive.[5−19] Instead, these MOs are usually numerically represented as linear combinations of atomic orbitals (LCAOs).[20,21]

In this paper, we show that we can do better than basis expansion by finding the exact analytical forms of all the MOs. In particular, we apply our recently developed method for solving SEs to single-electron diatomic molecules and derive the exact expressions for σ, π, δ orbitals and so on. For each MO, our formula is in an exactly factorized form, i.e., casting the wave function into a product of multiple factors resembling the exact formula of an AO. This is in contrast to the conventional LCAOs and other basis expansion methods that decompose the wave function into a sum of infinite terms. Our representation of MOs reveals their intrinsic analytical structure and furthermore has proven to have computational advantages, as shown in our recent work.[22] Importantly, the newly obtained analytical form of MOs could give us new insight into the nature of chemical bonds.

Methods

We start by reviewing our newly proposed method for solving one-dimensional (1D) SEs, as has been implemented for finding the exact analytical solutions of 1D hydrogen atom and H2+ with soft Coulomb potentials.[23,24] In particular, for a nicely behaved potential in 1D, by formulating the corresponding ground state wave function as ψ = Ceβ, we can transform the SE for ψ(x) into a Riccati equation[25] for , where energy enters as a parameter. In doing this, we reduce the second-order ordinary differential equation (ODE) into first-order at the sacrifice of linearity. The equation is then solved by expanding u into a Taylor series, which combined with the boundary conditions ultimately leads to an algebraic equation that determines the energy.[23,24] In this work, we extend our approach to real-world molecules in 3D. As will be shown, the increased dimensionality leads to an increase in the number of algebraic equations and unknown variables, yet the basic structure remains similar.

Without loss of generality, let us consider the generic single-electron diatomic molecular problem in atomic units:Here R’s are the electron–nuclear distances, R is the nuclear separation, and Z’s are the nuclear charges (we assume Z1 ≥ Z2). Although eq appears to be a coupled equation in terms of the three Cartesian coordinates, it has been shown to be separable in spheroidal coordinates (also called confocal elliptic coordinates).[4−7] In particular, , , and ϕ is the angle of rotation of the electron about the z axis; see Figure for an illustration. ξ and η are analogous to the radial distance r and the cosine of the polar angle θ in the spherical coordinate system; in the limit R → 0, ξ → 2r/R and η → cos θ. The wave function can be then factorized into ψ = M̃(η)Ñ(ξ)e. Here, analogous to the hydrogen atom, the equation for ϕ is an eigenvalue equation, for which one can define a magnetic quantum number m = 0, ±1, ±2, ... characterizing the z component of the electronic orbital angular momentum. Differing from the hydrogen atom, here ξ and η do not obey eigenvalue equations; instead, they satisfy two decoupled ODEs with two undetermined separation constants:Here Ẽ ≡ E – Z1Z2/R is the electronic energy, and à is related to the total angular momentum and the Runge–Lenz vector.[26−29] Similar to eigenvalue equations, though, only when Ẽ and à take special values can eqs –3 have solutions. It is worth noticing that these ODEs have singularities at η = ±1 and ξ = 1, respectively, indicating that M̃(η) and Ñ(ξ) cannot be finitely differentiable at these boundary points. In light of this and inspired by previous works,[5,8,9] we make the following change of variables, M̃(η) = (1 – η2)|M(η) and Ñ(ξ) = (ξ2 – 1)|N(ξ), which after straightforward derivation gives explicit equations for M(η) and N(ξ):Here A = à + m(m + 1).

Figure 1

Illustration of the spheroidal coordinates for a diatomic molecule with arbitrary nuclear charges. Here, by definition, ξ ≥ 1, −1 ≤ η ≤ 1, and 0 ≤ ϕ ≤ 2π.

To solve for the ground state, we apply our previous technique to eqs –5. As with the ground state of the hydrogen atom, here we assume that the wave function has an absence of nodal points and m = 0.[30] Then writing the wave function in exponential forms, M(η) = C1eβ and N(ξ) = C2eβ, and denoting and , we deduce the following Riccati equations:Equations –7 will be solved by performing Taylor expansions. To have a finite radius of convergence, it is preferable to transform ξ onto a finite interval [0,1) by introducing . The resulting ODE for q readsWe then expand u(η) and v(q) into Taylor series:Here we assume that the validity of these expansions extends to the entire domain of η and q, i.e., the radius of convergence is greater than 1. This assumption is valid for most cases of interest but has exceptions, which will be discussed later. Comparing terms η and q order by order in eqs and 8, we arrive at recursive relations for u’s and v’s, so that each u and v can be represented as a function of u0, v0, Ẽ, and A.

To determine these unknown variables, one has to invoke the boundary conditions. Assuming that the derivative terms in the ODEs are finite, one can readily find that they are eliminated at the boundary points η = ±1 and q = 0, 1; hence, eqs and 8 reduce to four algebraic equations regarding u(±1), v(0), and v(1), which by eqs –10 can be further rewritten in terms of u’s and v’s. Therefore, we end up with four coupled algebraic equations for four unknowns, which can be solved by a multidimensional Newton’s iteration approach.[31,32] Here we note that physics guarantees the existence and uniqueness of the solution.

Once u0, v0, Ẽ, and A are obtained, all the u’s and v’s are readily accessible by the recursive relations. It follows that β1 and β2 can be obtained by explicit integration:Here in eq 11, α1 = −u0 > 0, and . For heteronuclear cases, β1 is dominated by −α1η for large R; for homonuclear cases, α1 is strictly zero, and β1 reduces to F1 for all R. In eq 12, α2 = −∑∞v, γ = −∑∞kv, and F2(ξ) = ∑∞vh(ξ–1) with . Using the boundary condition at ξ = 1, we can prove and , whose limiting behaviors for R → 0 (united atom) and R → ∞ have been analyzed accordingly.[32]

The ground state wave function is then given bywhere C = C1C2 is the normalization constant and F(η, ξ) = F1(η) + F2(ξ). Cast in an exactly factorized form, the wave function manifests its analytical structure in the most concise and informative manner, which is much more physically meaningful than an LCAO type of basis representation. Compared with the hydrogen atom, we recognize both familiar features and new structures. The primary similarity is in the exponential decay, where the analogous decay pattern in eq is through e–α and the rate of decay is closely related to the energy. The major difference appears in the additional factors, among which we call special attention to F(η, ξ), which we define as the modulator function in the sense that it modulates the exponential decay. In contrast to the other terms in eq , F(η, ξ) can only be written as a series expression. Yet, it can have qualitatively different behavior for small and large nuclear separation and deserves some further discussion.

In particular, in deriving eq , our assumption about the Taylor expandability of u(η) and v(q) implies that they are free from singularities within the unit circle on the complex plane. It turns out, however, that the assumption for u(η) is violated for large R when Z1 ≈ Z2. Nevertheless, eq can still hold in such cases if one modifies the definition of F1. This can be achieved by moving the singularities outside the unit circle through a Mobius transformation; the resulting formula of F1 can be found in the Supporting Information.[32]

Importantly, the new analytical structures identified in eq appear to be generally applicable to the ground solution of SEs for Coulomb systems, as they have been observed in our previous works on 1D model problems.[23,24] Moreover, the exact formula also sheds light on simple approximate formulas. In fact, one can achieve a high accuracy by approximating F(η, ξ) as simple elementary functions and parametrizing the variational wave function with as few as three parameters.[32]

Our method for finding the exact ground state wave function can be extended to target all the excited states, for which one shall additionally factorize the nodal points, i.e., M(η) = C1∏(η – a)eβ, N(ξ) = C2∏(ξ – b)eβ. Here L and K are the number of nodes in M(η) and N(ξ), respectively; a and b denote the nodal positions. This is in the same spirit as our previous work on 1D model problems.[22] Substituting the factorized formulas into eqs –5, we can derive analogous equations to 6 and 8, which again can be solved by Taylor expansion. It is worth mentioning that the K + L nodal points now become unknowns. Substituting each of these variables into the ODEs leads to a new algebraic equation, which along with the boundary conditions is sufficient to determine all the unknowns.[32]

After repeating essentially the same steps for finding the ground state, we can deduce the exact formula for a generic eigenstate characterized by quantum numbers KLm as the following:Here and . The formulas of F and α1 are formally the same as those of the ground state, although u’s and v’s take different values. One can readily see that eq reduces to eq for the ground state (m = 0 without nodal points).

Equations –14 are the key results of this paper. In contrast to the wave function formulas proposed in the literature that involve an infinite summation as part of the factorization,[8,9]eq is an exact and complete factorization, elucidating the analytical structure as much as possible. For example, it clearly shows that the nodal surfaces of MOs are hyperboloids and ellipsoids with the nuclei as foci, validating the argument in the literature[6,7] by specifying the exact positions. As an additional remark, there have also been attempts at solving the H2+ problem by transforming the SE into Riccati equations, such as the Riccati–Padé method (RPM).[33−36] Yet, RPM uses the Padé approximation to represent u and v rather than targeting their exact formulas or the exact wave function. In addition, RPM gives no knowledge of the nodal points.

Results and Discussion

Next, we demonstrate the usefulness of our exact formulas by showing some intuitive examples. Starting with the H2+ problem, we compute potential energy curves of representative low-lying states; see Figure . Here, instead of labeling states with (KLm), we use the united-atom designation,[9,10] with (nlm) based on the types of reduced AOs in the limit R → 0 . One can easily work out the relations n = K + L + m + 1 and l = L + m. As with AOs, s, p, d, ... are used to reflect the information on l; by contrast, σ, π, δ, ... are used to specify the value of m (corresponding to m = 0, 1, 2, ...), which also shows the types of MOs.

Figure 2

Energy curves of representative states of H2+, labeled by united atom designation. Curves in the energy range from −0.18 to −0.04 Hartree are shown in the main plot, while the lowest bonding and antibonding states are shown in the lower inset. Bonding/antibonding states are drawn in solid/dashed lines. Interestingly, the 4fσ state has a local minimum around 21 Bohr; and the 4dσ state has a local maximum around 44 Bohr (see the enlarged plot in the upper inset).

Chemists are used to characterizing MOs as bonding or antibonding by judging (i) whether a buildup of charge occurs at the bond midpoint[1,2,4] and (ii) whether the attractive forces between atoms are strengthened or weakened by occupying the MO.[37−39] For H2+, by (i) an MO with a mirror/nodal plane between the nuclei is a bonding/antibonding orbital; by (ii), a bonding state is supposed to produce an energy well, while an antibonding state shall monotonically decrease its energy upon dissociation.[4] (i) and (ii) are consistently true for some lowest-lying states, such as the 1sσ and 3dσ states. However, by computing the exact energy curves to an extended range and to a high precision, we find that this is not generally true for other states. For example, Figure shows that the 4fσ state, which is an antibonding state by (i), develops a local minimum at R ≈ 21 Bohr, while the 2sσ state, which is a bonding state by (i), is monotonically decreasing. Even more surprisingly, bonding states such as 4dσ can have a local maximum in the large-R region, corresponding to a transition state. In fact, such examples are widespread as shown in Table , where we have tabulated the positions and energies of extreme points along the energy curves of some low-lying states.

Figure 3

(a) Energy curves of the 1sσ state for different Z1, fixing Z2 = 1. Energy at the dissociation limit of each curve has been set to zero. (b) Binding energy (energy difference between the dissociation limit and the minimum, if it exists) as a function of Z1 and Z2, could be negative. For the blank area, there exists no minimum in the energy curve.

Table 1

Extreme Points along Energy Curves of Representative Low-Lying States of H2+, Including Minima and Maxima (If They Exist) Denoted as Rmin and Rmax, Respectivelya

state	Rmin	Emin – ED	Rmax	Emax – ED
1sσ	2.0	–1.03 × 10–1
2pσ	12.5	–6.08 × 10–5
3dσ	8.8	–5.00 × 10–2
4fσ	20.9	–5.66 × 10–3
4dσ	17.8	–3.26 × 10–3	43.5	1.65 × 10–4
5gσ	23.9	–2.27 × 10–3
2pπ	7.9	–9.51 × 10–3	25.8	1.44 × 10–4

Their energies, Emin and Emax, are shown relative to their respective dissociation limit ED. Some states, such as 2sσ, 3pσ, and 3sσ, have neither maximum nor minimum and are not listed in the table. All values are in atomic units. Here we note that the shallow minimum of 2pσ has also been reported in ref (4).

From the LCAO perspective, a bonding/antibonding orbital has also been associated with a decrease/increase of energy relative to the separated atoms.[1,2] Yet, we find that this is not consistent with (i) either. For example, the 4fσ state at finite R, which is antibonding by (i), has an energy lower than its dissociation limit; see Figure .

For heteronuclear molecules, binding curves become more sophisticated. In Figure a, we compare ground state energy curves with different Z1, fixing Z2 = 1. As Z1 increases, we see a weakened binding behavior. When Z1 reaches 1.2, an obvious transition state emerges in the energy curve, separating the local minimum from the energetically more favorable dissociation limit. This minimum fades away when further increasing Z1. If we allow Z2 to change, all possible combinations of Z1 and Z2 that lead to a minimum (binding) encircle the colored region in Figure b.

These observations thus call into question the validity of traditional bonding vs antibonding concepts. Furthermore, as the interplay between the electronic energy Ẽ and the nuclear Coulomb repulsion yields energy curves with so many complicated features for systems as simple as single-electron diatomic molecules, one could likely find unexpected intermediates or transition states on the energy surfaces (particularly for excited states) of other molecules.

Besides accurate energies, perhaps more importantly, our formulas can accurately describe important features of MOs at a much lower computational cost. For example, eq gives accurate hyperboloids (3dσ) and ellipsoids (2sσ) as nodal surfaces as shown in Figure a and b. By contrast, conventional basis expansion methods cannot capture the nodal shapes even qualitatively with commonly used basis sets; see Figure c and d. More demonstrations of the computational advantage of our method over conventional basis expansion can be found in the Supporting Information. Importantly, if we decompose MOs using LCAOs in the infinite separation limit, the resulting AOs are sp hybrids rather than pure 2s or 2p orbitals.[32]

Figure 4

Contour plots of wave functions with a nuclear separation of 14 Bohr: (a) 2sσ state by our method; (b) 3dσ state by our method; (c) 2sσ state by basis expansion; (d) 3dσ state by basis expansion. For (a) and (b), 80 nonzero u’s and 64 v’s are used. For (c) and (d), the basis set of aug-cc-pV5Z (160 basis functions) is used. Apparently, the basis expansion here gives qualitatively incorrect nodal surfaces.

The modulator function F is an essential term in our factorization. In Figure , we show that F for the ground state behaves qualitatively differently for small and large R, for homonuclear as well as heteronuclear cases. In the limit R → 0, in particular, F approaches a constant because the wave function reduces to the ground state of a hydrogenic atom, given by Ce–α. This is manifested in Figure a and c, where the overall scale is small. When R is large, however, F changes rapidly in the internuclear region; see Figure b and d. This shows that F can be used as an indicator for distinguishing compact from dissociated molecules, which is consistent with our observations for H2+ in 1D,[24] and could be useful for tackling the delocalization error in density functional theories, particularly for improving the recently developed localized orbital scaling correction (LOSC) functional.[40−45] Of particular interest is the behavior of F for heteronuclear cases such as HeH2+. For small R as in Figure c, F is smooth and delocalized over the two nuclei, resembling the united atom limit. For large R as in Figure d, interestingly, we find that F is localized around the lighter atom (hydrogen), although the wave function is localized near the heavier atom (helium).

Figure 5

Modulator F shown as a function of x and z for the ground state of H2+ (upper panels) and HeH2+ (lower panels). In (a) and (c), R = 0.5 Bohr; in (b) and (d), R = 5 Bohr. Here we present results in the circled rather than the squared region.

Conclusions

In this paper, we have obtained the exact analytical forms of diatomic MOs, as given by the solutions of a single-electron SE for a diatomic molecule. We show that the best way of representing the ground MO is in our factorized form in eq involving a power prefactor, an exponentially decaying term, and a modulator on the exponential, while the best way of representing excited state MOs involves additional factors accounting for the nodal surfaces and the magnetic quantum number as in eq . Our factorized formulas are formally intuitive and physically informative and unify the exact formulas of AOs and MOs. The usefulness of our exact formulas has been demonstrated in several aspects. Importantly, our new findings about bonding/antibonding MOs have revealed the limitation of these concepts particularly in stretched molecules.