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.

I. Some details of transforming Schrödinger equation to algebraic equations S1 II. Some details of the Riccati equations for excited states S3

III. Some details of Mobius transformation and determination of coefficients S4
IV. Details of few-parameter ground state variational solutions S6 V. Some supplemental results S8 A. Evolution of MOs and their nodal surfaces as bond is stretched S8 B. Computational advantage of our method over LCAO expansion S9 C. Some additional analysis on the modulator function S10 References S14

I. SOME DETAILS OF TRANSFORMING SCHRÖDINGER EQUATION TO AL-GEBRAIC EQUATIONS
In the main text, we have shown that the Schrödinger equation for ψ can be transformed into the following equations for M and N , For the ground state, M (η) = C 1 e β 1 (η) , N (ξ) = C 2 e β 2 (ξ) and m = 0. The resulting equations for β 1 and β 2 read the Riccati equations reduce a 2nd-order ODE to two 1st-oder ODEs at the sacrifice of linearity of the equation. Furthermore, by introducing q = 1 − 1 ξ , the Riccati equations can be cast into the following forms: Here One can see that f i 's and g i 's are polynomials of η and q, i.e., f i = j f j i η j and g i = j g j i q j . With the Taylor expansion of u(η) = ∞ k=0 u k η k and v(q) = ∞ k=0 v k q k , one can compare the η k−1 (k ≥ 1) terms in Eq. (S5) and q k terms in Eq. (S6), which lead to In deriving Eq. (S8), we have used g 0 1 = g 0 2 = 0. Solving for u k and v k in Eq. (S7) and Eq. (S8) leads to the following recursive relations: S2 by which each u k or v k becomes a function of four unknowns, u 0 , v 0 ,Ẽ, A. To determine these unknowns, we invoke the boundary conditions at η = ±1, q = 0 and 1: Here we have been assumed that u, v and their derivatives are finite. By Eq. (S13), we can , reducing the number of unknowns to 3.
Then the remaining three equations can be written as correlated algebraic equations of u 0 , E and A once we plug in the Taylor expansion formula of u and v at the boundary points.
Truncating the Taylor series to finite order, i.e., considering only finite number of u k 's and v l 's, Eqs. (S11), (S12) and (S14) can be solved by the standard multi-dimensional Newton's iteration approach. [1] Although with a given truncation length, there could exist multiple solutions, one can find a unique convergent series when enlarging the truncation length, and the exact A,Ẽ and u 0 are given by this limiting solution.

II. SOME DETAILS OF THE RICCATI EQUATIONS FOR EXCITED STATES
For a generic state, the Riccati equations are given by Assuming finite derivatives at the nodes of the wave functions, one can write M (η) = P (η)e β 1 (η) and N (ξ) = Q(ξ)e β 2 (ξ) , where P (η) = L l=1 (η − a l ) and Q(ξ) = K k=1 (ξ − b k ). Here L and K are the number of nodal points of M (η) and N (ξ), respectively. Substituting these exponential expressions into Eq. (S15) and Eq. (S16) and after simplification, we arrive Here we have denoted u = dβ 1 dη , v = dβ 2 dξ . Next, we introduce q = 1 − 1 ξ in Eq. (S18) to transform the domain of ξ into a finite interval [0, 1). The resulting equation reads Eq. (S17) and Eq. (S19) can be solved analogously as for the ground state; here one shall treat each node (a l or b k ) as an additional variable, and introduce a nodal condition accordingly by plugging η = a l or ξ = b k into the formulas.

