Vibronic coupling in benzene cation and anion: Vibronic coupling and frontier electron density in Jahn-Teller molecules

Vibronic coupling constants of Jahn-Teller molecules, benzene radical cation and anion, are computed as matrix elements of the electronic part of the vibronic coupling operator using the electronic wave functions calculated by generalized restricted Hartree-Fock and state-averaged complete active space self-consistent-ﬁeld methods. The calculated vibronic coupling constants for benzene cation agree well with the experimental and theoretical values. Vibronic coupling density analysis, which illustrates the local properties of the coupling, is performed in order to explain the order of magnitude of the coupling constant from view of the electronic and vibrational structures. This analysis reveals that the couplings of the e 2 g (cid:1) 2 (cid:2) and e 2 g (cid:1) 3 (cid:2) modes in which the large displacements locate on C–C bonds are strong in the cation. On the other hand, they are greatly weakened in the anion because of the decrease of electron density in the region of the C–C bonds, which originates from the antibonding nature of the singly occupied molecular orbital of the anion. However, the difference of the electronic structure has a little inﬂuence on the vibronic coupling of the e 2 g (cid:1) 4 (cid:2) mode. These results indicate that the vibronic coupling depends not only on the direction of the nuclear


I. INTRODUCTION
Vibronic interaction or electron-vibration ͑phonon͒ interaction is one of the most investigated problems in molecular physics [1][2][3][4] since it plays an important role not only in solid state physics but also in chemical reaction theory, for instance, Jahn-Teller ͑JT͒ effect, 5 superconductivity, electron transfer reaction, and so on.Therefore, vibronic coupling constant ͑VCC͒ which rules the magnitude of the interaction has been calculated by many researchers.
For benzene cation, Johnson 6 has obtained VCC by the quadratic fit to the adiabatic potential energy curve.Döscher et al. 7 have calculated the derivative of the adiabatic potential energy at the Jahn-Teller origin R 0 ͑see Fig. 1͒.Applegate and Miller 8 calculated from the distortion vector.In the degenerate electronic state, symmetry breaking 9 in the electronic wave function, which makes the coupling matrix symmetry broken, is inevitable unless the state-averaged calculation is employed.In our previous work, 10 this has been resolved by calculating the vibronic coupling integrals as matrix elements of the electronic part of vibronic coupling operator over wave functions obtained by state-averaged calculation.Furthermore, we have proposed a concept of vi-bronic coupling density which illustrates the local property of vibronic coupling.An analysis with the vibronic coupling density was found to be useful in order to elucidate the magnitude of the interaction.
In the present work, we focus on the vibronic coupling of Jahn-Teller molecules, C 6 H 6 + ͑Bz + ͒ and C 6 H 6 − ͑Bz − ͒, in their electronic ground states, 2 E 1g and 2 E 2u , respectively.͑see Figs. 2 and 3͒.][24][25][26] Though Bz − has also been an example of Jahn-Teller effect in experiment 15,[27][28][29] and theory, [30][31][32] vibronic structure has not been observed.Since the vibrational structures of Bz + and Bz − could be almost the same, comparing their VCCs we can discuss the effect of their electronic structures on the coupling.
Since the electronic ground states of Bz + and Bz − with D 6h symmetry are degenerate, the electronic motion couples to the vibrational degrees of freedom due to the JT effect.Such JT-active mode is deduced from the symmetric products of electronic states, Among 3N −6=30 molecular vibrations, ⌫ vib = 2a 1g a 2g 2b 2g e 1g 4e 2g a 2u 2b 1u 2b 2u 3e 1u 2e 2u , ͑2͒ the vibrational modes which couple to the electronic states are the two a 1g and four e 2g modes.VCC of single Slater determinant is equal to the sum of the integrals over molecular orbital which is called orbital vibronic coupling integral ͑OVCI͒.The concept of OVCI was first proposed by Bersuker more than two decades ago, [33][34][35] and has been widely used in solid state and molecular physics.Interesting researches using the OVC have been piled up in many field, for instance, charge transfer reaction, 36 mixed-valence systems, 37 superconductivity, 38 Peierls distortion, 39 and so on. 4In the case of the coupling between the degenerate electronic state and degenerate modes, VCC is equal to the OVCI over frontier molecular orbital. 10Thus, VCCs of Bz + and Bz − are equal to OVCIs over e 1g ͓highest occupied molecular orbital ͑HOMO͔͒ and e 2u ͓lowest unoccupied molecular orbital ͑LUMO͔͒, respectively, as discussed in Sec.III.Hence, VCCs of Bz + and Bz − can directly be affected by the orbital patterns of e 1g and e 2u .Vibronic coupling density analysis will be suitable for such an analysis because it has the advantage of illustrating the local properties of the vibronic coupling from the electronic structure.
This paper is organized as follows: in Sec.II, we will describe the vibronic Hamiltonian, Sec.III is devoted to the computational details, and Sec.IV deals with the result of calculation using generalized restricted Hartree-Fock ͑GRHF͒ and state-averaged complete active space selfconsistent-field ͑CASSCF͒ wave functions.We will also discuss the relation between the magnitude of vibronic coupling and the orbital patterns of the frontier orbitals with the vibronic coupling density analysis in this section.Finally, we conclude this work in Sec.V.

