First principles study of the vibronic coupling in positively charged C$_{60}^+$

TL;DR

This study uses DFT calculations to analyze vibronic coupling and Jahn-Teller effects in positively charged C60, revealing significant stabilization energies and the importance of dynamical effects in understanding its electronic structure.

Contribution

The paper provides the first detailed DFT-based analysis of vibronic coupling parameters and Jahn-Teller effects in C60+, highlighting the role of dynamical effects and comparing with C60-.

Findings

Jahn-Teller deformation with D5d symmetry identified

Stabilization energies of 110 and 129 meV calculated

Dynamical Jahn-Teller effect is crucial in C60+

Abstract

Vibronic coupling parameters for C$_{60}^{+}$ were derived via DFT calculations with hybrid B3LYP and CAM-B3LYP functional, based on which the static Jahn-Teller effect were analyzed. The global minima of adiabatic potential energy surface (APES) shows a D$_{5d}$ Jahn-Teller deformation, with stabilization energies of 110 and 129 meV (with B3LYP and CAM-B3LYP respectively), which are two times larger than that in C$_{60}^-$, suggesting the crucial role of the dynamical Jahn-Teller effect. Present results enable us to assess the actual situation of dynamical Jahn-Teller effect in C$_{60}^{+}$ and excited C$_{60}$ in combination with the established parameters for C$_{60}^{-}$.

Tables

Table 1. Table 1: The frequencies ω Γ subscript 𝜔 Γ \omega_{\Gamma} (cm -1 ), orbital vibronic coupling parameters v Γ subscript 𝑣 Γ v_{\Gamma} (10 -4 a.u.), and stabilization energies E Γ subscript 𝐸 Γ E_{\Gamma} (meV) for the a g subscript 𝑎 𝑔 a_{g} , g g subscript 𝑔 𝑔 g_{g} , and h g subscript ℎ 𝑔 h_{g} modes. The data calculated with B3LYP for LUMOs are taken from Ref. Liu  et al. ( 2018 ) . g Γ = V Γ / ℏ ​ ω Γ 3 subscript 𝑔 Γ subscript 𝑉 Γ Planck-constant-over-2-pi superscript subscript 𝜔 Γ 3 g_{\Gamma}=V_{\Gamma}/\sqrt{\hslash\omega_{\Gamma}^{3}} is dimensionless vibronic coupling parameter.

	B3LYP							CAM-B3LYP
$\mu$	${\omega }_{\Gamma }$	$v_{\Gamma }$		$g_{\Gamma }$		$E_{\Gamma }$		${\omega }_{\Gamma }$	$v_{\Gamma }$		$g_{\Gamma }$		$E_{\Gamma }$
		1	2	1	2	1	2		1	2	1	2	1	2
HOMO
$a_{g1}$	497	$-0.121$		$-0.112$		0.389		506	$-0.118$		$-0.106$		0.355
$a_{g2}$	1498	1.040		0.184		3.159		1527	1.369		0.236		5.268
$g_{g1}$	481	$-0.375$		$-0.366$		3.984		502	$-0.405$		$-0.370$		4.265
$g_{g2}$	584	$-0.101$		$-0.074$		0.196		586	$-0.207$		$-0.150$		0.815
$g_{g3}$	768	0.923		0.446		9.466		787	1.083		0.504		12.400
$g_{g4}$	1092	$-1.281$		$-0.365$		9.019		1092	$-1.583$		$-0.451$		13.781
$g_{g5}$	1335	0.541		0.114		1.076		1356	0.739		0.152		1.947
$g_{g6}$	1540	1.477		0.251		6.029		1576	1.571		0.258		6.509
$h_{g1}$	266	0.690	$-0.080$	1.635	$-0.190$	44.099	0.593	271	0.727	$-0.090$	1.674	$-0.206$	47.120	0.714
$h_{g2}$	439	$-0.508$	$-0.290$	$-0.568$	$-0.324$	8.776	2.860	450	$-0.558$	$-0.363$	$-0.601$	$-0.391$	10.063	4.269
$h_{g3}$	726	$-0.043$	0.977	$-0.023$	0.514	0.023	11.869	745	$-0.061$	1.053	$-0.031$	0.532	0.043	13.086
$h_{g4}$	786	0.923	$-0.059$	0.431	$-0.027$	9.038	0.037	802	0.930	$-0.139$	0.421	$-0.063$	8.821	0.198
$h_{g5}$	1125	$-0.101$	$-0.495$	$-0.028$	$-0.135$	0.053	1.269	1148	$-0.074$	$-0.565$	$-0.020$	$-0.150$	0.028	1.591
$h_{g6}$	1269	0.965	0.272	0.220	0.062	3.790	0.301	1300	0.986	0.025	0.216	0.005	3.768	0.002
$h_{g7}$	1443	2.860	1.487	0.537	0.279	25.745	6.960	1480	3.425	2.010	0.618	0.363	35.072	12.081
$h_{g8}$	1607	2.721	$-1.542$	0.434	$-0.246$	18.790	6.034	1663	3.217	$-1.701$	0.488	$-0.258$	24.523	6.854
LUMO
$a_{g1}$	497	$-0.264$		$-0.245$		1.849		506	$-0.253$		$-0.228$		1.629
$a_{g2}$	1498	$-2.380$		$-0.422$		16.543		1527	$-2.921$		$-0.503$		23.971
$h_{g1}$	266	0.192		0.455		3.415		271	0.209		0.481		3.884
$h_{g2}$	439	0.450		0.503		6.886		450	0.456		0.491		6.735
$h_{g3}$	726	0.754		0.396		7.069		745	0.849		0.429		8.512
$h_{g4}$	786	0.554		0.259		3.256		802	0.575		0.260		3.367
$h_{g5}$	1125	0.766		0.209		3.038		1148	0.827		0.219		3.402
$h_{g6}$	1269	0.578		0.132		1.360		1300	0.513		0.113		1.019
$h_{g7}$	1443	2.099		0.394		13.867		1480	2.553		0.461		19.492
$h_{g8}$	1607	2.043		0.326		10.592		1663	2.325		0.352		12.808