III. SOME DETAILS OF MOBIUS TRANSFORMATION AND DETERMINA-TION OF COEFFICIENTS
In some cases u(η) has singularities in the unit circle on the complex plane so that our Taylor expansion is no longer valid on the entire domain of [0, 1). Then we shall move these singularities outside the unit circle through a Mobius transformation as has been practiced in Ref. [2]. Since this transformation is valid only for variables defined on [0, 1], it cannot be applied directly to our η, which is defined on [−1, 1]. Then our solution is to decompose u(η) into an even and an odd function of η, u(η) = s(p) + 2ηw(p) with p = η 2 . The resulting ODEs for s(p) and w(p), deducible from Eq. (S17), become coupled equations. Take the ground state for an example, they read 2(1 − p) ds dp Eqs. (S20) and (S21) are equivalent to Eq. (S5). By choosing proper parameter t > 0, one can move the singularity of u out of the unit circle by the following Mobius transformation, Then performing Taylor expansion in terms of p , we have convergent series s = ∞ k=0 s k (p ) k and w = ∞ k=0 w k (p ) k in the unit circle. Transforming p back into η, we can obtain the expression of u(η). Then by integration, We can further separate a linear term from β 1 , i.e., β 1 (η) = −α 1 η + F 1 (η). Then α 1 = −u(η)| η=0 = −u 0 = −s 0 , and Since β 1 and α 1 are t-invariant, F 1 must also be t-invariant even though its definition formally involves t.
For some excited states, singularities of v(q) can also arise in the unit circle of the q plane when R is small. Then we can introduce another Mobius transformation for q, where d ∈ [0, 1) is a parameter. Expanding v into a Taylor series of q leads to a convergent series in the unit circle, v = ∞ k=0 v k (q ) k . Then transforming q back into ξ and integrating S5 over ξ, we obtain the formula for β 2 as To recast Eq. (S26) as . α 2 and γ can be further simplified using the following relations, Setting q = 1 in Eq. (S19) and the derivative of Eq. (S19), we have by which we can solve α 2 = − 1 2Ẽ R 2 and γ = R(Z 1 +Z 2 ) 2α 2 − (K + m + 1). α 2 and γ are apparently d-invariant. Then it follows that F 2 is also d-invariant although it formally involves d. In the special choice of d = 0, one recovers the expression in the main text. Table S1 summarizesẼ, α 1 , α 2 and γ for the ground state in the limit R → 0 and R → ∞.

IV. DETAILS OF FEW-PARAMETER GROUND STATE VARIATIONAL SOLU-TIONS
The exact ground state wave function is given by a product of a power prefactor, an exponentially decaying term modulated by an F function. Inspired by this form, a simple S6 Table S1: Limiting behaviors ofẼ, α 1 , α 2 and γ for the ground state. When Z 1 = Z 2 , α 1 is strictly 0 for any R, which is not shown in the table. Term trial wave function for variational approximation is to neglect the modulator, where the subscript "V 3" indicates the number of variational parameters. The resulting estimate for the energy is given by the minimizer of the variational formula E = ψ|Ĥ|ψ ψ|ψ with respect to λ i 's. Improvement over this energy can be achieved by approximating the modulator by simple elementary functions that leads to the correct asymptotic behavior, giving rise to 4-parameter trial wave functions: In particular, the modulator has been approximated by ln cosh(λ 4 η) and λ 4 η 2 , respectively.
In Figure S1, we compare variational results with the exact ground state energy curve for H + 2 and HeH 2+ . Here (and also in the main text) the exact energy curve is accurate to 10 −10 Hartree. For H + 2 , λ 3 is strictly zero by symmetry, thus reducing the number of variational parameters to 3 even though it is under the name of "V4". As shown in Figure S1(a)(b), scheme "V 4" performs the best among the schemes tested, whose error is on the magnitude of 10 −4 Hartree, while scheme "V 3" gives a qualitatively wrong energy curve because of the unphysical plateau in the wave function. This shows the necessity of the modulator function. For HeH 2+ , to demonstrate the necessity of splitting the term −α 1 η from β 1 , we also consider the following 2-parameter trial wave function, As shown in Figure S1(c)(d), here scheme "V 3" is sufficient to yield an accurate enough energy curve, while scheme "V 2" gives qualitatively wrong energy because of an unphysical plateau in the wave function. In summary, one can construct trial wave functions with as few as 3 parameteres that give accurate ground state energy and wave function. 2 (ψ 2s,He + ψ 2pz,He ) with ψ 2s,He and ψ 2pz,He being the atomic orbitals of He. We have set the z-axis as the bonding axis and placed the helium nucleus at z = 0. The nodal surfaces have been highlighted by the red lines. As R increases, the nodal surface (which is an ellipsoid) gradually enlarges its radius of curvature near the helium position and expands its encircled region until becoming a paraboloid when R = ∞, which is nothing but the nodal surface of the sp hybrid atomic orbital of He.

A. Evolution of MOs and their nodal surfaces as bond is stretched
for the roots of these functions, we end up with η = η 0 or ξ = ξ 0 for some η 0 and ξ 0 . By definition, such nodal surfaces are hyperboloids or ellipsoids. To see how such nodal surfaces evolve into the nodal surfaces of atomic orbitals in the dissociation limit, particularly for a heteronuclear diatomic molecular ion, we illustrate in Figure S2 using the example of 2sσ orbital of HeH 2+ for increasing R. It is worth noticing that the 2sσ molecular orbital reduces to the sp hybrid atomic orbital rather than pure 2s or 2p orbital of He. This is reasonable because 2s and 2p orbitals of He become degenerate in the infinite nuclear separation. Such hybridization also occurs for homonuclear molecules such as H + 2 , see Figure 4 of the main text.

