Vibronic coupling in C$_{60}^-$ anion revisited: Precise derivations from photoelectron spectra and DFT calculations

TL;DR

This study revises the vibronic coupling constants of C$_{60}^-$ using high-resolution photoelectron spectra and DFT calculations, finding weaker couplings than previously reported and resolving earlier discrepancies between theory and experiment.

Contribution

The paper provides a more accurate derivation of vibronic coupling constants for C$_{60}^-$, reducing previous overestimations and aligning theoretical calculations with experimental data.

Findings

Weaker Jahn-Teller mode couplings than earlier reports.

30% reduction in total stabilization energy.

DFT calculations agree with revised experimental constants.

Abstract

The vibronic coupling constants of C$_{60}^-$ are derived from the photoelectron spectrum measured by Wang {\it et al}. [X. B. Wang, H. K. Woo, and L. S. Wang, J. Chem. Phys., {\bf 123}, 051106 (2005).] at low temperature with high-resolutions. We find that the couplings of the Jahn-Teller modes of C$_{60}^-$ are weaker than the couplings reported by Gunnarsson {\it et al}. [O. Gunnarsson, H. Handschuh, P. S. Bechthold, B. Kessler, G. Gantef{\"{o}}r, and W. Eberhardt, Phys. Rev. Lett., {\bf 74}, 1875 (1995).]. The total stabilization energy after $h_g$ and $a_g$ modes is reduced with respect to the previous derivation of Gunnarsson {\it et al}. by 30 \%. The computed vibronic coupling constants using DFT with B3LYP functional agree well with the new experimental constants, so the discrepancy between theory and experiment persistent in the previous studies is basically solved.

Tables

Table 1. Table 1: Absolute values of dimensionless vibronic coupling constants obtained in the present work.

	Frequency	PES					B3LYP333 The percentage 20 %, 25 %, and 30 % indicate fractions of the Hartree–Fock exact exchange taken in the exchange-correlation functional.
	(cm-1)	Wang111 (1),(2),and (3) are derived from the PES of Wang et al. (Ref. Wang et al., 2005)			Gunnarsson222 (4),(5) are derived from the PES of Gunnarsson et al. (Ref. Gunnarsson et al., 1995)		6-311G(d)			6-311+G(d)	cc-pVTZ
							20 %	25 %	30 %	20 %	20 %	25 %
		(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
$a_g\left(1\right)$	496	0.505	0.505	0.500	0.141	0.505	0.287	0.272	0.269	0.346	0.289	0.286
$a_g\left(2\right)$	1470	0.100	0.200	0.300	0.424	0.200	0.415	0.445	0.460	0.455	0.430	0.450
$h_g\left(1\right)$	273	0.500	0.500	0.490	0.740	0.820	0.436	0.437	0.444	0.426	0.442	0.452
$h_g\left(2\right)$	437	0.525	0.520	0.515	0.860	0.690	0.498	0.504	0.508	0.479	0.498	0.494
$h_g\left(3\right)$	710	0.465	0.460	0.455	0.390	0.350	0.418	0.464	0.476	0.412	0.403	0.414
$h_g\left(4\right)$	774	0.310	0.310	0.300	0.490	0.490	0.259	0.241	0.243	0.252	0.273	0.283
$h_g\left(5\right)$	1099	0.285	0.280	0.280	0.320	0.300	0.211	0.233	0.241	0.211	0.212	0.217
$h_g\left(6\right)$	1250	0.220	0.230	0.235	0.190	0.160	0.126	0.169	0.178	0.126	0.125	0.124
$h_g\left(7\right)$	1428	0.490	0.470	0.435	0.320	0.430	0.398	0.414	0.433	0.392	0.398	0.415
$h_g\left(8\right)$	1575	0.295	0.285	0.260	0.350	0.410	0.338	0.335	0.345	0.330	0.333	0.343

Table 2. Table 2: Stabilization energies (meV) obtained in the present work.

	Frequency	PES					B3LYP333 The percentage 20 %, 25 %, and 30 % indicate fractions of the Hartree–Fock exact exchange taken in the exchange-correlation functional.
	(cm-1)	Wang111 (1),(2),and (3) are derived from the PES of Wang et al. (Ref. Wang et al., 2005)			Gunnarsson222 (4),(5) are derived from the PES of Gunnarsson et al. (Ref. Gunnarsson et al., 1995)		6-311G(d)			6-311+G(d)	cc-pVTZ
							20 %	25 %	30 %	20 %	20 %	25 %
		(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
$a_g\left(1\right)$	496	7.8	7.8	7.7	0.6	7.8	2.5	2.3	2.2	3.7	2.6	2.5
$a_g\left(2\right)$	1470	0.9	3.6	8.2	16.4	3.6	15.7	18.0	19.3	18.9	16.9	18.5
$h_g\left(1\right)$	273	4.2	4.2	4.1	9.3	11.4	3.2	3.2	3.3	3.1	3.3	3.5
$h_g\left(2\right)$	437	7.5	7.3	7.2	20.0	12.9	6.7	6.9	7.0	6.2	6.7	6.6
$h_g\left(3\right)$	710	9.5	9.3	9.1	6.7	5.4	7.7	9.5	10.0	7.5	7.2	7.6
$h_g\left(4\right)$	774	4.6	4.6	4.3	11.5	11.5	3.2	2.8	2.8	3.1	3.6	3.8
$h_g\left(5\right)$	1099	5.5	5.3	5.3	7.0	6.1	3.0	3.7	4.0	3.0	3.0	3.2
$h_g\left(6\right)$	1250	3.9	4.1	4.3	2.8	2.0	1.2	2.2	2.5	1.2	1.2	1.2
$h_g\left(7\right)$	1428	21.3	19.6	16.8	9.1	16.4	14.0	15.2	16.6	13.6	14.0	15.2
$h_g\left(8\right)$	1575	8.5	7.9	6.6	12.0	16.4	11.2	11.0	11.6	10.7	10.9	11.5
$E_s$		8.7	11.4	15.9	17.0	11.4	18.2	20.3	21.5	22.6	19.5	21.0
$E_{\mathrm{JT}}$		65.0	62.3	57.7	78.4	82.1	50.2	54.4	57.8	48.4	49.9	52.6
$E_s+E_{\mathrm{JT}}$		73.7	73.7	73.6	95.4	93.5	68.4	74.7	79.3	71.0	69.4	73.6

Table 3. Table 3: The lowest vibronic levels (cm -1 ) and the statistical weights p k ​ J subscript 𝑝 𝑘 𝐽 p_{kJ} (%) at 70 K and 90 K. J 𝐽 J is the magnitude of the vibronic angular momentum. To calculate the vibronic levels, the set of VCCs (1) in Table 1 is used.

Level	$J$	Energy	Weight
			70 K	90 K
1	1	-962.85	97.75	92.48
2	3	-713.97	1.37	4.04
3	2	-683.40	0.52	1.77
4	1	-672.75	0.25	0.90
Sum			99.89	99.19

Table 4. Table 4: Residuals of the experimental and simulated spectra. The calculation of the residual is performed for all sets of VCCs in Table 1 at 70 K and 90 K.

	Set (1)		Set (2)		Set (3)
	70 K	90 K	70 K	90 K	70 K	90 K
$R\times {10}^{-4}$	8.35	8.07	8.41	8.16	8.76	8.41

Table 5. Table 5: Comparison between the obtained stabilization energies of C − 60 superscript subscript absent 60 {}_{60}^{-} with previous results (meV).

	Freq.	PES		LDA		GGA	MNDO	B3LYP
	(cm-1)	(3)111Results given in Tables 2.	Gun Gunnarsson et al. (1995)	Man Manini et al. (2001)	Bre Breda et al. (1998)	FreFrederiksen et al. (2008)	Var Varma et al. (1991) 222The stabilization energies of the $a_g$ modes are not reported.	Sai Saito (2002) 222The stabilization energies of the $a_g$ modes are not reported.	Laf Laflamme Janssen et al. (2010)	(10)111Results given in Tables 2.
$a_g\left(1\right)$	496	7.7	0.6	0.2	1.5	1.5	-	-	1.8	2.6
$a_g\left(2\right)$	1470	8.2	16.4	2.7	9.0	11.0	-	-	16.4	16.9
$h_g\left(1\right)$	273	4.1	11.4	2.7	6.0	2.8	1.8	3.6	3.5	3.3
$h_g\left(2\right)$	437	7.2	24.0	6.3	15.6	7.0	0.6	6.6	6.5	6.7
$h_g\left(3\right)$	710	9.1	7.8	5.5	6.6	6.1	0.6	6.6	7.1	7.2
$h_g\left(4\right)$	774	4.3	10.8	2.4	3.0	2.4	0.0	3.0	3.1	3.6
$h_g\left(5\right)$	1099	5.3	7.2	2.6	3.6	2.6	3.6	3.0	3.0	3.0
$h_g\left(6\right)$	1250	4.3	3.0	1.5	1.8	1.9	0.0	1.8	1.3	1.2
$h_g\left(7\right)$	1428	16.8	10.2	9.0	9.6	9.0	20.4	13.2	13.8	14.0
$h_g\left(8\right)$	1575	6.6	13.8	8.2	4.8	8.8	6.6	10.2	10.6	10.9
$E_s$		15.9	17.0	2.9	10.5	12.5	-	-	18.2	19.5
$E_{\mathrm{JT}}$		57.7	88.2	38.2	51.0	40.6	33.6	48.0	48.9	49.9
$E_s+E_{\mathrm{JT}}$		73.6	105.2	41.1	61.5	53.1	-	-	67.1	69.4

Table 6. Table 6: The computed vibronic levels (cm -1 ) and statistical weights p k ​ J subscript 𝑝 𝑘 𝐽 p_{kJ} (%) using the sets of VCCs (4), (5) which are derived from the experimental spectrum of Gunnarsson et al . Gunnarsson  et al. ( 1995 ) The statistical weights are calculated at 200K.

Level	Set (4)			Set (5)
	$J$	Energy	Weight	$J$	Energy	Weight
			200 K			200 K
1	1	-1067.6	36.56	1	-1134.4	37.00
2	3	-845.2	17.22	3	-914.7	17.77
3	2	-785.9	8.03	2	-850.2	7.98
4	1	-771.3	4.34	1	-829.9	4.14
5	3	-697.5	5.95	3	-742.5	5.15
6	2	-616.0	2.37	2	-687.8	2.48
7	5	-613.1	5.10	5	-685.5	5.37
8	1	-590.0	1.18	1	-667.0	1.28
9	3	-580.3	2.56	3	-650.1	2.65
10	2	-535.1	1.32	2	-595.0	1.27
11	4	-531.4	2.06	4	-594.0	2.02
12	1	-525.4	0.74	1	-590.6	0.74
13	3	-485.2	1.29	3	-543.5	1.23
14	5	-464.8	1.75	5	-515.2	1.58
15	4	-445.7	1.11	4	-495.6	1.00
16	2	-404.6	0.52	2	-456.7	0.47
17	7	-371.4	1.22	7	-446.6	1.31
18	4	-362.7	0.61	3	-436.7	0.57
19	3	-362.1	0.53	4	-432.7	0.63
20	5	-339.6	0.71	5	-414.8	0.77
21	0	-330.5	0.06	0	-400.2	0.06
22	2	-317.4	0.28	4	-388.3	0.46
23	4	-316.8	0.44	2	-387.9	0.29
24	1	-291.3	0.14	1	-357.6	0.14
25	6	-275.4	0.53	6	-336.5	0.52
Sum			96.62			96.88

Table 7. Table 7: The splitting between the lifted one-electron t 1 ​ u subscript 𝑡 1 𝑢 t_{1u} levels Δ ​ ϵ Δ italic-ϵ \Delta\epsilon (meV).

	6-311G(d)			6-311+G(d)	cc-pVTZ
	20 %	25 %	30 %	20 %	20 %	25 %
$\Delta ϵ$	691	874	1059	686	693	876

Table 8. Table 8: Stabilization energies (meV) derived from the experimental C-C bond lengths of C 60 and C − 60 superscript subscript absent 60 {}_{60}^{-} .

	Frequency	Stabilization energies of $a_g$ modes
	(cm-1)	(meV)
		(3)	(10)	Narymbetov111The structure of C60 is taken from neutron diffraction at 5K (Ref. David et al., 1991).	Fujiwara111The structure of C60 is taken from neutron diffraction at 5K (Ref. David et al., 1991).	Fujiwara222The structure of C60 is taken from X-ray diffraction at 110K (Ref. Bürgi et al., 1992).
				7K	25K	90K
$a_g\left(1\right)$	496	7.7	2.6	8.0	0.1	58.8
$a_g\left(2\right)$	1470	8.2	16.9	0.0	7.5	7.0

Equations

$$\left[T_{1u}^2\right]=a_g\oplus h_g.$$

$$H$$

$$\hat{I}$$

$${\hat{C}}_{-1}$$

$${\hat{C}}_1$$

$$\left\{|m_{\mathrm{el}}⟩|\cdots {\mathbf{n}}_{\mu }\cdots ⟩|v_1{v}_2{⟩}_{a_g}\right\}.$$

$$|{\Psi }_{v_1{v}_{2}kJM}⟩$$

$$g_{\Gamma (\mu )}$$

$$S_{v^{\prime }v}(g)$$

$$\left\{|m_{\mathrm{el}}⟩|\cdots {\mathbf{n}}_{\mu }\cdots ⟩;\sum _{\mu =1}^8\sum _{m=-2}^2{n}_{\mu m}\le N\right\}.$$

$$E_{s,i}$$

$$E_{\mathrm{JT},\mu }$$

$$E_s$$

$$E_{\mathrm{JT}}$$

$$I(\Omega )$$

$$p_{v_i}$$

$$p_{kJ}$$

$$F(\Omega )={\int }_{-\infty }^{\infty }I({\Omega }^{\prime })\exp \left[\frac{-(\Omega -{\Omega }^{\prime }{)}^2}{2{\sigma }^2}\right]d{\Omega }^{\prime }.$$

$$R$$

$$V_{a_g(i)}$$

$$=$$

$$V_{h_g(\mu )}$$

$$=$$

$$V_{a_g(i)}=-\sum _{A=1}^{60}\left({\mathbf{R}}_A-{\mathbf{R}}_{0,A}\right)\cdot \frac{{\mathbf{u}}_A^{a_g(i)}}{\sqrt{M}}.$$

$$E(\mathbf{R})$$

$$V_{{\Gamma }^{\prime }{\gamma }^{\prime }}$$

Full text

I INTRODUCTION

Much attention has been paid to the Jahn–Teller effect of fullerene (C60) in various electronic states not only because the Jahn–Teller effect is an interesting problem in molecular physicsChancey and O’Brien (1997) but also because it is expected to play an important role in the mechanism of the superconductivity in alkali-doped fullerides. Gunnarsson (2004) Thus, the strength of the electron-vibration coupling (vibronic coupling) of C60 which characterizes the Jahn–Teller effect has been one of the important topics. The vibronic coupling constants (VCCs) have been estimated experimentally

Gunnarsson et al. (1995); Winter and Kuzmany (1996); Hands et al. (2008)

and theoreticallyVarma et al. (1991); Schluter et al. (1992); Faulhaber et al. (1993); Antropov et al. (1993); Breda et al. (1998); Manini et al. (2001); Saito (2002); Frederiksen et al. (2008); Laflamme Janssen et al. (2010).

In the experimental studies of vibronic coupling in fullerene, a landmark is the photoelectron spectroscopy (PES) of C${}_{60}^{-}$ in gas phase by Gunnarsson et al.Gunnarsson et al. (1995) As C${}_{60}^{-}$ is one of the most studied systems, in addition to this experimental work, computational works have been performed by many authors. However, discrepancy between the coupling constants of the experimental and theoretical works have been reported. Gunnarsson et al. (1995); Manini et al. (2001); Saito (2002) The theoretical stabilization energies as estimated by density functional theory (DFT) calculation were always obtained smaller than that derived from the experiment of Gunnarsson et al. Besides uncertainties intrinsic to the DFT method, whose predictions depend on the used exchange-correlation functional, one should note that the derivation of vibronic coupling constants in Ref. Gunnarsson et al., 1995 is not perfect either. First, the thermal excitations were not included in the simulation, although the vibrational temperature of C${}_{60}^{-}$ was estimated about 200 K in the experiment.

Second, not all vibronic coupling constants have been estimated from the spectrum because of the low resolution.

For the same reason, the computed VCCs of totally symmetric modes were used to simulate the PES.

Recently, Wang et al. remeasured photoelectron spectra of C${}_{60}^{-}$. Wang et al. (2005) In their experiment, the vibrational temperature of C${}_{60}^{-}$ is between 70 K and 90 K and the resolution is about 16 meV, i.e., much smaller than the resolution of 40 meV in the experiment of Gunnarsson et al. Accordingly, the spectrum of Wang et al. is narrower and has more structures, therefore, it is expected to yield more reliable coupling constants.

In this work, we simulate the photoelectron spectra of Wang et al. Wang et al. (2005) and Gunnarsson et al.Gunnarsson et al. (1995) and give new derivations of the VCCs of C${}_{60}^{-}$. We also compute the VCCs of C${}_{60}^{-}$ using the DFT method and compare them with the experimental values.

II THE SOLUTION OF THE JAHN–TELLER PROBLEM OF C${}_{60}^{-}$

The equilibrium geometry of neutral fullerene is taken as the reference nuclear configuration. At this reference structure, the ground electronic state of C${}_{60}^{-}$ is $T_{1u}$. According to the selection rule, the $T_{1u}$ electronic state couples with two $a_{g}$ and eight $h_{g}$ vibrational modes:

[TABLE]

We consider the linear $T_{1u}\otimes(2a_{g}\oplus 8h_{g})$ Jahn–Teller Hamiltonian. The Hamiltonian is written as follows:

[TABLE]

where $Q_{\Gamma(\mu)m}$ is the mass-weighted normal coordinate of $m$ element of the $\Gamma(\mu)$ mode ($\Gamma=a_{g},h_{g}$), $P_{\Gamma(\mu)m}$ is the conjugate momentum of the normal coordinate $Q_{\Gamma(\mu)m}$, $\omega_{\Gamma(\mu)}$ is the frequency of the $\Gamma(\mu)$ mode, $V_{\Gamma(\mu)}$ is the VCC of the $\Gamma(\mu)$ mode, and $\hat{I}$ and $\hat{C}_{-m}$ are the $3\times 3$ unit matrix and a matrix whose elements are Clebsch–Gordan coefficients, respectively. The normal modes and frequencies of C60 are used for C${}_{60}^{-}$ , so the higher vibronic terms which mix the normal modes of fullerene are neglected. As a $T_{1u}$ electronic basis set $\{|m_{\rm el}\rangle;m_{\rm el}=-1,0,1\}$ and normal coordinates of the $h_{g}$ modes $\{Q_{h_{g}(\mu)m};m=-2,-1,0,1,2\}$, we use complex basis which transform as spherical harmonics $\{Y_{1m_{\rm el}};m_{\rm el}=-1,0,1\}$ and $\{Y_{2m};m=-2,-1,0,1,2\}$, respectively, under the rotations. O’Brien (1971); Auerbach et al. (1994) Then $\hat{I}$ and $\hat{C}_{-m}$ are written as Edmonds (1974)

[TABLE]

This type of the Jahn–Teller problem was investigated by O’Brien O’Brien (1969) and the vibronic coupling constants defined by her are often used. Thus we introduce the coefficient $\sqrt{5/2}$ in front of the vibronic term to make $V_{h_{g}}$ the same as O’Brien’s coupling constants.

Since the linear $T_{1u}\otimes(2a_{g}\oplus 8h_{g})$ Jahn–Teller Hamiltonian (2) commutes with squared vibronic angular momentum $\mathbf{J}$ and the $z$ component of $\mathbf{J}$,Bersuker and Polinger (1989) the eigenstate of the Hamiltonian (2) is the simultaneous eigenstate of the vibronic angular momentum $J$, the $z$ component of the vibronic angular momentum $M$. Here, the vibronic angular momentum $\mathbf{J}$ is the sum of the vibrational angular momentum $\mathbf{L}$ and the “energy spin” $\mathbf{S}$ describing the threefold orbital degeneracy ($S=1$).Bersuker and Polinger (1989) In the case of linear vibronic coupling, the eigenstate of $H$ is the product of the $T_{1u}\otimes(8h_{g})$ Jahn–Teller part and the $a_{g}$ vibrational part. As a vibronic basis a set of the products of electronic states and vibrational states of the $a_{g}$ and $h_{g}$ modes is used:

[TABLE]

Here, $\mathbf{n}_{\mu}$ means a set of vibrational quantum numbers of the $h_{g}(\mu)$ mode, $\mathbf{n}_{\mu}=\left\{n_{\mu m}\right\}$, $v_{1},v_{2}$ are vibrational quantum numbers of the $a_{g}(1)$ mode and the $a_{g}(2)$ mode respectively. Then the eigenstate $|\Psi_{v_{1}v_{2}kJM}\rangle$

of the Hamiltonian (2) which belongs to the eigenvalue

$\sum_{i=1}^{2}\left[\hslash\omega_{a_{g}(i)}v_{i}-V_{a_{g}(i)}^{2}/(2\omega_{a_{g}(i)}^{2})\right]+E_{kJ}$

is represented as a linear combination of the vibronic basis with constants $C^{\rm JT}_{m_{\rm el},\mathbf{n}_{1}\cdots\mathbf{n}_{8};kJM}$.

[TABLE]

where $E_{kJ}$ is an eigenvalue of the $T_{1u}\otimes(8h_{g})$ Jahn–Teller Hamiltonian, $k$ distinguishes energy levels with the same $J$ and $M$, the dimensionless VCC of the $\Gamma(\mu)$ mode $g_{\Gamma(\mu)}$ ($\Gamma=a_{g},h_{g}$) is defined as

[TABLE]

the Franck–Condon factor of the $a_{g}$ mode $S_{v^{\prime}v}(g)$ is written as

[TABLE]

$l_{\rm min}=0$ for $v\leq v^{\prime}$, and $l_{\rm min}=v-v^{\prime}$ for $v>v^{\prime}$. The origin of the energy is the lowest energy of C${}_{60}^{-}$ without vibronic couplings.

To obtain the vibronic states, we diagonalize the linear $T_{1u}\otimes 8h_{g}$ Jahn–Teller Hamiltonian numerically using Lanczos method. We use a truncated vibronic basis set,

[TABLE]

Here, $N$ is the maximum number of the vibrational excitations in the vibronic basis set (8). We treat the vibronic states which $J$s are from [math] to $7$. Frequencies $\omega_{a_{g}(i)},\omega_{h_{g}(\mu)}$ are taken from the experimental frequencies of Raman scattering in solid state C60.Bethune et al. (1991)

Lastly, we introduce stabilization energies which we use to show our results. The stabilization energy of each mode is defined as

[TABLE]

and the total stabilization energies of the $a_{g}$ modes and $h_{g}$ modes are

[TABLE]

They represent the depth of the potential energy surface from the energy of undistorted fullerene monoanion.

III SIMULATION OF THE PHOTOELECTRON SPECTRUM

The photoelectron spectrum is simulated within the sudden approximation. Hedin and Lundqvist (1969) We assume that each C${}_{60}^{-}$ is in a thermal equilibrium state, hence we use a Boltzmann’s distribution to calculate the statistical weight. With these assumptions, the intensity of the transition which appears at the binding energy $\hslash\Omega$ is written as follows:

[TABLE]

where, $p_{v_{i}}$ and $p_{kJ}$ are the statistical weights of the $a_{g(i)}$ mode and the Jahn–Teller part, respectively,

[TABLE]

$Z_{a_{g}}$ and $Z_{\rm JT}$ are corresponding statistical sums, and $E_{0}$ is the gap between the ground electronic energies of C60 and C${}_{60}^{-}$. The envelope function is represented by using the Gaussian function with the standard deviation $\sigma$:

[TABLE]

For a decent simulation of experimental PES one should include in Eq. (16), in principle, also the contributions from the rotational spectrum of C${}_{60}^{-}$. However, due to a large momentum of inertia of fullerene and restrictive selection rules for the transitions between different rotational levels Child and Longuet-Higgins (1961) our estimations gave an expected enlargement of the transition band of only several wave numbers. This is negligible compared the full width at half maximum (FWHM) given by the envelope function (16).

To evaluate the agreement between the simulated spectrum and the experimental spectrum, we calculate the residual of theoretical spectrum $F_{\rm calc}(\Omega)$ and the experimental spectrum $F_{\rm exp}(\Omega)$. The residual $R$ is defined by the equation:

[TABLE]

Here, $f$ is the parameter to vary the height, $\Omega_{\rm shift}$ is the parameter to shift the experimental spectrum, $\Omega_{i}$ is a sampling point. The minimum and maximum of $\Omega_{i}$ is $\Omega_{\rm min}$ and $\Omega_{\rm max}$, and the gap between adjacent sampling points $\Delta\Omega$ is constant. Then $M$ is represented as $M=(\Omega_{\rm max}-\Omega_{\rm min})/\Delta\Omega$ and $\Omega_{i}=\Omega_{\rm min}+i\Delta\Omega$. In the calculation of the residual $R$, $\Omega_{\rm min}$, $\Omega_{\rm max}$, and $\Delta\Omega$ are $-$200 cm*-1*, 1600 cm*-1*, and 0.5 cm*-1* respectively. We avoid the truncation of the zero phonon line of Gunnarsson et al.Gunnarsson et al. (1995) VCCs are varied in order to make $R$ as small as possible within the accuracy of the experiment. The accuracy is determined from the range of the vibrational temperature of C${}_{60}^{-}$ in the experiment of Wang et al. Wang et al. (2005) In their experiment, the vibrational temperature is between 70 K and 90 K. Although the shapes of the simulated spectra at 70 K and 90 K are different from each other, we cannot distinguish them from the experiment of Wang et al. In terms of the residual $R$, the difference between $R$ at 70 K and $R$ at 90 K is practically indistinguishable.

IV DFT CALCULATION OF VIBRONIC COUPLING CONSTANTS

The linear vibronic coupling constant of the $a_{g}(i)$ mode is a diagonal matrix element of the first derivative of the electronic Hamiltonian with respect to the normal coordinate at the reference geometry. Bersuker and Polinger (1989)

[TABLE]

where $\psi$ is the ground electronic state. By applying the Hellmann–Feynman theorem Feynman (1939) to Eq. (18) and then transforming it into the formula with the vibrational vector, we obtain

[TABLE]

Here, $A$ indicates a carbon atom in C60, $\mathbf{R}_{A}$ is the Cartesian coordinate of $A$, $\mathbf{R}$ is the set of all $\mathbf{R}_{A}$, $H_{\rm el}(\mathbf{R})$ is the electronic Hamiltonian at the structure $\mathbf{R}$, $\mathbf{R}_{0}$ is the reference nuclear configuration, $E(\mathbf{R})$ is the ground electronic energy $\langle\psi|H_{\rm el}(\mathbf{R})|\psi\rangle$, $M$ is the mass of carbon atom, $\mathbf{u}_{A}^{a_{g}(i)}$ is the vibrational vector of the $a_{g}(i)$ mode. Similarly, absolute value of the coupling constant of the $h_{g}(\mu)$ mode is written as

[TABLE]

The equilibrium geometry $\mathbf{R}_{0}$, the vibrational vectors $\mathbf{u}^{a_{g}(i)}$, $\mathbf{u}^{h_{g}(\mu)m}$, and the gradient of the electronic energy $\left(\partial E(\mathbf{R})/\partial\mathbf{R}_{A}\right)_{\mathbf{R}_{0}}$, entering the Eqs. (20) and (22), are obtained from ab initio calculations.

Note that the vibronic coupling constants (19), (21) are not equal to the gradients of the frontier levels (see Appendix).

We compute the VCCs of C${}_{60}^{-}$ using the DFT method. As exchange-correlation functional, the hybrid functional of Becke Becke (1993) (B3LYP) is used. To find VCCs which are close to the experimental results the fraction of the Hartree–Fock exchange energy are varied from the original fraction 20% to 30% by 5%. We use the triple zeta basis sets, 6-311G(d), 6-311+G(d), and cc-pVTZ.

The structure optimization and the calculation of the vibrational modes are performed for the neutral fullerene. The electronic wavefunction of C${}_{60}^{-}$ are obtained from the variational calculation of an unrestricted Slater determinant. As far as the method based on the single determinant is used, the spatial symmetry of the wavefunction is broken and the degeneracy of the singly occupied degenerate level is lifted. Sato et al. (2006, 2009) However, in the case of the cyclopentadienyl radical, it was demonstrated that the splitting of the total electronic energies estimated by the unrestricted B3LYP method is only 0.4 meV. Sato et al. (2006, 2009) It is expected that the splitting of the $T_{1u}$ ground electronic energies of C${}_{60}^{-}$ is tiny and the symmetry of the electronic state may not be broken significantly. Thus we treat the wavefunction as a $T_{1u}$ wavefunction. We calculate the energy gradient $\left.\partial E(\mathbf{R})/\partial\mathbf{R}_{A}\right|_{\mathbf{R}\rightarrow\mathbf{R}_{0}}$ with the coupled perturbed Kohn–Sham method. In the calculation of the dimensionless VCCs (6) and the stabilization energies (10), we use the experimental frequencies.Bethune et al. (1991) To compute electronic structures we use the Gaussian 03 program.Frisch et al. (2004)

V Derivation of the vibronic coupling constants of $a_{g}$ modes from

the structures of C60 and C${}_{60}^{-}$

We also derive the stabilization energies of the $a_{g}$ modes from the experimental bond lengths of C60 and C${}_{60}^{-}$. The structures of C60 and C${}_{60}^{-}$ with $I_{h}$ symmetry are determined by the C-C bond lengths of the edges between two hexagons (6:6) and a hexagon and pentagon (6:5). We use average 6:6 and 6:5 C-C bond lengths of TDAE-C60 for C${}_{60}^{-}$ and fullerite for C60. The data of TDAE-C60 are obtained from the results of X-ray diffraction at 7K by Narymbetov et al. Narymbetov et al. (1999) and at 25K and 90K by Fujiwara et al.Fujiwara et al. (2005). The average bond lengths of fullerite are taken from the results of neutron diffraction at 5K by David et al.David et al. (1991) and X-ray diffraction at 110K by Bürgi et al.Bürgi et al. (1992) To remove the thermal expansion of the C-C bond lengths, we use sets of bond lengths of C60 and C${}_{60}^{-}$ which are measured at close temperature. That is, the bond lengths of C60 measured at 5K is used with the bond lengths of C${}_{60}^{-}$ measured at 7K and 25K, and the bond lengths of C60 at 110K is used with the bond lengths of C${}_{60}^{-}$ at 90K. The vibronic coupling constants of the $a_{g}$ modes $V_{a_{g}(i)}(i=1,2)$ are obtained from the equation

[TABLE]

To perform the calculation, we use the vibrational vector defined in the calculations with the B3LYP method and the cc-pVTZ basis.

VI RESULTS AND DISCUSSIONS

VI.1 Simulation of the PES of Wang et al.

We simulate the photoelectron spectrum measured by Wang et al. Wang et al. (2005) at 70 K. The basis set in Eq. (LABEL:Eq:vibrostate) includes up to 6 vibrational excitations ($N=6$). The experimental and simulated spectra are shown in Fig. 1. The transition between the ground states of C${}_{60}^{-}$ and C60 (the 0–0 line) 0- is chosen as the origin of these spectra. From the spectrum of Wang et al., we obtain several sets of VCCs listed as (1), (2), and (3) in Tables 1 and 2.

We extracted $\sigma=80$ cm*-1* by fitting the FWHM of the 0–0 line (188 cm*-1*). The increase or decrease of the $\sigma$ makes the agreement between the simulated and experimental spectra worse.

To assess the thermal population of the excited vibronic states, we calculate statistical weights of the excited Jahn–Teller levels $p_{kJ}$ at 70 K and 90 K. The vibronic levels are obtained using the set of VCCs (1). In the calculation of the distribution function $Z_{\rm JT}$, we include all excited vibronic levels whose weights are larger than $10^{-7}$. The computed weights are shown in Table 3. Although these statistical weights are computed using the set (1), rest of the sets of VCCs (2), (3) give similar results. The statistical weights of the ground vibronic level at 70 K and 90 K are more than 90 %. This indicates that the transition from the ground vibronic level is dominant in the PES of Wang et al. We focus, therefore, on the ground vibronic level to discuss the effect of the size of the basis (8) on the calculated vibronic states. The ground vibronic level is $-962.65$ cm*-1* when we use the basis set with $N=5$. Compared with the gap between the ground and first excited vibronic levels with $J=1$, the change of the ground vibronic level due to the increase of the size of the vibronic basis set is only about 0.07 %. Therefore, we regard our basis set as large enough to simulate the spectrum of Wang et al.

The differences between several sets of VCCs are in the constants of $a_{g}(2)$, $h_{g}(7)$, and $h_{g}(8)$ modes. If we increase the dimensionless VCC (the stabilization energy) of the $a_{g}(2)$ mode from 0.1 to 0.3 (0.9 to 8.2 meV) and at the same time decrease the dimensionless VCCs of $h_{g}(7),h_{g}(8)$ modes, the shape of the PES does not vary significantly (see Fig. 2). This is due to a poor resolution of the peaks of $a_{g}(2),h_{g}(7)$, and $h_{g}(8)$ modes and essentially the same problem arose in the analysis of Gunnarsson et al.Gunnarsson et al. (1995) In the latter case, the stabilization energy of $a_{g}(2)$ is varied from 0 to 45 meV, i.e., in a range larger than ours. Owing to the narrow peaks of the spectrum of Wang et al., we can derive the VCCs with less ambiguity.

We compute the residuals (LABEL:Eq:R) of the experimental and the simulated spectra for all sets of VCCs at 70 K and 90 K. The values are shown in Table 4.

The differences of the residuals of the different sets of VCCs are almost within the ambiguity of the vibrational temperature. Therefore, we conclude that these sets of VCCs cannot be distinguished from the experiment of Wang et al.

Although we obtain several sets of VCCs, the stabilization energies $E_{s}+E_{\rm JT}$ are similar to each other (see (1), (2), and (3) in Table 2). On the other hand, present stabilization energies are smaller than the stabilization energy of Gunnarsson et al. Gunnarsson et al. (1995) by 30 % (see Table 5), i.e., the Jahn–Teller coupling is weaker than previously expected. We find that the distributions of $E_{s,i}$ and $E_{{\rm JT},\mu}$ also differ from each other (Tables 2, 5). In Ref. Gunnarsson et al., 1995, the stabilization energy of $h_{g}(2)$ was found the strongest, while our results show that the strongest is the stabilization energy of $h_{g}(7)$.

Hands et al. estimated the Jahn–Teller stabilization energy $E_{\rm JT}$ of 57.94 meV within the single-mode $T_{1u}\otimes h_{g}$ Jahn–Teller model from the visible and near-infrared spectrum. Hands et al. (2008) The present Jahn–Teller stabilization energy agrees well with their value.

VI.2 Simulation of the PES of Gunnarsson et al.

As a preliminary calculation, we compute the vibronic levels using the data from Ref. Gunnarsson et al., 1995, that is, the same VCCs and the same size of the vibronic basis ($N=5$). The statistical weight of the ground state at 200 K is obtained ca. 39 %. This result indicates that not only the ground level but also excited levels must be considered in order to simulate the spectrum of Gunnarsson et al. Gunnarsson et al. (1995)

We simulate this spectrum at 200 K with the FWHM of 283 cm*-1* ($\sigma=120$ cm*-1*). The size of the vibronic basis set is $N=7$. The experimental and the simulated spectra are shown in Fig. 3. As was mentioned also by Gunnarsson et al.,Gunnarsson et al. (1995) we obtain several sets of VCCs that give close stabilization energies. These sets of dimensionless VCCs and stabilization energies are (4), (5) in Table 1 and 2 respectively. In comparison with the original stabilization energy of Gunnarsson et al., present stabilization energies are smaller by 10 meV. However, stabilization energies are still larger than those obtained from the spectrum of Wang et al. by 20 meV. The inconsistencies of these VCCs come from the difference between the shapes of the spectra of Wang et al. and Gunnarsson et al., due to different vibrational temperature and resolution.

Simulating the spectrum of Gunnarsson et al., we encounter two problems. First, the spectrum is too broad and, second, the vibrational temperature is too high. In fact, the statistical weights of the ground vibronic level at 200 K are about 37 % in both cases (see Table 6), hence, we must consider many excited vibronic states. To represent the excited vibronic states with enough accuracy, we expect that the vibronic basis must be larger than the present one. Furthermore, as the vibrational temperature is further increased, the weight of each vibronic level and the shape of the spectrum varies easily. Although the range of the vibrational temperature is not reported, we increased the temperature by 20 K which is the uncertainty range of vibrational temperature in the case of Wang et al. Wang et al. (2005) The statistical weights of the ground vibronic level decreased from ca. 37 % to ca. 31 % with this increase of the temperature. This change of the weight affects the shape of the spectrum. Therefore it is difficult to perform an accurate simulation and to estimate VCCs from the spectrum of Gunnarsson et al.

Given better the experimental conditions of Wang et al. Wang et al. (2005) allowing for more accurate simulations, we may conclude that the VCCs extracted from these experiments should be considered more reliable than those obtained by Gunnarsson et al. Gunnarsson et al. (1995)

VI.3 DFT calculations of the vibronic coupling constants

We compute the vibronic coupling constants of C${}_{60}^{-}$ using the DFT method described in Sec. IV.

For DFT calculations with pure functionals, it is well known that in the high-symmetry geometry the occupied level belonging to one degenerate representation moves upwards in energy relative to the empty levels belonging to the same degenerate manifold. Bruyndonckx et al. (1997) However in our case the situation is opposite (the occupied level is the lowest one) because of the Hartree–Fock exchange contribution contained in the B3LYP functional. The splitting between the $t_{1u}$ Kohn–Sham levels are about 1 eV (see Table 7) and the variations of the total electronic energies for different occupation schemes of $t_{1u}$ orbitals are less than 0.2 meV. Moreover, the vibronic coupling constants do not depend on the choice of the electronic states significantly. The variation of the total stabilization energy is ca. 1 meV.

The dimensionless VCCs and the stabilization energies are shown in Table 1 and 2. Although we use several basis sets, the VCCs do not depend on the basis set significantly. On the other hand, the VCCs vary with the increase of the fraction of the Hartree–Fock exchange energy in the exchange-correlation functional. Increasing this fraction leads to larger VCCs and stabilization energies. We find that the stabilization energies of high frequency modes, $a_{g}(2),h_{g}(7)$, and $h_{g}(8)$ are the strongest. Compared with other DFT calculations, we may conclude that the stabilization energies of $a_{g}$ modes agree well with the previous calculations, while the present stabilization energy $E_{s}+E_{\rm JT}$ is larger than the previous results.

In comparison with present simulation of the experimental PES, the DFT calculations with the energy functionals including fractions of 20 % and 25 % of the Hartree–Fock exchange energy give close values. Although the stabilization energy $E_{s}+E_{\rm JT}$ obtained using the original B3LYP functional is slightly smaller than the experimental value, the result obtained with it is also close to the experimental one. The distribution of the computed stabilization energy of each $h_{g}$ mode $E_{{\rm JT},\mu}$ qualitatively agrees with the experimental results. The stabilization energy of the $h_{g}(7)$ mode is obtained smaller and that of the $h_{g}(8)$ mode is obtained larger than the experimental values. The slight difference between theoretical and experimental results should originate from still inaccurately computed vibrational vectors. Indeed, it was shown that a small mixing of the vibrational vectors in fullerene affects the values of VCCs significantly. Gunnarsson et al. (1995) Besides the present computational results, the LDA calculation by Manini et al.Manini et al. (2001) and the GGA calculation by Frederiksen et al.Frederiksen et al. (2008) give similar relative values for the coupling constants of $h_{g}$ modes as the present simulations of PES of Wang et al. Wang et al. (2005) On the other hand, these calculations give smaller absolute values of the VCCs and of the total stabilization energies than the presently obtained. Contrary to these calculations, the B3LYP calculation by Saito Saito (2002) and Laflamme Janssen Laflamme Janssen et al. (2010) give close values of VCCs to the present results. However, despite of using the same B3LYP functional, our calculations of VCCs differ from the ones in Refs. Saito, 2002 and Laflamme Janssen et al., 2010 since we used here the derivative of the total energy, which does not coincide with the derivative of the Kohn–Sham orbital energy (see the appendix).

The LDA calculation of Breda et al.Breda et al. (1998) gives the distribution of relative strengths of VCCs which is similar to the results of Gunnarsson et al. and do not agree with the values derived here from PES of Wang et al. Varma et al.Varma et al. (1991) computed VCCs using MNDO method, however, the distribution of the stabilization energies is different from the present simulations of experiment and the theoretical values obtained here.

On the contrary, the theoretical stabilization energy of the $a_{g}(1)$ mode is too small and that of the $a_{g}(2)$ mode is too large compared to $E_{s,i}$ derived from experiment (Tables 1, 2). To find the correct order of the corresponding VCCs, we derive the coupling constants of the $a_{g}$ modes from the experimental bond lengths as described in Sect. V. The obtained stabilization energies $E_{s,i}$ are shown in Table 8.

Unfortunately, as we can see, $E_{s,i}$ depend very strongly on the set of the C-C bond lengths, thus we cannot draw a conclusion about their relative strength. This comes from the fact that we use the structural data measured by different techniques (X-ray and neutron scattering) on different systems (TDAE-C60 and fullerite) for C${}_{60}^{-}$ and C60, respectively. In both cases, fullerenes should be deformed due to the environment compared with the free C60 molecule. The distortions caused by the crystal fields of C60 in TDAE-C60 and fullerite are different from each other. Furthermore, the expected changes in the bond lengths in C60 and C${}_{60}^{-}$ are within the experimental accuracy of the structural data. Thus, the origin of the discrepancy of the relative strength of VCCs of the $a_{g}$ modes between the simulations and DFT calculations remains unclear.

Note that although the values of $E_{s,i}$ are different from the experimental results, their sum, as well as, the distribution of $E_{{\rm JT},\mu}$ and $E_{s}+E_{\rm JT}$ are close to the experimental values. Therefore, we may conclude that the present theoretical method gives improved values of vibronic coupling constants.

VII CONCLUSION

In this work, we simulated the PES of Wang et al. and derived the vibronic coupling constants of C${}_{60}^{-}$. We obtain several sets of VCCs, because the frequencies of $a_{g}(2)$, $h_{g}(7)$, and $h_{g}(8)$ modes are close to each other. Considering the ambiguity of the vibrational temperature in the experiment, these sets of VCCs cannot be distinguished. Thus, to obtain more accurate coupling constants, it is desired to perform an observation of a PES of C${}_{60}^{-}$ in still better experimental conditions. Although we find several sets of VCCs from the spectrum, the stabilization energies are similar to each other. In comparison with the total stabilization energy derived by Gunnarsson et al., Gunnarsson et al. (1995) our value is smaller by 30 %. We also calculated the VCCs using the DFT method. Even though the experimental and theoretical orders of $E_{s,i}$ disagree with each other, the distribution of $E_{{\rm JT},\mu}$ and the total stabilization energy $E_{s}+E_{\rm JT}$ agrees well with the experimental values. Thus we may conclude that the problem of the discrepancy between the experimental and calculated coupling constants, persistent in the previous studies, is basically solved in the present work. As an extension of the present work we expect that the theoretical approach used here could be successfully applied for the calculation of VCC of C${}_{60}^{n-}$ anions in AnC60 fullerides as well as of their multiplet splitting parameters.

Appendix A Vibronic coupling constants and the gradient of Kohn–Sham levels

The total energy $E(\mathbf{R})$ in the DFT is written as follows:

[TABLE]

Here, $\Gamma\gamma$ is the irreducible representation of the Kohn–Sham orbital, $\mu$ is the quantum number other than $\Gamma\gamma$, $\epsilon_{\mu\Gamma\gamma}$ is the Kohn–Sham level, $\sum_{\mu\Gamma\gamma}$ is taken over occupied levels, $n(\mathbf{r};\mathbf{R})$ is the ground electronic density, $E_{xc}\left[n(\mathbf{r};\mathbf{R})\right]$ is the exchange-correlation energy functional, $V_{xc}(\mathbf{r},\mathbf{R})$ is the exchange-correlation potential, and $V_{nn}(\mathbf{R})$ is the Coulomb potential energy between nuclei. The vibronic coupling constant of $\Gamma^{\prime}\gamma^{\prime}$ mode is

[TABLE]

For the totally symmetric modes, all the derivatives in the right-hand side of Eq. (25) are not zero. For the Jahn–Teller active modes, the sum of the gradient of the completely occupied Kohn–Sham levels belonging to the same $\Gamma$ is zero due to the symmetry reasons, and the gradient of the Coulomb potential between the nuclei is zero also because of the symmetry. However, the second and third terms which include the derivative of $n(\mathbf{r};\mathbf{R})$ and $V_{xc}(\mathbf{r};\mathbf{R})$ with respect to $Q_{\Gamma^{\prime}\gamma^{\prime}}$, respectively, are nonzero because $n(\mathbf{r};\mathbf{R})$ is not a totally symmetric function. Therefore, in general, the vibronic coupling constant is not equal to the gradient of the frontier Kohn–Sham level.

Acknowledgements.

We would like to thank Dr. X. B. Wang and Prof. L. S. Wang for sending us unpublished data. N.I. would like to thank the research fund for study abroad from the Research Project of Nano Frontier, graduate school of engineering, Kyoto University. T.S. and N.I. are grateful to the Division of Quantum and Physical Chemistry at the University of Leuven for hospitality. Theoretical calculations were partly performed using Research Center for Computational Science, Okazaki, Japan. This work was supported by a Grant-in-Aid for Scientific Research, priority area ‘Molecular theory for real systems’ (20038028) from the Japan Society for the Promotion of Science (JSPS). This work was also supported in part by the Global COE Program ‘International Center for Integrated Research and Advanced Education in Materials Science’ (No. B-09) of the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, administered by the JSPS. Financial support from the JSPS–FWO (Fonds voor Wetenschappelijk Onderzoek-Vlaanderen) bilateral program is gratefully acknowledged.