Table 2. Table 2: Jahn-Teller energies of C + 60 superscript subscript absent 60 {}_{60}^{+} (meV) for D 5 ​ d subscript 𝐷 5 𝑑 D_{5d} and D 3 ​ d subscript 𝐷 3 𝑑 D_{3d} minima of the APES, respectively.

Functional	Method	$D_{5d}$	$D_{3d}$	Ref
B3LYP	(I)	110	30	Present
CAM-B3LYP	(I)	129	39	Present
LDA	(I)	69	22	Manini et al. (2001)
B3LYP	(II)	121	-	Lykhin, Ahmadvand, and Varganov (2019)
LDA	(III)	74	27	Ramanantoanina et al. (2013)
OPBE	(III)	74	28	Ramanantoanina et al. (2013)
B3LYP	(III)	80	32	Ramanantoanina et al. (2013)
PBE	(III)	74	28	Muya et al. (2013)
LDA	(IV)	72	20	Ramanantoanina et al. (2013)
OPBE	(IV)	74	21	Ramanantoanina et al. (2013)
B3LYP	(IV)	94	25	Ramanantoanina et al. (2013)

Equations

$$[h_u\otimes h_u]=a_g\oplus g_g\oplus 2h_g.$$

$$H$$

$$H_0$$

$$H_{\text{JT}}$$

$$E_g^{\text{JT}}=-\frac{V_g^2}{2{\omega }_g^2},E_{nh}^{\text{JT}}=-\frac{V_{nh}^2}{2{\omega }_h^2},(n=1,2).$$

$$E_a$$

$$q_{\theta }=\sqrt{\frac{5}{8}}{Q}_{\theta }+\sqrt{\frac{3}{8}}{Q}_{ϵ},q_{ϵ}=\sqrt{\frac{3}{8}}{Q}_{\theta }-\sqrt{\frac{5}{8}}{Q}_{ϵ},$$

$$V_{1h}=\frac{F_{Hb}}{4\sqrt{5}},V_{2h}=\frac{F_{Ha}}{12\sqrt{5}}.$$

$$V_{h_g}^{-}$$

$$V_{\Gamma }=-v_{\Gamma },$$

$${\text{q}}_{h_g}^{D_{5d}}$$

$${\text{q}}_{g_g}^{D_{5d}}$$

$${\text{q}}_{h_g}^{D_{3d}}$$

$${\text{q}}_{g_g}^{D_{3d}}$$

$$U^{D_{5d}}$$

$$U^{D_{3d}}$$

$$E_{\text{JT}}^{D_{5d}}$$

$$E_{\text{JT}}^{D_{3d}}$$

$$\frac{V_i}{{\omega }_i^2}=c_i\sqrt{M}\Delta {\text{R}}_{\text{JT}}.{\text{e}}_i,$$

$${\hat{C}}_{ga}=\left(\begin{array}{ccccc}\frac{1}{2}\sqrt{\frac{3}{2}}& 0& 0& 0& 0\\ 0& \frac{1}{2}\sqrt{\frac{3}{2}}& 0& 0& 0\\ 0& 0& -\frac{1}{\sqrt{6}}& 0& 0\\ 0& 0& 0& -\frac{1}{\sqrt{6}}& 0\\ 0& 0& 0& 0& -\frac{1}{\sqrt{6}}\end{array}\right)$$

$${\hat{C}}_{gx}=\left(\begin{array}{ccccc}0& 0& -\frac{\varphi }{4}\sqrt{\frac{5}{2}}& 0& 0\\ 0& 0& \frac{{\varphi }^{-2}}{4}\sqrt{\frac{5}{6}}& 0& 0\\ -\frac{\varphi }{4}\sqrt{\frac{5}{2}}& \frac{{\varphi }^{-2}}{4}\sqrt{\frac{5}{6}}& 0& 0& 0\\ 0& 0& 0& 0& \frac{1}{2}\sqrt{\frac{5}{6}}\\ 0& 0& 0& \frac{1}{2}\sqrt{\frac{5}{6}}& 0\end{array}\right)$$