II. VIBRONIC HAMILTONIAN
The model Hamiltonian employed here is 10 where r is a set of electronic coordinates, Q i a normal coordinate of the ith mode, H e an electronic Hamiltonian, U a sum of an electron-electron, electron-nuclear, nuclearnuclear potential operator, and i a frequency of the ith mode.Only the linear vibronic coupling is considered.The operator involved in the third term, is called the electronic operator of the ith mode, and Q i V i describes the vibronic coupling.We treat the model Hamil-tonian within the space spanned by the electronic states ͉E͘ and ͉E⑀͘,

͑6͒
The coupling matrix can be reduced using Wigner-Eckart theorem as where the integrals ͗E͉V ͑j͒ ͉E͘ and so on are called VCIs.VCC is defined as which is the quantity we will calculate in this article.On the other hand, the interaction matrix between vibrational a 1g mode k is written as is a VCC for a totally symmetric a 1g mode.The VCC can be obtained from a calculation of the VCI.Finally, we introduce here some dimensionless quantities.For a vibrational mode i, the normal coordinate Q i is measured by ͱ ប / i , where q i is a dimensionless normal coordinate.Dimensionless coupling constant D i can be defined as

III. METHOD OF CALCULATION
As a reference structure R 0 of the JT system, we take the structure of neutral benzene, which is called parent system throughout this article.Since the parent system does not give rise to any JT distortion, its geometry has D 6h symmetry.
The electronic operator V i is a sum of one-electron operators v i ͑a͒ and the derivative of the nuclear-nuclear repulsion potential U nn with respect to Q i , where indices ␣ and a denote nucleus and electron and Z ␣ charge of the nucleus ␣.It should be noted that ‫ץ‬U nn / ‫ץ‬Q i is zero except for the a 1g modes.
For the HF methods, the wave functions for E states are written as a single Slater determinant, for Bz + , and for Bz − , and ␣, ␤ are spin functions.
The vibronic coupling matrix for the HF wave function is given by

͑18͒
For the e 2g modes, the VCI over Slater determinants ͉͘ and ͉⑀͘ can be decomposed into OVCIs for the mode over molecular orbitals, where m runs over the occupied molecular orbitals with E 1g , E 2g , E 1u , and E 2u symmetries, and ͑m͒ and ⑀͑m͒ denote the degenerate pair of molecular orbitals.n ͑m͒ and n ͑m͒ are occupation numbers of ͑m͒ and ⑀͑m͒ .Note that the orbitals with the E 1g , E 2g , E 1u , and E 2u symmetries can couple to the e 2g modes since Furthermore, it should be noted here that Therefore, the vibronic coupling matrix in the HF methods is equal to the OVCI matrix, for Bz + , and for Bz − .As for the a 1g modes, there is no such cancellation, where m and mЈ run over the occupied spin orbitals of ͉͘ and ͉⑀͘, respectively.
Here, we define the scaled dimensionless vibronic coupling constant DЈ, where is a scaling parameter which satisfies the relation, where ⌬E is the energy difference between the minimum structure of the radical optimized within D 2h symmetry and the structure of the conical intersection optimized within D 6h symmetry.Therefore, D j Ј can be calculated from V j , j , and

