Given the remarkable stability of P 4 , it is not unreasonable to search for its heaver valence isoelectronic species. Despite many attempts, it is unclear whether the Bi 4 molecule has been positively identified in the laboratory, apart from mass spectrometry in several studies, [1][2][3][4][5] most recently those of Duncan and coworkers. 6,7 As of yet, no vibronic spectra have been unambiguously assigned to this molecule, though many have tried. [8][9][10][11][12] Although Bi 4 has proven evasive to spectroscopists, it has been explored by theoretical chemists using several varieties of density functional theory (DFT). Various characteristics of Bi 4 have been studied, including its electronic structure, [13][14][15][16] electron affinity, 14,[16][17][18] ionization potential, 14,16,18 , and binding energy, 14,19 . The preferred geometry 13,14,16,[18][19][20][21] and harmonic vibrational frequencies 14,19 of Bi 4 have been computed mostly by DFT studies, with one study 13 utilizing dynamical correlation methods. While DFT and static correlation methods have been fruitful in many endeavors, they are not the highest level of theory that may be utilized, and high precision is of the utmost importance for this system.
Several spectroscopic studies [8][9][10][11][12] initially claimed to have identified three vibronic emission bands of Bi 4 . Arrington and Morse 22 later (2008) observed these same emission bands but identified the source as Bi 3 via time-of-flight mass spectrometry. An excellent ab initio analysis of the electronic structure of Bi 4 at the MRCISD/CASSCF level of theory was provided a) Electronic mail: ccq@uga.edu by Zhang and Balasubramanian as early as 1992. 13 In the latter study, the authors give further evidence that these emission bands should not be assigned to Bi 4 on the grounds that only one of the low energy Bi 4 electronic excitations is formally dipole allowed, and its expected excitation energy is qualitatively different from the excitations observed. Additionally, Wakabayashi et al. analyzed small bismuth clusters using laser induced fluorescence spectroscopy and their results also support the reassignment of the supposed Bi 4 emission bands to Bi 3 . 23,24 Several theoretical studies have explored the energetically low-lying constitutional isomers of Bi 4 and their corresponding harmonic vibrational frequencies, although the highest level of theory used in these studies was MRCISD with a triple zeta basis set. From a Lewis-bonding perspective, 25 one would expect that the lowest energy isomer would be one that has three bonds per bismuth to find pairs for its three unpaired valence electrons. Jia et al. noted that due to the relativistic contraction of the bismuth 6s orbital, 14 covalent bonding in Bi n species will be dominated by the 6p orbitals. While the Lewis-bonding perspective yields one constitutional isomer of Bi 4 , other energetic local minima isomers are theoretically possible.
Lauher 26 lists four possible geometric isomers of metal cluster tetramers labeled as tetrahedral (T d ), butterfly (C 2v and D 2h ), and square planar (D 4h ). Yuan et al. predicted the three lowest energy isomers of Bi 4 , utilizing the BPW functional with a relativisitic effective core potential (RECP) and a double numeric plus d polarization functions (DN+d) basis set, to be, T d , C 2v (+27.4 kcal mol -1 ), and D 4h (+39.9 kcal mol -1 ). 18 Akola and co-workers mention only the T d and C 2v structures as low energy isomers in their computations using the PBE functional with an RECP and a plane-wave basis , 19 shedding some doubt on the existence of the square planar isomer. Several studies [27][28][29][30] have shown that the C 2v isomer can exist for P 4 analogues if a ligand is inserted into a bond, thus giving experimental validation for the existence of this isomer in group 5 tetramers. Gausa et al. found the lowest energy isomer of the Bi - 4 anion to be the C 2v structure, 16 begging the question of why the T d isomer is not a minimum for the anion. Jia et al. determined the lowest energy isomer of neutral Bi 4 to be the T d structure through PBE/DN+d with RECP computations. However, they computed the D 2d structure of anionic Bi - 4 to be the lowest in energy, again not finding an anionic T d local minimum. 14 Considering the longstanding misassignment of spectra previously mentioned, no isomers of Bi 4 have been spectroscopically confirmed.
Theoretical treatment of bismuth clusters comes with its own set of challenges. The considerable number of electrons to be modeled, along with the presence of relativistic effects within the immense core has resulted in many low level theoretical treatments. One method of solving both of these problems has been the application of effective core potentials (ECPs). 31 This computational tool accounts for relativistic core effects while allowing for high-level correlation treatments, such as coupled-cluster (CC) theory, of the valence electrons. In the case of bismuth, large core ECPs reduce the number of explicitly treated electrons to those occupying the valence 6s and 6p-orbitals. The Bi 4 molecule will thus be analyzed with the aid of ECPs in this study.
Optimized equilibrium geometries and harmonic vibrational frequencies for Bi 4 isomers of T d and C 2v symmetry were computed in this study. These two isomers were found to be the only local minima of those mentioned by Lauher. 26 Additionally, the energy difference between these isomers and the dissociation energy of Bi 4 (T d ) into two Bi 2 molecules were obtained, and Natural Bond Orbital 32 (NBO) analysis was performed on the T d and C 2v isomers.
The equilibrium geometries were obtained by utilizing Peterson's 33 correlation consistent quadrupleζ basis set designed for use with relativistic psuedopotentials 31 (cc-pVQZ-PP). The large core psuedopotentials used in this work were developed by Metz and coworkers 31 and encompass all electrons that are not present in the 6s or 6p orbitals of bismuth. This ECP78MDF encompasses 78 electrons per bismuth and allows for treatment of the remaining five electrons per bismuth with higher levels of theory. Coupled cluster theory with full single, double, and perturbative treatment of triple excitations 34,35 [CCSD(T)] was utilized in the geometry optimizations with the Cfour 2.0 36 (henceforth Cfour) suite of electronic structure codes. The T d and C 2v isomer geometry optimizations had convergence criteria for the RMS energy gradient of 10 -8 . Molpro 2010 37 was utilized to check for multireference character in the wavefunction space of each isomer through its T 1 38 and D 1 39 diagnostics. NBO 32 analysis was performed on both isomers utilizing the Q-Chem 5.0 40 quantum chemistry software package interfaced with NBO 6.0 41 augmented by Natural Resonance Theory (NRT). 42 Harmonic vibrational frequencies were obtained using a 5-point finite difference method also utilizing Cfour at the CCSD(T)/cc-pVQZ-PP level of theory as well, and with the same convergence criteria as the geometry optimizations. After the harmonic frequencies were obtained they were inspected to ensure that no imaginary vibrational modes were present and that each optimized geometry represented a genuine minimum on the potential energy surface (PES). Full cubic and semidiagonal quartic anharmonic corrections to the harmonic vibrational frequencies were computed using second order vibrational perturbation theory (VPT2) 43,44 through finite differences, as implemented within Cfour and at the same level of theory as that for the harmonic vibrational frequencies.
Focal point analysis [45][46][47][48] was employed in determining the dissociation energy of Bi 4 → 2Bi 2 and the relative energy between the T d and C 2v isomers. The single point energies, computed in the range of cc-pVnZ-PP (n = D, T, Q, 5), and with a reference geometry optimized at the CCSD(T)/cc-pVQZ-PP level of theory, were utilized to extrapolate to the complete basis set (CBS) limit. Correlation energy methods up to coupled cluster single, double, triple, and quadruple excitations (CCSDTQ) were utilized. Extrapolation to the CBS limit for the restricted Hartree-Fock energies and correlation energies were accomplished by three-point 49 and two-point 50 extrapolation schemes, respectively, with the following functions.
The determination of the correlation energy due to full triples and perturbative quadruples excitation for basis set sizes up to the extrapolated CBS limit was accomplished by assuming additivity from the perturbative triples to the full triples and full triples correlation energy to the perturbative quadruples energy computed with the cc-pVTZ-PP basis set and utilizing this same energy difference for all larger basis set sizes.
The CCSD(T) and CCSDT(Q) computations were both run in Cfour 36 utilizing the ECC 36 and NCC [51][52][53] modules, respectively.
Diagonal Born-Oppenheimer corrections 54,55 (DBOCs) and relativistic corrections were carried out at the HF/cc-pVTZ-PP and CCSD(T)/cc-pVTZ-PP levels of theory, respectively. The scalar relativistic corrections were computed through second order perturbation theory with mass-velocity and Darwin terms (MVD2). 56,57 Bi Bi Bi Bi r Bi-Bi
Bi-Bi = 2.965 Å FIG. 1: The equilibrium geometry of the T d isomer of Bi 4 predicted at the CCSD(T)/cc-pVQZ-PP level of theory.
A. Geometries
The T d is the more thoroughly studied of the two isomers of Bi 4 presented in this paper. The T d isomer has a valence electron configuration of 1a
. Numerous studies 13,14,[18][19][20]58 report an equilibrium geometry for this isomer at low levels of theory, and these structures will be compared to the geometries predicted in this research. When examining the methods used among previous studies it may be seen that the geometries computed in this research were obtained at a more rigorous level than any of its predecessors. Many of the prior studies utilize plane-wave basis sets and DFT to recover the correlation energy and these methods are less reliable than the methods used in this work. Balasubramanian and Zhang 13 used multireference methods (CASSCF and MRCISD) to analyze excited states of the T d isomer, and so the T 1 and D 1 diagnostics for this system were computed in the present research, and found to be 0.016 and 0.030, respectively. The value for the D 1 diagnostic is in an intermediary range 59 indicating that there may be some multireference character to the T d isomer. Optimizing the geometry and computing the vibrational frequencies with a high level of dynamic correllation, however, should recover multireference character. This effect is supported by the focal point energies discussed in the Energetics section of this paper, as there is fairly good convergence of the electron correlation methods. A serendipitous aspect of optimizing the T d isomer is that due to the high amount of symmetry present, the only coordinate to be optimized is the Bi-Bi bond distance. This Bi-Bi bond distance is visualized in Figure 1 as r Bi-Bi . The (most reliable) optimized r Bi-Bi bond length was computed to be 2.965 Å in this study. When this bond length is compared to the bond lengths computed in prior studies 13,14,[18][19][20] it may be seen that the r Bi-Bi predicted in this research is the smallest. Gao et al. 20 compute a value of 3.02 Å for a difference from this work of 0.05 Å. It should be noted that Gao et al. utilized the PW91 functional with plane wave basis sets, so it is among the lower levels of theory discussed here. Zhang and Balasubramanian 13 computed the r Bi-Bi value with the highest deviation from our study, that deviation being 0.15 Å. They utilized MRCISD with a triple zeta equivalent basis set including d polarization functions and a relativistic ECP that spanned the same electrons as the ECP used in the present research. Bond lengths can affect many properties of a molecule including the frequencies reported in this work and this effect will be discussed later.
When discussing the C 2v isomer, two terms will be used to describe pairs of symmetric atoms. These terms are the 'spine' and 'wing' atoms which will hence-forth refer to atoms 1 and 2, and atoms 3 and 4, respectively, of Figure 2.
While performing optimizations on the C 2v isomer it was necessary to impose a specific electronic state to keep the geometry from relaxing to the T d isomer, even though the C 2v is indeed a local minimum. It appears that the T d and C 2v isomers lie on two different potential energy surfaces and while it was not rigorously studied in this work, the interaction between these surfaces (especially in geometries intermediate to the T d and C 2v isomers) would be a fascinating avenue for future study. The valence electronic configuration of the
The C 2v isomer has been studied far less in previous research than the T d isomer. A straightforward explanation for this discrepancy is that the T d isomer is the lowest energy isomer, and when Bi 4 is generated the T d isomer is expected to be the most abundant of the two. Additionally, the C 2v isomer is on a separate electronic PES from the T d isomer. The C 2v structure was found to exhibit strong multireference character with a T 1 diagnostic of 0.027 and a D 1 diagnostic of 0.087. In light of these obstacles, only one other study 19 has reported a geometry for the neutral C 2v isomer.
As with the T d isomer the bond lengths reported by this work for the C 2v constiutional isomer are shorter than those of the previous work. 19 Akola et al. (henceforth Akola) utilized the PBE functional with a plane wave basis. Akola report an r 1 of 3.19 Å and an r 2 of 2.98 Å, both deviating from this research by 0.11 Å. The geometric data show very little difference in the angles between Akola and this work, aside from the torsional angle of the Akola study being slightly closer to planarity, implying that the geometries are approximately similar, the geometry in this work being more compact.
It should be noted that the geometry of D 2d and D 4h structures were optimized but harmonic vibrational frequency computations found an imaginary vibrational mode for each of these structures, and they are not reported as valid local minima. Both isomers preferred to form a bond between two central bismuth atoms and descend along their PES to the C 2v isomer.
The geometry and vibrational frequency of Bi 2 were computed in this research to affirm the validity of the later energetics. The bond length was computed to be 2.627 Å and the harmonic vibrational frequency was computed to be 197.6 cm -1 utilizing CCSD(T)/cc-pVQZ-PP with the same large core ECP used in the T d and C 2v Bi 4 optimizations. These predictions are in reasonable agreement with those computed by Peterson and Yousaf 60 utilizing CCSD(T)/CBS (45) with the cc-pwCVXZ-PP and a small core psuedopotential. Peterson and Yousaf predict the bond length to be 2.618 Å and the harmonic vibrational frequency to be 197.9 cm -1 . Gerber and Broida report the experimental bond length to be 2.6597 Å61 and Effantin et al. 62 deduce the Bi 2 harmonic vibrational frequency to be 173 cm -1 from experiment. 62 A bond length of 2.637 Å, in closer agreement with experiment, was obtained with the higher level CCSDT(Q)/cc-pVQZ-PP method in this research. While many studies have reported the geometry of the T d isomer, only three others 14,19,58 have reported its vibrational frequencies. This should come as no surprise because Bi 4 is a system with massive nuclei and so its vibrational potential energy wells will be shallow, making the computation of vibrational frequencies challenging. VPT2 corrections were made to the computed harmonic frequencies in this study in order to obtain fundamental frequencies for the T d structure. The previous three studies report only harmonic frequencies.
The harmonic vibrational frequencies reported by the other two studies are generally lower than those reported in this work. The three vibrational modes are computed to be 176, 101, and 131 cm -1 in this work (listed in Table I) as compared to 165, 94, and 124 cm -1 (Jia et al.) 14 and 163, 93, and (120,122) cm -1 (Akola). 19 Akola reported a splitting of 2 cm -1 in ω 2 , which is likely a slight geometric distortion due to numerical imprecision. Liang et al. report only one IR active vibrational frequency at 115 cm -1 which likely corresponds to ω 3 . 58 This trend of lower harmonic frequencies in prior works could be due to the aforementioned disparity in geometries, as the earlier theoretical geometries have longer bond distances.
The only IR-active vibrational mode for T d Bi 4 is ν 3 . This vibrational feature could be used to confirm the presence of the T d structure of the Bi 4 molecule. According to two mass spectrometry studies, 5,7 utilizing heat and laser vaporization respectively, Bi 2 and Bi 4 are produced in greatest abundance with very small amounts of Bi 3 and Bi 5 . Many studies report vibrational frequencies for Bi 3 . 19,[22][23][24]58,63 However, only one fundamental, computed by Choi et al., 63 is close enough to interfere with the detection of the T d vibrational band (ν 3 ). This vibrational frequency has a value of approximately 124 cm -1 . According to Choi et al. this fundamental is not IR active (infrared intensity of 0). Additionally, Choi et al. report very low infrared intensities for the rest of the Bi 3 frequencies, the largest intensity being 0.01 km mol -1 for modes with frequencies of about 161 and 163 cm -1 . In conjuction with the low abundance of Bi 3 and Bi 5 , it is unlikely that any other vibrational frequencies will overlap with ν 3 for T d Bi 4 and its detection will confirm the presence of the T d structure.
Anharmonic corrections are computed for the T d isomer and found to very slightly lower each harmonic frequency. The corrections display a linear increase in absolute value as the harmonic frequency energies get larger. The harmonic frequencies computed in this work appear to be a good approximation for the frequencies of the T d isomer as the largest anharmonic correction is -1 cm -1 . TABLE II: Harmonic and anharmonic vibrational frequencies for the C 2v Bi 4 (in cm -1 ) structure at the CCSD(T)/cc-pVQZ-PP level of theory. Harmonic IR intensities are reported in parentheses (in km mol -1 ).
The C 2v isomer vibrational frequencies have been reported in only one other study. 19 Akola reported vibrational frequencies of 142, 76, 58, 89, 124, and 104 cm -1 for modes one through four, while this work predicts them to be 161, 93, 54, 111, 154, and 117 cm -1 (listed in Table II). There is a broad range of deviation between this study and Akola, with most of Akola's frequencies being smaller aside from ω 3 . This ω 3 vibrational mode corresponds to a 'flapping' motion of the C 2v structure. This disparity in the trend could stem from the difference in torsional angle between the two geometries, as the Akola geometry is slightly closer to planarity than the geometry of this work. Another possible contributing factor could be that the geometry of this work is wider set, as the 'wing' atoms are predicted to be 4.5 Å apart with a torsion angle (τ 4213 ) of 57.8°and Akola predict this distance to be 4.6 Å with a τ 4213 of 57.7°. Our structure is virutally identical to Akola.
Similar to the T d isomer, the anharmonic corrections for the C 2v isomer are small and negative. The largest anharmonic correction is -1 cm -1 .
In addition to geometries and harmonic frequencies, the energy difference between the two isomers and the dissociation energy of T d Bi 4 into two Bi 2 molecules have been computed. Both of these values were obtained by employing the focal point analysis method discussed in the Theoretical Methods section.
From Table III it may be seen that after extrapolating out to the CCSDT(Q)/CBS level of theory the energy difference between the T d and C 2v isomers is +43.59 kcal mol -1 . This energy can then be augmented by a DBOC of -0.08 kcal mol -1 and an anharmonic ZPVE correction of -0.12 kcal mol -1 (the relativistic correction is already accounted for in the ECP) yielding a final energy difference of +43.39 kcal mol -1 . The absolute DBOCs for the T d and C 2v structures are 0.38 kcal mol -1 and 0.30 kcal mol -1 , respectively. Stability analysis was performed on the T d structure and a nearby 1 A 1 state (contained within a 3 T 1 state) was identified that may be contributing to a Jahn-Teller distortion, leading to a large DBOC for both structures. This large energy difference suggests that the T d isomer would exist in vast excess compared to the C 2v isomer.
The dissociation energy of Bi 4 into two Bi 2 molecules is shown in Table IV. This dissociation energy is extrapolated out to the same level of theory as the conformational energy change and the result is an energy difference of +79.24 kcal mol -1 . A large DBOC of -0.36 kcal mol -1 (the relativistic correction again being negligible) can be applied alongside an anharmonic ZPVE correction of -0.53 kcal mol -1 to yield a corrected dissociation energy of 78.35 kcal mol -1 .
The vertical electronic excitation energies for both isomers were computed and found to be 170 kcal mol -1 for the electronic transition from the T d electron configuration to the C 2v electron configuration at the T d equilibrium geometry computed in the present research. An energy difference of 13 kcal mol -1 was found for the electronic transition from the C 2v electron configuration to the T d electron configuration at the C 2v equilibrium geometry. Thus there must be a crossing of these potential energy surfaces somewhere between these minima. While this phenomenon is not explored further in the present research, a future study could be done that examines the transition between the two isomers using multireference methods to compute both diabatic and adibatic potential energy surfaces between the isomers.
Natural Bond Orbital (NBO) analyses were performed on both the T d and C 2v isomers. The results obtained for the T d isomer showed that the six bonds were almost doubly occupied (1.98) with four lone pairs being essentially doubly occupied (2.00). This bonding pattern matches exactly the conventional view of a group 5 T d tetramer, such as P 4 .
The C 2v structure, however, yields a fascinating result. The two bismuth atoms at the tips of the 'wings' (atoms 3 and 4 in Fig. 2 of the C 2v structure) have a "non-Lewis" bond connecting them with an occupancy of 0.54. We attribute this "non-Lewis" 32 interaction to long-bonding, a phenomenon recently explored by Landis and Weinhold. 64 This long-bonding orbital has primarily atomic p-orbital character. Landis and Weinhold stress the necessity of lowering the energy of the long bond through an intermediary electron density donating source. In the C 2v isomer of Bi 4 , a Lewis bond between the 'spine' atoms, atoms 1 and 2 in Fig. 2, appears to be the primary donator of electron density. This deduction was made due to the aforementioned Lewis bond being geometrically in the correct position to donate electron density to the 'long-bond' and the occupancy of the orbital denoted in Figure 3(a) being 1.60. These 'wing' atoms are 4.488 Å apart, approximately 1.5 Å longer than the T d bond distance of 2.965 Å. The NBO described by the in-phase porbitals on atoms 3 and 4 has an occupancy of 1.80, and it is reasonable to say that this NBO also donates a small amount of electron density to the long bond. This is due to its geometric and energetic proximity to the long-bond, as well as the other higher energy NBOs not having enough collective occupancy to fully explain the in-phase NBO's occupancy deficiency.
This C 2v structure 'long-bond' appears to be the energy lowering factor that distinguishes the C 2v struc- ture as a local minimum where neither the D 4h nor D 2d structures are local minima. The extreme distance between the non-bonding bismuth atoms and the lack of an electron density donator for the D 4h and D 2d structures do not allow them to take advantage of this energy lowering effect. While no computations were performed on the Bi 4 anion in this work, it may be speculated that the C 2v isomer which is reported to be the global minimum of the Bi 4 anion by Gausa et al., 16 is so because more electron density may then be donated into the long-bond, thus lowering the energy of this isomer. This does not explain, however, why the D 2d isomer was reported as a global maximum for this anion in the study by Jia et al., 14 and more work should be done to assert the validity of the C 2v structure 'long-bond'. The T d isomer is likely not the global minimun for the Bi 4 anion because its extra electron density will be donated into the antibonding orbital that pushes the d isomer into the C 2v isomer by splitting a Bi-Bi bond.
Two constitutional isomers of Bi 4 of T d and C 2v symmetry were explored in this study. These two structures were found to be local minima while the D 4h and D 2d isomers were determined to be first-order saddle points. The optimized geometry of the T d isomer can be fully described by the Bi-Bi interatomic radius, which is computed to be 2.965 Å and this leads to a more compact structure than reported in prior studies. The harmonic vibrational frequencies of the T d isomer are computed to be, in general, larger than those of prior studies. The present research is the first to predict fundamental vibrational frequencies of Bi 4 . The optimized geometry of the C 2v isomer computed in this research was found to be more compact than the geometry from a prior study, 19 while the two theoretical geometries have nearly identical conformational structures. The C 2v isomer is found to have higher harmonic vibrational frequencies in all modes but one, the ω 3 mode. The T d and C 2v isomer energy difference was computed as +43 kcal mol -1 , with the T d isomer being lower in energy. The dissociation of Bi 4 (T d ) into two Bi 2 molecules was computed to have an energy of +79 kcal mol -1 . Through NBO analysis, the C 2v isomer of Bi 4 was found to exhibit 'long-bonding' between the furthest apart 'wing' atoms. This 'long-bonding' is facilitated by the σ-bonding orbital between the 'spine' atoms of the C 2v isomer.
Finally, we point to the obvious way to spectroscopically observe the bismuth tetramer. The t 2 symmetry ν 3 fundamental predicted at 130 cm -1 has a moderate IR intensity, namely 0.20 km mol -1 . This Bi 4 fundamental is well separated from the previously observed Bi 2 (173 cm -1 ) and Bi 3 (161, 163 cm -1 ) vibrational frequencies.
We acknowledge Dr. Kevin Moore and Professor Kirk Peterson for their helpful comments and insight. A special thanks goes to Professor Frank Weinhold for his insight on NBO and long-bonding. This research was supported by the U.S. National Science Foundation under grant number CHE-1661604.