$${\hat{C}}_{gy}=\left(\begin{array}{ccccc}0& 0& 0& \frac{{\varphi }^{-1}}{4}\sqrt{\frac{5}{2}}& 0\\ 0& 0& 0& -\frac{{\varphi }^2}{4}\sqrt{\frac{5}{6}}& 0\\ 0& 0& 0& 0& \frac{1}{2}\sqrt{\frac{5}{6}}\\ \frac{{\varphi }^{-1}}{4}\sqrt{\frac{5}{2}}& -\frac{{\varphi }^2}{4}\sqrt{\frac{5}{6}}& 0& 0& 0\\ 0& 0& \frac{1}{2}\sqrt{\frac{5}{6}}& 0& 0\end{array}\right)$$

$${\hat{C}}_{gz}=\left(\begin{array}{ccccc}0& 0& 0& 0& \frac{1}{4}\sqrt{\frac{5}{2}}\\ 0& 0& 0& 0& \frac{5}{4\sqrt{6}}\\ 0& 0& 0& \frac{1}{2}\sqrt{\frac{5}{6}}& 0\\ 0& 0& \frac{1}{2}\sqrt{\frac{5}{6}}& 0& 0\\ \frac{1}{4}\sqrt{\frac{5}{2}}& \frac{5}{4\sqrt{6}}& 0& 0& 0\end{array}\right)$$

$${\hat{C}}_{1h\theta }=\left(\begin{array}{ccccc}\frac{\sqrt{5}}{16}& \frac{3\sqrt{3}}{16}& 0& 0& 0\\ \frac{3\sqrt{3}}{16}& -\frac{\sqrt{5}}{16}& 0& 0& 0\\ 0& 0& -\frac{{\varphi }^{-2}}{4}& 0& 0\\ 0& 0& 0& \frac{{\varphi }^2}{4}& 0\\ 0& 0& 0& 0& -\frac{\sqrt{5}}{4}\end{array}\right)$$

$${\hat{C}}_{1hϵ}=\left(\begin{array}{ccccc}\frac{3\sqrt{3}}{16}& -\frac{\sqrt{5}}{16}& 0& 0& 0\\ -\frac{\sqrt{5}}{16}& -\frac{3\sqrt{3}}{16}& 0& 0& 0\\ 0& 0& -\frac{\sqrt{3}\varphi }{4}& 0& 0\\ 0& 0& 0& \frac{\sqrt{3}{\varphi }^{-1}}{4}& 0\\ 0& 0& 0& 0& \frac{\sqrt{3}}{4}\end{array}\right)$$

$${\hat{C}}_{1h\xi }=\left(\begin{array}{ccccc}0& 0& -\frac{{\varphi }^{-2}}{4}& 0& 0\\ 0& 0& -\frac{\sqrt{3}\varphi }{4}& 0& 0\\ -\frac{{\varphi }^{-2}}{4}& -\frac{\sqrt{3}\varphi }{4}& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right)$$

$${\hat{C}}_{1h\eta }=\left(\begin{array}{ccccc}0& 0& 0& \frac{{\varphi }^2}{4}& 0\\ 0& 0& 0& \frac{\sqrt{3}{\varphi }^{-1}}{4}& 0\\ 0& 0& 0& 0& 0\\ \frac{{\varphi }^2}{4}& \frac{\sqrt{3}{\varphi }^{-1}}{4}& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right)$$

$${\hat{C}}_{1h\zeta }=\left(\begin{array}{ccccc}0& 0& 0& 0& -\frac{\sqrt{5}}{4}\\ 0& 0& 0& 0& \frac{\sqrt{3}}{4}\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ -\frac{\sqrt{5}}{4}& \frac{\sqrt{3}}{4}& 0& 0& 0\end{array}\right)$$

$${\hat{C}}_{2h\theta }=\left(\begin{array}{ccccc}\frac{9}{16}& -\frac{\sqrt{15}}{16}& 0& 0& 0\\ -\frac{\sqrt{15}}{16}& -\frac{9}{16}& 0& 0& 0\\ 0& 0& \frac{\varphi }{4}& 0& 0\\ 0& 0& 0& -\frac{{\varphi }^{-1}}{4}& 0\\ 0& 0& 0& 0& -\frac{1}{4}\end{array}\right)$$

$${\hat{C}}_{2hϵ}=\left(\begin{array}{ccccc}-\frac{\sqrt{15}}{16}& -\frac{9}{16}& 0& 0& 0\\ -\frac{9}{16}& \frac{\sqrt{15}}{16}& 0& 0& 0\\ 0& 0& -\frac{{\varphi }^{-2}}{4\sqrt{3}}& 0& 0\\ 0& 0& 0& \frac{{\varphi }^2}{4\sqrt{3}}& 0\\ 0& 0& 0& 0& -\frac{1}{4}\sqrt{\frac{5}{3}}\end{array}\right)$$

$${\hat{C}}_{2h\xi }=\left(\begin{array}{ccccc}0& 0& \frac{\varphi }{4}& 0& 0\\ 0& 0& -\frac{{\varphi }^{-2}}{4\sqrt{3}}& 0& 0\\ \frac{\varphi }{4}& -\frac{{\varphi }^{-2}}{4\sqrt{3}}& 0& 0& 0\\ 0& 0& 0& 0& \frac{1}{\sqrt{3}}\\ 0& 0& 0& \frac{1}{\sqrt{3}}& 0\end{array}\right)$$