⌬E.
In order to obtain the optimized geometry R 0 and vibrational structure of the parent system, RHF method is employed for the parent system, neutral benzene.At the geometry R 0 , we employed state-averaged CASSCF method using GAUSSIAN 03 ͑Ref.40͒ and GRHF method using CADPAC ͑Ref.41͒ to determine the wave functions of Bz + and Bz − .All calculations were performed using the 6-31G͑d , p͒ basis set.The VCI was evaluated using these wave functions.
In the vibrational analysis, positive directions of the normal coordinates are defined in Figs. 4 and 5.All quantities are given in a.u.except for bond lengths and wave numbers throughout this article.

A. Geometrical and vibrational structures
The optimized symmetries of the parent system C 6 H 6 and that of the conical intersection for the radical ions are D 6h .In Table I, the bond lengths of the neutral benzene and the radical ions with the D 6h symmetry are tabulated.It is found that the optimized geometries of the energy minimum for benzene and the conical intersection for the radical ions are almost the same.The bond lengths differ within 0.02 Å for C-C and 0.01 Å for the C-H bonds.Thus we take the geometry of neutral benzene as that of the JT crossing point of radical ions in order to save computational resources.
Vibrational frequencies are summarized in Table II.Vibrational structures are also the same.Therefore, we took the geometrical and vibrational structures of neutral benzene as those of the JT crossing structure of the radical ions throughout this study.

B. Vibronic coupling constant
The calculated VCCs of Bz + and Bz − using RHF, GRHF, and CASSCF methods are tabulated in Tables III and IV, respectively.We confirmed that these methods yielded appropriate wave functions with the correct symmetry. 10Note that the RHF wave function is not variationally optimized for the radicals since the RHF calculation was performed for the parent system, neutral benzene.For the a 1g modes, the RHF wave function yields quite large value, comparing with the variationally optimized wave functions.This is because all the occupied orbitals contribute to the VCCs of the a 1g modes, and the errors are accumulated, while only frontier e 1g or e 2u orbitals contribute to those of e 2g modes.

Benzene cation
In all the calculations of V for Bz + shown in Table III, it is found that the order of the coupling is e 2g ͑3͒ Ͼ e 2g ͑2͒ Ͼ e 2g ͑1͒ Ͼ e 2g ͑4͒, and this agrees with the experimental result.In other words, as long as we are interested in a qualitative aspect of the e 2g modes, we can employ the RHF wave function for the VCC calculation.
The calculated D and DЈ are tabulated in Tables V and VI, respectively.Table VI demonstrates good agreement of the present results with the calculation and experiment except for the fact that the order of e 2g ͑2͒ and e 2g ͑3͒ is interchanged.

Benzene anion
In all the calculations for Bz − shown in Table IV, V of the e 2g ͑2͒ mode has the smallest value among the four e 2g modes.The calculated VCCs of Bz − become smaller, compared with those of Bz + , and the order of VCCs of e 2g ͑1͒, e 2g ͑3͒, and e 2g ͑4͒ depends on the method of calculation.
The calculated D and DЈ are tabulated in Tables VII and VIII, respectively.e 2g ͑1͒ mode has the largest value and e 2g ͑4͒ has the smallest.This result is the same as Bz + .Our result agrees with all previous works in the point that e 2g ͑4͒ has the smallest values and the result of Lipari et al. 16 in the point that the e 2g ͑1͒ has larger coupling than e 2g ͑2͒.