B. Computational advantage of our method over LCAO expansion
We have shown in the main text that the conventional basis expansion methods cannot capture the nodal surface shapes even qualitatively with commonly used basis sets. Here we further demonstrate that our method can achieve greater accuracy with much smaller  Figure S3: Comparing computational efficiency between conventional basis expansion (computed with pySCF) and our method, as indicated by their energy error for the ground state of (a) H + 2 and (b) HeH 2+ for a chemically interesting range of R. Here E exact is obtained using our method with extremely large number of u k 's and v k 's. For each line presented in the figure, we have marked the total number of basis functions (for basis expansion) or terms used in the series summation (for our method). For H + 2 , half of the u k 's vanish due to even symmetry and only those non-vanishing u k 's are counted.
computational cost. As shown in Figure S3(a), for the ground state of H + 2 with bond length of chemical interest (R < 8 Bohr), our solution gives very accurate energy (for most cases error less than 10 −7 Hartree) with a total number of 60 u k 's and v k 's [3]. This is several orders of magnitude better than the aug-cc-pV5Z results with 160 basis functions. In fact, one has to use an unprecedentedly large Gaussian basis set to reach our accuracy, as one can estimate from the error trend when enlarging basis set. Similar situation is also observed for heteronuclear molecules such as HeH 2+ , see Figure S3(b). Furthermore, we have found that our computational advantage is even more pronounced for excited states. For example, for the 2sσ state of H + 2 , we can achieve much better accuracy than Gaussian basis set with as few as 30 u k 's and v k 's, as shown in Figure S4.

C. Some additional analysis on the modulator function
The number of u k 's or v k 's (let us call them cut-offs here) needed for an accurate evaluation of the wave function and modulator function are determined by converging the total energy to a given threshold (set to 10 −10 Hartree in the present work). For chemically interesting R (typically within 10 Bohr for the ground state, could be larger for excited states), a cut-off of less than 50 is sufficient in most cases. For example, a cut-off of 20 is sufficient to give S10  Figure S4: (a) Comparing computational efficiency between conventional basis expansion (computed with pySCF) and our method, as indicated by their energy error for the excited 2sσ state of H + 2 for a chemically interesting range of R. Note that the basis expansion results have a spurious avoided crossing in the range of R = 3.5 − 4.5 Bohr where a nonadiabatic transition takes place between 2sσ and 3dσ. To make a fair comparison, here we show the diabatic curves with 2sσ symmetry for regions away from the nonadiabatic transition. (b) Comparing the exact energy curves of 2sσ and 3dσ states with the basis expansion results (using aug-cc-pV5Z as basis), which have a spurious avoided crossing. [4] accurate modulator functions for HeH 2+ , see Figure S5 for an illustration. For larger R, however, the cut-off (especially for u k ) grows with R. For example, 500 u k 's are needed for R = 50 Bohr. Yet such large R is of little significance for chemical problems.
In our previous work on 1D H + 2 with soft Coulomb potential, we have obtained the exact analytical ground state wave function for the following Schrödinger equation: Here we compare our 3D results with 1D particularly for the modulator function F . Note that F is a singlevariable function of λ = 1− √ 4a 2 +R 2 y in 1D, where y = a 2 + (x + R/2) 2 + a 2 + (x − R/2) 2 is analogous to our Rξ in 3D. [2] To make an easy comparison, here we reduce F (ξ, η) also to a single-variable function by plotting it along a line parallel to the bonding axis with distance x 0 and parametrizing this line by λ = 1 − √ 4x 2 0 +R 2 Rξ , in analogy to that in 1D.
The comparison is shown in Figure S6. As can be seen, the overall trend is similar in 1D and 3D; F is a mildly varying function of λ for small R and grows rapidly in the small-λ region for large R. This validates that our 1D model of H + 2 captures the essential 3D physics during molecular dissociation. Although F is a decreasing function in 1D but not in 3D for small R, this is largely due to the artifact of the soft Coulomb potential and only of minor S11 importance.   Figure S6: Comparing modulator F of H + 2 with different bond length R (in unit of Bohr) as a function of λ for 3D (solid lines) and 1D with soft Coulomb potential (dashed lines). For the 3D case, F is plotted along a line parallel to the bonding axis with distance x 0 = 1 Bohr; for the 1D case, the soft Coulomb parameter is a = 1 Bohr. We have shifted each curve here such that F is zero when λ = 0.