$${\hat{C}}_{2h\eta }=\left(\begin{array}{ccccc}0& 0& 0& -\frac{{\varphi }^{-1}}{4}& 0\\ 0& 0& 0& \frac{{\varphi }^2}{4\sqrt{3}}& 0\\ 0& 0& 0& 0& \frac{1}{\sqrt{3}}\\ -\frac{{\varphi }^{-1}}{4}& \frac{{\varphi }^2}{4\sqrt{3}}& 0& 0& 0\\ 0& 0& \frac{1}{\sqrt{3}}& 0& 0\end{array}\right)$$

$${\hat{C}}_{2h\zeta }=\left(\begin{array}{ccccc}0& 0& 0& 0& -\frac{1}{4}\\ 0& 0& 0& 0& -\frac{1}{4}\sqrt{\frac{5}{3}}\\ 0& 0& 0& \frac{1}{\sqrt{3}}& 0\\ 0& 0& \frac{1}{\sqrt{3}}& 0& 0\\ -\frac{1}{4}& -\frac{1}{4}\sqrt{\frac{5}{3}}& 0& 0& 0\end{array}\right)$$

Full text

I Introduction

Recent experimental confirmation of the presence of C${}_{60}^{+}$ as interstellar materials Campbell et al. (2015, 2016) renewed the interests in C60 cations, and burst various spectroscopic and theoretical investigations on C${}_{60}^{+}$ and related systems Spieler et al. (2017); Yamada, Ross, and Ito (2017); Strelnikov et al. (2018); Kaiser et al. (2018); Lykhin, Ahmadvand, and Varganov (2019); Cordiner et al. (2019). It is known that C60 in its charged and excited states exhibits complex Jahn-Teller (JT) effect due to its high symmetry ($I_{h}$) Chancey and O’Brien (1997); Bersuker (2006). In particular, JT effect in C${}_{60}^{+}$ is one of the most involved cases because of the five-fold degenerate highest occupied molecular orbitals (HOMOs) of C60. Toward understanding JT effect in C${}_{60}^{n+}$ cations, spectroscopic Brühwiler et al. (1997); Canton et al. (2002); Kern et al. (2013) and theoretical Ceulemans and Fowler (1990); Moate et al. (1996); De Los Rios, Manini, and Tosatti (1996); Manini, Gattari, and Tosatti (2003); Lijnen and Ceulemans (2005); Hands et al. (2006a, b, 2007); Ceulemans et al. (2012) investigations have been piled up.

Realistic description of JT effect in positively charged C60 relies on the combination of an adequate model and accurate enough vibronic coupling parameters. Derivations of vibronic coupling parameters have been addressed Manini et al. (2001); Saito (2002); Ramanantoanina et al. (2013), and comprehensive sets of parameters have been estimated by density functional theory (DFT) calculations at local density approximation (LDA) level Manini et al. (2001); Ramanantoanina et al. (2013). Nevertheless, in studies about C${}_{60}^{-}$, it has been shown that LDA tends to underestimate the coupling parameters Iwahara et al. (2010), while hybrid B3LYP functional is found to give closer parameters to the experimental data. Furthermore, a good agreement between B3LYP and GW approximation calculations Faber et al. (2011) supports the accuracy of hybrid functional in studies of C60. Besides, one recent study showed that B3LYP with long-range interaction correction, CAM-B3LYP, could improve the accuracy of vibronic parameters in C${}_{60}^{-}$ with respect to experimental data, indicating CAM-B3LYP could give vibronic parameters much closer to the real situationHuang and Liu (2020). Therefore, it is desired to derive the coupling parameters at a better level than LDA for accurate description of C60 cations.

In this work, we derived orbital vibronic coupling parameters for C${}_{60}^{+}$ via DFT calculations with both B3LYP and CAM-B3LYP hybrid functionals. These obtained vibronic coupling parameters were compared with the previous data obtained by LDA calculations. Based on these parameters, the adiabatic potential energy surface (APES) was analyzed, and the symmetry of JT deformed C${}_{60}^{+}$ as well as static JT energies were established.

II Vibronic Hamiltonian

The highest occupied molecular orbitals (HOMOs) of C60 with $I_{h}$ symmetry Kroto et al. (1985) are characterized by five-fold degenerate $h_{u}$ irreducible representation. According to selection rule, these orbitals linearly couple to mass-weighted normal vibrational modes involved in the symmetric product of $h_{u}$ representation Jahn and Teller (1937):

[TABLE]

Among them, $g_{g}$ and $h_{g}$ modes are JT active, while $a_{g}$ is not because is does not change the symmetry of molecules. Thus, taking the equilibrium structure of neutral C60 as the reference, $H\otimes(a\oplus g\oplus 2h)$ JT Hamiltonian for C${}_{60}^{+}$ is expressed as Ceulemans and Fowler (1990); Chancey and O’Brien (1997); Bersuker (2006)

[TABLE]