Vibronic coupling density
Vibronic coupling density j ͑r͒ ͑Ref.10͒ for the e 2g ͑j͒ mode is defined by where ͑r͒ is the frontier electron density of the molecular orbital , and v ͑j͒ ͑r͒ the one-electron operator defined in Eq. ͑13͒, = e 1g for Bz + , and = e 2u for Bz − .Using j , the VCC is written as

͑29͒
The vibronic coupling density enables us to analyze the calculated VCC in terms of the electronic structure ͑r͒ and the derivative of the potential with respect to the normal coordinate v ͑j͒ ͑r͒.In Fig. 4, the JT-active e 2g modes are shown.It should be noted that the larger components of the vibrational modes lie on the carbon atoms for e 2g ͑1͒ and e 2g ͑3͒, while the hydrogen atoms in e 2g ͑2͒ or e 2g ͑4͒ have large contribution.

Benzene cation
Figure 6͑a͒ shows the contour plot of the electron density ͑r͒ of e 1g orbital calculated from the GRHF wave function.
In Fig. 7͑a͒, the derivative of the potential with respect to the e 2g ͑1͒ mode is shown.It is found that the large values are distributed around C1 and C4 atoms ͑see Fig. 2͒, because the large displacement lies on these atoms.However, there is little electron density around these atoms.Consequently, the vibronic coupling density of this mode 1 is not so large, as shown in Fig. 7͑b͒.Figure 8͑a͒ shows the derivative of the potential with respect to the e 2g ͑2͒ mode.The derivative of the potential has a large value in the region surrounded by H2-C2-C3-H3 and H6-C6-C5-H5.This is because the large displacement of this mode lies on the hydrogen atoms.Since the electron density is localized on the C2-C3 and C5-C6 bonds, the vibronic coupling density for this mode 2 gives rather large near the C2-C3 and C5-C6 bonds, as shown in Fig. 8͑b͒.
The derivative of the potential with respect to the e 2g ͑3͒ mode is shown in Fig. 9͑a͒.This mode yields the largest VCC among the four.The derivative plot shows a distribution on the C2-C3 and C5-C6 bonds.Figure 10͑a͒ shows the derivative of the potential with respect to the e 2g ͑4͒ mode.The derivative of the potential has a large value on the area near H1 and H4 atoms.This is because the large displacement of this mode lies on these atoms.The vibronic coupling density of e 2g ͑4͒ mode ͓Fig.10͑b͔͒ illustrates low density near H1 and H4 atoms.The distribution of the derivative potential does not coincide with the electron density distribution .This is the reason why the VCC of this mode is the smallest.
It is interesting to note that it is not necessary for the VCC of the mode which has large component on the carbons to be large.In other words, the direction of the nuclear displacement plays an essential role.

Benzene anion
Figure 6͑b͒ shows the contour plot of the electron density ͑r͒ of the e 2u orbital calculated from the GRHF wave function.
Figures 8͑c͒ and 9͑c͒ show the vibronic coupling densities of e 2g ͑2͒ and e 2g ͑3͒ modes, respectively.Large reduction of the electron density in the region of C2-C3 and C5-C6 bonds leads to the reduction of the vibronic coupling density in this region.Thus, VCCs of Bz − of these modes are decreased, compared with those of Bz + .
Figures 7͑c͒ and 10͑c͒ show the vibronic coupling densities of e 2g ͑1͒ and e 2g ͑4͒ modes, respectively.Since these modes have large potential derivative near the midpoint of C-H bonds, large reduction of the electron density on C-C bonds has little effect on the vibronic coupling density.Thus, VCCs of Bz − of these modes are as large as those of Bz + .
It is interesting to note that the interactions of e 2g ͑2͒ and e 2g ͑3͒ are larger in Bz + ͑ e 1g ͒ than Bz − ͑ e 2u ͒, but that of e 2g ͑4͒ are larger in Bz + .This indicates that the magnitude of the vibronic coupling depends not only on the direction of the nuclear displacement but also on the phase pattern of the orbital.