where $\omega_{\Gamma}$ ($\Gamma=a_{g},g_{g},h_{g}$) are vibration frequencies, $q_{\Gamma}$ are mass-weighted normal coordinates Inui, Tanabe, and Onodera (1990), $V_{\Gamma}$ is vibronic coupling parameter for $\Gamma$ mode, and $\hat{C}_{\Gamma\gamma}$ ($\gamma$ = a, x, y, z, $\theta$, $\epsilon$, $\xi$, $\eta$, $\zeta$) are the Clebsch-Gordan coefficients, which are taken from Ref. Liu et al. (2018) and listed in Appendix A. A coefficient $\sqrt{5}/2$ is multiplied to vibronic couplings terms of $h_{g}$ modes so that JT energy becomes:

[TABLE]

There are two vibronic couplings to one $h_{g}$ mode because $h_{g}$ representation appears twice in selection rule, Eq. (1). The basis of vibronic Hamiltonian matrices are $H_{u}$ electronic states of C${}_{60}^{+}$ in the order of $|H_{u}\theta\rangle$, $|H_{u}\epsilon\rangle$, $|H_{u}\xi\rangle$, $|H_{u}\eta\rangle$, $|H_{u}\zeta\rangle$. For $h_{g}$ and $h_{u}$ representations, $d$ orbital type basis are used, and hence, $\theta,\epsilon,\xi,\eta,\zeta$ transform as $(2z^{2}-x^{2}-y^{2})/\sqrt{6}$, $(x^{2}-y^{2})/\sqrt{2}$, $\sqrt{2}yz$, $\sqrt{2}zx$, $\sqrt{2}xy$, respectively, under rotation. Although C60 has two $a_{g}$, six and eight sets of $g_{g}$ and $h_{g}$ modes, respectively, and indices distinguishing them are not explicitly written in Eq. (3) for simplicity. Phase factors of mass-weighted normal modes are the same as those in the Supplemental Materials of Ref. Liu et al. (2018). Since the equilibrium geometry of C60 is chosen as the reference structure of C${}_{60}^{+}$, vibronic coupling parameters of totally symmetric modes are also nonzero. JT energy by totally symmetric deformation is

[TABLE]

In many other literatures, like the work of A. CeulemansCeulemans and Fowler (1990), linear combinations of real $d$-type functions, $d_{z^{2}},d_{x^{2}-y^{2}}$, are used, which could be transformed into irreducible representation in this work by

[TABLE]

where $Q_{\theta}$ and $Q_{\epsilon}$ are normal coordinates in Ref. Ceulemans and Fowler (1990), while $q_{\theta}$ and $q_{\epsilon}$ are normal coordinates in this work, resulting in the relation between coupling parameters $V_{1h}$ and $V_{2h}$ in this work and $F_{Hb}$ and $F_{Ha}$ defined in Ref. Ceulemans and Fowler, 1990:

[TABLE]

The modification introduced here is to treat the JT Hamiltonian in a framework consistent to the standard one in C${}_{60}^{-}$ Chancey and O’Brien (1997).

III Results

III.1 Orbital vibronic coupling parameters

Vibronic coupling constants of C60 have been most intensively investigated in the case of C${}_{60}^{-}$, and coupling constants have been derived by various methods. By definition, vibronic coupling parameters for JT active $h_{g}$ modes, $V_{h_{g}}^{-}$, can be derived by Iwahara et al. (2010)

[TABLE]

where $E_{z}^{-}(\textbf{q}_{h_{g}})=\langle T_{1u}z|\hat{H}^{-}|T_{1u}z\rangle$, $|T_{1u}z\rangle$ is $t_{1u}z$ electronic state of C${}_{60}^{-}$, $\hat{H}^{-}$ is Hamiltonian for C${}_{60}^{-}$, and $\textbf{q}_{h_{g}}=\textbf{0}$ indicates equilibrium structure of C60. Because of symmetry, contributions from occupied orbitals are zero, and only partially filled $t_{1u}$ orbital level contribute to vibronic couplings. In the case of C${}_{60}^{-}$, the nature of $t_{1u}$ orbitals do not differ from that of neutral C60: although $t_{1u}$ orbitals are mixed with the other $t_{1u}$ orbitals, the mixing is very small due to high symmetry, and orbitals are separated from each other by large orbital energy gaps. Consequently, gradients of total energy can be approximated by those of orbital energy levels with respect to normal modes of neutral C60. Indeed, in the case of C${}_{60}^{-}$, these two approaches give very close results Laflamme Janssen et al. (2010); Iwahara et al. (2010); Faber et al. (2011); Liu et al. (2018).

Similar situation is expected in C${}_{60}^{+}$: the nature of the $h_{u}$ orbital does not change by adding one hole to the same molecular structure. Thus, orbital vibronic coupling parameters of neutral C60 can be used to express vibronic coupling parameters of C${}_{60}^{+}$. Since C${}_{60}^{+}$ has one hole in $h_{u}$ HOMOs, it is convenient to perform particle-hole transformation Fetter and Walecka (2003), under which, the sign of orbital vibronic coupling parameter for one electron in HOMOs of C60 should be inverted for that in the case of one hole in HOMOs. Therefore,

[TABLE]

where $v_{\Gamma}$ is one orbital vibronic coupling parameter for C60 and $V_{\Gamma}$ is the parameter for C${}_{60}^{+}$.

Orbital vibronic coupling parameters were calculated using frozen-phonon approach. Orbital energy levels of distorted C60 are fitted to the eigenvalues of JT Hamiltonian matrix (Eq. (4)). In the present case, ${h_{g}(\mu)\epsilon}$ and ${g_{g}(\mu)a}$ deformations are used because diagonalizing the model Hamiltonian is easier. For DFT calculations, a triple-zeta basis set [6-311G(d)] was employed for both B3LYP and CAM-B3LYP functionals within Gaussian Frisch et al. (2016). Some fittings are shown in Fig. 1 (see Supplemental Materials for other fittings). Black points indicate DFT levels originating from HOMOs, with gray lines for energy level calculated from model Hamiltonian. Derived orbital vibronic coupling parameters are shown in Table 1. From this table, JT stabilization energies for different JT active modes are improved about 30 % with CAM-B3LYP compared with these with B3LYP. One should note that there is almost no nonlinear splitting due to vibronic effect in HOMO levels, indicating weak quadratic or higher vibronic couplings as in the case of C${}_{60}^{-}$ Liu et al. (2018). This guarantees the validity of linear vibronic model (Eq. (3)) for the description of JT effect of C60 cations.

III.2 Static Jahn-Teller effect

Vibronic coupling lifts degeneracy with the deformation keeping one of the highest subgroup symmetries Liehr (1963), resulting in six $D_{5d}$ and ten $D_{3d}$ minima Ceulemans and Fowler (1990), as there are six $C_{5}$ and ten $C_{3}$ axes in C60. Thus, based on present basis, using symmetry adapted deformations Hands et al. (2006b), deformations for $D_{5d}$ and $D_{3d}$ minima are expressed by

[TABLE]

and

[TABLE]

respectivelyLiu et al. (2018). Substituting these symmetrized deformations, Eqs. (III.2) and (III.2), into potential terms of model Hamiltonian (kinetic energy term is ignored), and then diagonalizing model Hamiltonian, we obtain the lowest adiabatic potential energies as

[TABLE]

for $D_{5d}$ and $D_{3d}$ deformations, respectively. Furthermore, these global minima energies could be expressed in terms of stabilization energies (Eq. (5)), asCeulemans and Fowler (1990)

[TABLE]

JT stabilization energies of C${}_{60}^{+}$ are obtained using these equations and the calculated vibronic coupling parameters. When treating C${}_{60}^{+}$, we have to sum up contributions from all $g_{g}$ and $h_{g}$ modes. Manini et al. also have derived vibronic coupling parameters and JT stabilization energies in the same way with LDA. Manini et al. (2001) Besides, there is always the stabilization due to the totally symmetric modes (Eq. (6)).

JT stabilization energies of C${}_{60}^{+}$ have been calculated with different methods with various functionals. These methodologies to derive JT energies are classified into four. We denote the present method (I). In the second method, JT stabilization energy is directly obtained from the energy difference between high- and low-symmetric structures. This method was employed by Lykhin et al. Lykhin, Ahmadvand, and Varganov (2019) with B3LYP. The third and the fourth methods, (III) and (IV), are called interaction mode approach Khlopin, Polinger, and Bersuker (1978); Bersuker and Polinger (1989) and intrinsic distortion path approach Bruyndonckx et al. (1997); Zlatar, Schläpfer, and Daul (2009); Zlatar et al. (2010); Zlatar and Gruden (2019), respectively. In both methods, vibronic coupling parameters or JT energies are extracted from optimized geometry with each subgroup of $I_{h}$. The interaction mode induces deformation along JT minima from high-symmetric coordinates and is expressed by a linear combination of normal modes of high-symmetric C60. The coefficients of the linear combinations contain the information of the vibronic coupling. By expanding the difference of JT deformed and high-symmetry geometries, $\Delta\textbf{R}_{\text{JT}}$, with the eigen modes of mass-weighted normal modes $\textbf{e}_{i}$ Inui, Tanabe, and Onodera (1990) at high-symmetry structure as

[TABLE]

the vibronic coupling parameters $V_{i}$ could be obained. The frequencies $\omega_{i}$ are obtained from first principles calculations and $M$ is mass of carbon, and the coefficients $c_{i}$ depend on the structure of the JT interaction matrix. This approach was used in Refs. Ramanantoanina et al. (2013); Muya et al. (2013). On the other hand, within the intrinsic distortion path analysis, the high-symmetric structure is expressed by the linear combination of the eigen vectors from low-symmetric structure. Combining the vibronic coupling parameters derived from the deformation, $V^{\prime}_{i}$, and the frequencies at the low-symmetric structure, $\omega^{\prime}_{i}$, the JT energy could be written as $E_{\text{JT}}=\sum_{i}V_{i}^{\prime 2}/(2\omega_{i}^{\prime 2})$. The last method was applied to C${}_{60}^{+}$ in Ref. Ramanantoanina et al. (2013).