V. CONCLUSION
We calculate vibronic coupling constants of benzene radical cation and anion.The couplings are computed as matrix elements of the electronic part of the vibronic coupling operator, employing GRHF and state-averaged CASSCF electronic wave functions.Our result on benzene radical cation agrees well with the experimental and theoretical values.Analysis by the vibronic coupling density reveals that the e 2g ͑2͒ and e 2g ͑3͒ vibrational modes including C-C stretching strongly couple with electron in the cation, while weakly couples with electron in the anion because of the reduction of electron density in the C-C bonds.In contrast, the e 2g ͑4͒ mode including C-H stretching moderately couples with electron in both cation and anion since the electron in the midpoint of C-C bonds has little coupling.This indicates that the magnitude of the vibronic coupling depends not only on the direction of the nuclear displacement but also on the phase pattern of the orbital.In other words, the difference on VCC between Bz + and Bz − originates from the frontier electron density.This suggests that we can control a VCC by modifying its electronic structure.Studies along this line are in progress.The results will be published elsewhere.

FIG. 1 .
FIG. 1. Cross section of the Jahn-Teller potential.Jahn-Teller crossing R 0 is the nuclear configuration of the molecule without Jahn-Teller distortion, and energy minimum R min is the molecular structure with the lowest energy.The D 6h symmetry at R 0 is lowered into the D 2h at R min because of the Jahn-Teller effect.For the structure R 0 Ј which is obtained after the conical inter- section search, the energy difference E e ͑R 0 Ј͒ − E e ͑R min ͒ = ប • D 2 / 2 is called Jahn-Teller stabilization energy ⌬E, where D is dimensionless vibronic coupling constant.

FIG. 7 .
FIG. 7. ͑a͒ Contour map of the one-electron vibronic coupling operator v ͑r͒ with respect to the e 2g ͑1͒ mode on the plane z = 1.0 a.u.͑b͒ Contour map of the vibronic coupling density e 2g ͑1͒ ͑r͒ for Bz + on the plane z = 1.0 a.u.͑c͒ Contour map of the vibronic coupling density e 2g ͑1͒ ͑r͒ for Bz − on the plane z = 1.0 a.u.

FIG. 9 .
FIG. 9. ͑a͒ Contour map of the one-electron vibronic coupling operator v ͑r͒ with respect to the e 2g ͑3͒ mode on the plane z = 1.0 a.u.͑b͒ Contour map of the vibronic coupling density e 2g ͑3͒ ͑r͒ for Bz + on the plane z = 1.0 a.u.͑c͒ Contour map of the vibronic coupling density e 2g ͑3͒ ͑r͒ for Bz − on the plane z = 1.0 a.u.
aThe energy is estimated using Koopmans' theorem.The energy of neutral benzene is −230.7139a.u. and the orbital energy of LUMO is 0.1491 a.u.

TABLE V .
Unscaled dimensionless vibronic coupling constants D of Bz + .The vibrational vectors and frequencies employed in these calculations were obtained with RHF/6−31G͑d , p͒ for neutral benzene.Negative signs are neglected.
a VCC is calculated from derivative of one-electron orbital energy at R 0 .b VCC is calculated from derivative of the adiabatic potential energy at R 0 .c VCC is calculated from distortion vector.d VCC is calculated from quadratic fit to the adiabatic potential.

TABLE VII .
Unscaled dimensionless vibronic coupling constants D of Bz − .The vibrational vectors and frequencies employed in these calculations were obtained with RHF/6−31G͑d , p͒ for neutral benzene.Negative signs are neglected.
a VCC is calculated from the derivative of one-electron orbital energy at R 0 .FIG. 6. Contour plot of the electron density ͑r͒ for ͑a͒ e 1g ͑HOMO͒ and ͑b͒ e 2u ͑LUMO͒ calculated by the GRHF/ 6-31G͑d , p͒ method on the plane z = 1.0.