JT stabilization energies in this work, as well as those from previous studies, are shown in Table 2, from which we could see that CAM-B3LYP could enhance JT stabilization energies for $D_{5d}$ and $D_{3d}$ minima by 17% and 30% respectively compared to that with B3LYP. And JT stabilization energies obtained with both B3LYP and CAM-B3LYP are larger than those from LDA or PBE-related functionals. In particular, data calculated with LDA from Manini et al.,Manini et al. (2001) is only about 60 % of the present data. The underestimation of JT energy by LDA method is consistent with the situation in C${}_{60}^{-}$ Iwahara et al. (2010). However, since with LDA methods (II) and (III) give similar data as the one by method (I), Manini et al. (2001) the difference of these methods would not be the origin of discrepancies seen in B3LYP data. Thus, a possible reason is that the deformed geometry with B3LYP functional in Ref. Ramanantoanina et al. (2013) is the one at a local minima. For B3LYP, Ref. Lykhin, Ahmadvand, and Varganov (2019) and present work give much close results, whereas the former is larger than the latter by 11 meV. So such a difference is expected to come from that JT stabilization energies contributed by totally symmetric modes, Eq. (6), are included in Ref. Lykhin, Ahmadvand, and Varganov (2019). However, JT stabilization energy contributed from $a_{g}$ modes in this work is about 3.5 meV (Table 1), which is still smaller than 11 meV. The underestimation of vibronic coupling of totally symmetric modes may come from the lack of contributions from occupied orbitals and change of frequencies. In general, only the partially filled frontier orbitals contribute to the vibronic coupling to the JT modes due to the symmetry Sato, Tokunaga, and Tanaka (2006), whereas all occupied orbitals do contribute to vibronic coupling to the totally symmetric modes Sato, Tokunaga, and Tanaka (2008). Indeed, in the case of a planar molecule, picene, orbital vibronic coupling parameters for totally symmetric modes differ by 10-20 % from those obtained by fitting the gradients of total system energy Sato, Iwahara, and Tanaka (2012). Frequencies can be changed due to the removal of electron by 5-15 % Matsuda et al. (2018). Another possible origin of this discrepancy is nonlinear vibronic coupling, however, such coupling is much weaker than linear vibronic coupling as in the case of C60 anions Liu et al. (2018).

Besides these works, we note that Kern et al Kern et al. (2013) have optimized the structure to simulate infrared (IR) absorption spectrum with BP86 functional, but JT energy was not derived.

IV Discussion

To fully reveal the molecular nature of C${}_{60}^{+}$, non-adiabatic dynamical Jahn-Teller effect is crucial. The most straightforward way would be exact diagonalizing the molecular Hamiltonian, which fully quantizes both nuclear and electronic coordinates, nevertheless, it is not practical. To partly overcome this difficulty, combining Jahn-Teller model Hamiltonian (3) with accurate enough vibronic coupling parameters is indispensable to derive low-energy states.

Vibronic coupling parameters of C${}_{60}^{+}$ have been derived by using DFT calculations with LDA Manini et al. (2001); Ramanantoanina et al. (2013) and B3LYP Saito (2002) functionals. As discussed in Sec. II, model Hamiltonian is described by one vibronic coupling parameter for each $a_{g}$ and $g_{g}$ mode, and two parameters for each $h_{g}$ mode. Both parameters for the $h_{g}$ modes, $V_{1h}$ and $V_{2h}$, have been derived only in Ref. Manini et al. (2001) and Ref. Ramanantoanina et al. (2013), while they have not in Ref. Saito (2002). Ramanantoanina et al. Ramanantoanina et al. (2013) shows that the magnitudes of derived coupling parameters obtained from the gradient of HOMO levels Manini et al. (2001) and those from adiabatic potential energy surface agree well with each other, which has been also seen in C${}_{60}^{-}$ Iwahara et al. (2010); Laflamme Janssen et al. (2010); Liu et al. (2018). Thus, the present orbital vibronic coupling parameters derived from C60 must be close to the parameters for C${}_{60}^{+}$ derived based on the definition.

The accuracy of LDA data have been discussed Ramanantoanina et al. (2013); Ponzellini (2014) based on the comparison between experimental photoelectron spectra (PES) Canton et al. (2002) and those from numerical simulation Manini, Gattari, and Tosatti (2003). Indeed, PES is very useful to establish vibronic coupling parameters, whereas the simulation of PES requires high accuracy both in theoretical simulation and experiments. In the case of C${}_{60}^{-}$, vibronic coupling parameters derived from broad PES at high-temperature Gunnarsson et al. (1995) have been proved to be overestimated by the simulation of high-resolution PES spectra Iwahara et al. (2010). Furthermore, it was also found that the error bar of vibronic coupling parameters derived from broad PES is very large Iwahara et al. (2010). Thus, the derivation of accurate coupling parameters is only possible if we have high-resolution PES spectra measured at low-temperature Wang, Woo, and Wang (2005); Huang et al. (2014). In the case of C${}_{60}^{+}$, as pointed out by Manini et al. Manini, Gattari, and Tosatti (2003), experimental PES is broad and fine structure of low-energy region due to vibronic coupling is completely smeared out, which prevents the direct comparison between theory and experiment. Moreover, in the case of PES of C${}_{60}^{+}$, theoretical ratio of the second strongest peak to the strongest one is smaller than those of experimental data, implying the underestimation of vibronic coupling by LDA. From this point, the present data larger than the LDA data by 40 % would give better agreement.

The quality of B3LYP calculations has been checked in C60 anions by comparing theoretical and experimental data in previous studies. Besides the good agreement between coupling parameters from B3LYP calculations and high-resolution PES Iwahara et al. (2010), the good quality of B3LYP calculations also has been confirmed from Néel temperature Iwahara and Chibotaru (2013), spin gap Liu, Iwahara, and Chibotaru (2018), and the explanation for the origin of temperature evolution of infrared spectra Matsuda et al. (2018) in Mott-insulating Cs3C60 using the same vibronic coupling parameters. Furthermore, vibronic coupling parameters from B3LYP calculations tend to give good description of inelastic electron tunneling spectra of other organic molecule Shizu, Sato, and Tanaka (2010). All these facts show that B3LYP values are closer to the reality in C${}_{60}^{+}$, but there still a mismatch when compared with experimental data. The application of CAM-B3LYP could eliminate such mismatch of vibronic parameters in fullerene system, as shown in recent study of C${}_{60}^{-}$Huang and Liu (2020).

Although the derivation of vibronic coupling parameters is the first step toward full description of the molecular states of C${}_{60}^{+}$, we believe this is a crucial step. Once calculations of accurate enough vibronic states become possible, it is possible to interpret various spectra such as scanning tunneling measurements of C60 Frederiksen et al. (2008), inverse PES Große et al. (2017) and angle resolved PES Latzke et al. (2019) to mention a few. Furthermore, present coupling parameters are derived based on the well-defined phase factor of normal modes which has been also used for orbital coupling parameters of LUMO Liu et al. (2018) and next LUMO Huang and Liu (2020). Therefore, by combining present coupling parameters with them, it is also possible to address complex vibronic problems of excited C60 Qiu, Chibotaru, and Ceulemans (2001), and also to analyze e.g. luminescence spectra Akimoto and Kan’no (2002) and relaxation process and thermally activated delayed luminescence Stepanov et al. (2002).

V Conclusions

In this work, orbital vibronic coupling parameters for $h_{u}$ HOMO level of C60 are derived using both B3LYP and CAM-B3LYP hybrid functional. We believe that these vibronic coupling parameters are high accurate and close to the real situation. With these obtained coupling parameters, JT stabilization energies of C${}_{60}^{+}$ are calculated, and JT structure at the minima of APES is confirmed to be $D_{5d}$, with the stablization energy 110 meV and 129 meV calculated with B3LYP and CAM-B3LYP, respectively. JT stabilization energies in C${}_{60}^{+}$ are about two times larger than that in C${}_{60}^{-}$, suggesting the crucial role of the dynamical JT effect to reveal the actual situation of C${}_{60}^{+}$.

Present coupling parameters have been derived within the same framework used for our studies on ground and excited C${}_{60}^{-}$. Thus, combining present data with that from other works, it is also possible to analyze the vibronic problems of excited C60.

acknowledgments

The authors thank Dr. Naoya Iwahara and Prof. Dr. Liviu Chibotaru for fruitful discussions. They also gratefully acknowledge funding by the China Scholarship Council (CSC). Dr. Dan. Liu is supported by ”the Fundamental Research Funds for the Central Universities”(G2019KY0517, G2019KY05104)

conflict of interest

There are no conflicts to declare.

VI appendix

VI.1 Clebsch-Gordan coefficients: $\hat{C}_{\Gamma\gamma}$

For the derivation of the vibronic Hamiltonian, the Clebsch-Gordan coefficients, $\hat{C}_{\Gamma\gamma}$ are taken from Ref. Liu et al. (2018), and listed below, in which $\phi=(1+\sqrt{5})/2$, $\Gamma$ = g, 1h, 2h, and $\gamma$ = a, x, y, z, $\theta$, $\epsilon$, $\xi$, $\eta$, $\zeta$.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

VII JT splitting of the HOMO levels

There are eight $q_{h\epsilon}$ and six $q_{ga}$ deformations, which are distinguished by the subindex, as $q_{h(1)\epsilon}$ corresponding to the first $q_{h\epsilon}$ deformation.

The DFT data with B3LYP hybrid functional and the defination of the phase factors of the normal modes are taken from Ref. Liu et al. (2018). The fitting for of the DFT HOMO levels to the model hamiltonian for $q_{h_{g}(i)\epsilon},~{}i=2,3,4,5,6,7,8$ are shown in Fig. 2, 3, 4, 5, 6, 7 and 8, while $q_{g_{g}(i)},~{}i=2,3,4,5,6$ are shown in Fig. 9, 10, 11, 12, and 13.

Considering CAM-B3LYP functional, Fig. 14, 15, 16, 17, 18, 19 and 20 depict the fitting of $q_{h_{g}(i)\epsilon},~{}i=2,3,4,5,6,7,8$, with Fig. 21, 22, 23, 24, and 25 for $q_{g_{g}(i)},~{}i=2,3,4,5,6$.