Electronic and vibrational linear and nonlinear polarizabilities of Li@C60 and [Li@C60]+

Electronic and vibrational nuclear relaxation (NR) contributions to the dipole (hyper)polarizabilities of the endohedral fullerene Li@C60 and its monovalent cation [Li@C60]+ are calculated at the (U)B3LYP level. Many results are new, while others differ significantly from those reported previously using more approximate methods. The properties are compared with those of the corresponding hypothetical noninteracting systems with a valence electron transferred from Li to the cage. Whereas the NR contribution to the static linear polarizabilities is small in comparison with the corresponding electronic property, the opposite is true for the static hyperpolarizabilities. A relatively small, but non‐negligible, NR contribution to the dc‐Pockels effect is obtained in the infinite frequency approximation. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011


Introduction
Since their discovery endohedral fullerenes have been extensively investigated because of their novel structure and properties. 1 Electrical properties have been of major interest owing to a variety of possible applications ranging from qubits for quantum computation 2 to organic photovoltaic devices. 3 Our interest in this connection lies in the theoretical determination of the linear and nonlinear optical properties, i.e., the (hyper)polarizabilities, of these materials. For that purpose we have chosen initially to study the prototypical metal endohedral fullerene, Li@C 60 , and its cation [Li@C 60 ] + .
It has been recognized for some time now 4 that large amplitude vibrational motions, due to weak interactions between the dopant Li atom and the fullerene cage, could give rise to large vibrational (hyper)polarizabilities. Thus, both the pure electronic and vibrational contributions need to be examined. So far both contributions have been treated only at a rudimentary level, to a large extent because of the large size of the fullerene cage and the open shell character of the neutral. Furthermore, only static properties have been considered as far as computations are concerned. The pure electronic first hyperpolarizability 5 and second hyperpolarizability 6 of the neutral have been calculated using an uncoupled Hartree-Fock scheme with molecular orbitals obtained from a restricted open-shell Hartree-Fock calculation. These studies also contain experimental measurements. Very recently Yaghobi et al. 7 used a modified Su-Schrieffer-Heeger (Huckel-type) Hamiltonian, 8 coupled with a sum-over-states procedure, to obtain the Li@C 60 electronic linear polarizability and second hyperpolarizability tensors.
Whereas all the theoretical studies mentioned here discuss the vibrational contribution to the (hyper)polarizabilities, only Whitehouse and Buckingham 4 made an attempt to calculate these properties for the type of system in which we are interested. They used a much simplified potential, in conjunction with a classical analysis -estimated to be valid above 20 K -to obtain a (temperature-dependent) expression for the vibrationally averaged dipole moment. Then, from the field-dependence of this expression, formulas for the vibrational linear polarizability and second hyperpolarizability were extracted and the former quantity was evaluated for the [Li@C 60 ] + cation. Their work indicated that the vibrational contribution could be many times larger than the electronic contribution for this system.
As concerns endohedral fullerenes in general, some calculations of their polarizabilities have been reported, usually employing DFT methods, 9 but reports on vibrational polarizabilities are rare. We mention here the work of Pederson et al., 10 who computed the vibrational polarizability of Kr@C 60 , as well as of C 60 itself, in the lowest order of perturbation theory (double-harmonic approximation). At that level, these contributions were found to be very small compared with the electronic property.
There have been many important advances in computational capability and theoretical methodolgy since the articles on Li doped C 60 noted above have appeared. While these articles establish a clear interest in the linear and nonlinear optical properties of the prototypical neutral and cation, the time is ripe for a significantly improved treatment. That is the aim of our current study.

Methods
The geometry of Li@C 60 and the singly charged cation were optimized at the DFT level using the B3LYP functional in the unrestricted version for the former and the restricted version for the latter. We note here that spin contamination was not an issue in the unrestricted DFT calculations, neither with nor without applied electric fields. In all cases, the expectation value of the total spin squared operator, S 2 for the converged Li@C 60 wavefunctions was found to be between 0.755 and 0.756, which is close to the exact value of 0.75 for a pure doublett. Due to the large size of the molecule, as well as the tight convergence requirement for calculating vibrational (hyper)polarizabilities, the rather small 6-31G basis set was employed in the latter calculations (some 6-31+G values are included for comparison). For the electronic properties, the 6-31G basis was found to be inadequate, thus they were additionally computed with the 6-31+G and 6-31+G* basis sets. In addition, several geometry optimisations were done with 6-31+G and 6-311G*, although no Hessians could be computed with these large basis sets.
The static electronic polarizabilites [α e αβ ≡ α e αβ (0; 0)], first hyperpolarizabilities [β e αβγ ≡ β e αβγ (0; 0, 0)], and second hyperpolarizabilities (γ e αβγ δ ≡ γ e αβγ δ (0; 0, 0, 0)) are defined by the Taylor series expansions for the dipole moments µ α (F) 11 or energies E(F), 11 in terms of the static field F: From these computed properties, i.e., the dipole moment and energy, the (hyper) polarizabilites were obtained using the Romberg differentiation procedure. Convergence difficulties in the DFT self-consistent orbital calculations severely limit the accuracy of the numerical differentiation as well as the range of fields that can be used. The problem increases for larger fields and thereby impacts the accuracy of the (hyper) polariabilities. Thus, the second hyperpolarizability could be obtained with sufficient statistical confidence from the dipole moments, but not from the energies. Generally, a third-order Romberg differentiation 12 with a minimal applied field strength of 0.001 au was applied. In some cases, only a second order Romberg differentiation was possible due to convergence problems, but in some other cases even a higher order was used to ensure the reliability of the values. In addition to the DFT computations, selected Hartree-Fock (HF) and second order Møller-Plesset (MP2) computations were performed with the 6-31G basis set to provide a comparison with traditional wavefunction theory methods. The computations were done with Gaussian03 13 and Gaussian09. 14 Some of the problems mentioned earlier could have been avoided by using an analytical derivative method, e.g., analytical response theory, as implemented for linear, quadratic and cubic response functions in time-dependent DFT by Jansik et al. 15 in the program package Dalton. 16 This method has also been extended using spin-restricted DFT for open-shell systems. 17,18 However, trials to compute the (hyper) polarizabilities of Li@C 60 using the smallest basis set (6-31G) were successful only for linear polarizabilities; convergence problems in determining the response vectors prevented calculation of the hyperpolarizabilities. Considering that these calculations were quite expensive, and finite fields were needed for the nuclear relaxation treatment of vibrational contributions, we decided to employ finite field techniques throughout. We mention, finally, that the polarizabilities obtained from spinrestricted analytical response theory were nearly identical to those calculated by unrestricted finite field methods, showing again that spin-contamination is not a problem in our case.
Although there has been a lot of progress in the last few years in the field of computing vibrational (hyper) polarizabilities based on vibrational seld-consistent field theory and correlated versions thereof (see e.g., refs. 19,20), the corresponding methods are computationally still much to expensive for the large systems of interest here. Thus, the vibrational contributions were mostly computed using the finite field approach pioneered by Bishop et al., 21 and later implemented by Luis et al. 22 In this approach, the molecular geometry is first optimized in the presence of a static electric field while strictly maintaining the Eckart conditions. 22 Then the difference in the static electronic properties due to the change in geometry induced by the field is expanded as a power series in the field. Each term in the expansion yields the sum of a static electronic (hyper) polarizability plus a nuclear relaxation (NR) vibrational term. For example, the change of the dipole moment and of the linear polarizability are given by: 21 with The argument R F implies structure relaxation in the field, and P nr means the nuclear relaxation part of P, while the subscript ω → ∞ invokes the so-called "infinite optical frequency (IOF)" approximation. In principle, this procedure allows one to obtain most of the major dynamic vibrational NR contributions in addition to the purely static ones of eqs. 5-7. The linear term in the electric field expansion of eq. (4) gives the dc-Pockels effect; the quadratic term gives the optical Kerr Effect; and the linear term in the expansion of beta yields dc-second harmonic generation (all in the IOF approximation). For laser frequencies in the optical region it has been demonstrated that the latter approximation is normally quite accurate. [23][24][25] In fact, this approximation is equivalent to neglecting terms of the order (ω v /ω) 2 with respect to unity (ω v is a vibrational frequency). In terms of Bishop and Kirtman perturbation theory 26-28 all vibrational contributions through first-order in mechanical and/or electrical anharmonicity, and some of second-order, are included in the NR treatment. 29 The remaining (higher-order) vibrational contributions can, in principle be computed as well using a related formulation. 30 However, that treatment requires computation of the field-dependent zero-point vibrationally averaged properties, which was not feasible for the systems studied here because of their large size and complicated potential energy surface (PES). Indeed, of the dynamic properties mentioned earlier, we were only able to obtain the dc-Pockels Effect due to instabilities for high fields that will be described later.
Generally, field strengths from 0.0001 to 0.0128 au were tried in the Eckart-constrained optimizations and the energies and dipole moments of the successfully optimized structures were subjected to a numerical Romberg differentiation. As in the case of the electronic properties, the numerical differentiation of the energies was too unstable to yield all the properties of interest. However, the numerical differentiation of field-induced dipole moments allowed us to obtain stable values for most of the components of the NR contribution to the static α , β, and γ and to the IOF approximation for the dc-Pockels first hyperpolarizability. Due to the high computational cost of these calculations, the 6-31G basis set had to be used. A few control calculations with the 6-31+G basis set showed that the influence of diffuse basis functions on the vibrational properties is not negligible, but smaller than on the electronic properties.

Geometry Optimization
Zhang et al. 31 have calculated the UB3LYP/6-311G* potential energy surface (PES) for motion of Li along five different rays passing through the center of an undistorted fullerene cage in Li@C 60 . The two most important rays, as far as the structure is concerned, were along the line from the cage center to the center of a C 6 hexagon (symmetry C 3v ) and along the line from the cage center to the center of a C 5 pentagon (symmetry C 5v ). Since localization due to cage distortion can be important, as they found, we carried out geometry optimizations for near C 3v and near C 5v symmetry at the same level while allowing the cage to fully relax. We also obtained the stationary point at the near icosahedral symmetry under the same conditions. In all cases, it was necessary to lower the actual symmetry to obtain the optimized structure, due to SCF convergence problems in the high symmetry calculations. Zhang et al. 31 report that similar difficulties ocurred in their calculations.
In the case of approximate I h symmetry the cage was slightly distorted yielding a C s optimized structure with the Li atom slightly (0.015 Å for UB3LYP/6-31G) removed from the cage center. As already well-known this stationary point is not a minimum; in fact, there are four imaginary frequencies. For the two minima (near C 5v and near C 3v ) the Li atom shifts about 0.1 Å off the ray that goes from the center of the cage to the center of the polygon and the symmetry is again reduced to C s . In both instances, the Li atom was located at about 1.5 Å from the center of the cage in the optimized structure. The eccentric position of Li in Li@C 60 has been interpreted in terms of dispersion and repulsion 32 interactions.
All of the above results agree semiquantitatively with Zhang et al. 31 as expected. The structure with approximate C 3v symmetry was found to be 3.4 kJ/mol (1.6 kJ/mol) more stable than the one with C 5v symmetry using UB3LYP/6-311G* (UB3LYP/6-31G). Our value is slightly higher than the one found by Zhang et al.
(2.6 kJ/mol). The reason for this difference may be due to cage relaxation and/or small deviations of the Li atom from the fixed ray they employed. It is also posible that the two minima are not directly related; as shown later, there seem to be several minima close by. The energy difference between the near-Ih and near C 3v symmetry structures was found to be 56.5 kJ/mol at the 6-31G/UB3LYP level. Since we are interested in the ground vibrational state all further investigations were focused on the most stable near-C 3v structure. The optimization using the 6-31+G basis was started from the 6-31G optimized structure, and the final optimized structure was very close to the starting one, as expected.
The geometry of the monovalent cation [Li@C 60 ] + was also determined for the near C 3v symmetry, using restricted B3LYP and the 6-31G basis set. A control optimization using B3LYP/6-311G* did not show any substantial structural differences. At the minimum, the Li atom is about 1.4 Å off the center of mass of the cage. The average C-C bond-length (1.4406 Å) is nearly the same as that of the neutral (1.4412 Å). However, the cation is somewhat more spherical. As a measure of the sphericality we use where I x is the principal component of the cage inertia tensor, with respect to the center of mass, in the xdirection. Our values are I = 1.4 g Å 2 /mol for the cation and I = 31.5 g Å 2 /mol for the neutral. For comparison, the average moment of inertia I = 1/3(I a + I b + I c ) is about 3050 g Å 2 /mol for both species. Finally, the coordinates of both near C 3v optimized structures can be obtained from the authors.

Electronic Properties
To assess the reliability of the level of theory chosen, we show in Table 1   The rather small dipole moment of the neutral depends strongly on the basis set and correlation treatment. Fortunately our interest lies in the (hyper) polarizabilities. Nonetheless, we can say that our dipole moment results are in qualitative agreement with the value computed by Campbell et al., 5 but much smaller than reported by Li and Tomanek. 33 The comparison of our calculated (hyper) polarizabilities with those of Campbell et al. will be made later. We note that the addition of diffuse functions to the 6-31G basis is always crucial. But whereas the further addition of polarization functions has a minor effect on α and γ , for β they partially (or totally) offset the effect of the diffuse functions. Comparison of HF and (U)MP2 shows that correlation has a very large effect on β, but not α The fact that (U)B3LYP yields values similar to (U)MP2 suggests that both account fairly well for correlation (even though the calculations are only at the 6-31G level). Due to the unavailability of γ at the (U)MP2 level, no conclusions can be drawn in this respect for the second hyperpolarizability. The cation properties are somewhat less sensitive to correlation than those of the neutral. Overall, we conclude that (U)B3LYP/6-31+G* is the minimal level required to obtain reliable electronic properties. Finally, the large change in hyperpolarizabilities upon going from the cation to the neutral is not unexpected in view of the additional electron in a (formerly) unoccupied orbital localized on the C 60 moiety. For γ the effect is more evident for the other two diagonal components shown in Table  2 (see later).
In Table 2 we show the computed static (U)B3LYP/6-31+G* electronic properties of Li@C 60 , [Li@C 60 ] + , C − 60 , C 60 , and Li. The several additional species were included for comparison of Li@C 60 with both noninteracting Li + C 60 and Li + + C − 60 , as well as for  34 in the current context and are not taken into account in the discussion later. For our comparisons the geometry of the endohedral fullerenes was optimized at the (U)B3LYP/6-31G level and the same geometry was retained for the noninteracting species. In principle, a BSSE correction should be applied to the Li-doped fullerenes, but is omitted since it would have no effect on our conclusions. As compared to the hypothetical noninteracting species the interaction between the Li + cation and either the C − 60 or C 60 cage leads to a moderate reduction of the diagonal polarizabilities and, in the case of [LiC 60 ] + also of the second hyperpolarizabilities. The reduction for α may be due to a contraction of the electron density caused by the attraction of the cation. Such an explanation will not suffice for γ since the effect of the interaction on the diagonal components is quite different for Li@C 60 , and the second hyperpolarizability is not simply related to the size of the electron distribtion. The first hyperpolarizabilities arise from asymmetry of the charge distribution and are, consequently, strongly enhanced in the endohedral species.
Campbell et al. 6 used an uncoupled approximation to coupledperturbed HF theory -or, as they prefer to call it, a "computationally expensive tight-binding approach" to compute the hyperpolarizabilities of Li@C 60 , using the 6-31G* basis set, and obtained for γ the values ∼(320, 540, −320)×10 3 au, for the x, y, and z diagonal components, respectively. While the x value is quite close to ours, the other two values do not agree even in sign. In ref. 5,Campbell et al. also computed the first hyperpolarizability using the same methodology. The values they obtained (β zzz ∼ 15,000 au, β yyy approximately −7000 au) are at least one order of magnitude larger than ours in the most similar geometry they considered (Li displaced about 1.5 Å from the center towards an hexagon). Campbells' values are based on orbitals obtained from a ROHF calculation, while ours are computed at the UB3LYP level. The large differences between the two results confirm the unreliability of the HF method for the hyperpolarizabilities of Li@C 60 , as found here for the UHF values (cf. Table 1). In the same approximation, but now using RHF, they obtain ∼50.000 au for the diagonal component in C 60 . 6 This value is about 2.5 times smaller than ours, which may be mostly due to the different basis sets, but also due to different approximations in both approaches, as well as differences in geometry etc. Finally, we note that Jansik et al. 15 computed values for the (hyper) polarizabilities of C 60 in I h symmetry with analytic response theory using larger basis sets, specifically tailored for the purpose of computing hyperpolarizabilities, and obtained α av = 547.0 au and γ av = 118 × 10 3 au with B3LYP using the cc-pVDZ+spd (their notation) basis set. These values are quite comparable to ours, taking into account the differences in symmetry (I h versus C s ), geometry, and basis set (it was assumed in our case that the geometry of C 60 is sufficiently spherical so that the nondiagonal terms of γ do not deviate appreciably from the relation γ iijj = 1/3γ iiii ).
We also investigated the influence of the position of the Li atom along the dipole axis on the electric properties of Li@C 60 . The  Table 3.
For reference purposes the surface of the cage is at a distance of ∼3.4 Å. Interestingly, the diagonal z-component of the polarizability, first hyperpolarizability, and dipole moment change little (µ, α) or moderately (β) (for µ we are assuming that the value of the change is reliable even though the value of the dipole moment itself is not) between r = 0.729Å and r = 2.0 Å, but γ zzzz is altered much more drastically, even undergoing a sign change at small distances from the cage center. This is consistent with γ being due to electron density that is distant from the surface of the cage. The large gradient in γ zzzz for Li@C 60 could possibly be used in a potential nonlinear "flip-flop" device. This would require a mechanism such as an STM electric field to shift the equilibrium position of the Li atom between different regions. The magnitude of such a shift has been investigated by Delaney and Greer 2 who found that it is difficult to move the Li atom very far because of the large screening effect of the fullerene cage. In the calculations reported below we find that a shift from the field-free position of about z ∼ 0.03 au will result when a 0.0128 au field is applied. According to Table 3, this shift is much too small to change γ zzzz appreciably.

Nuclear Relaxation Contribution to Vibrational nlo Properties
In this subsection, we present the nuclear relaxation (NR) contributions to the vibrational (hyper) polarizabilities of Li@C 60 and [Li@C 60 ] + . As mentioned in Sec. II our treatment requires a geometry optimization in the presence of a finite field. A problem can arise when there are multiple minima on the PES separated by low energy barriers. The finite field method works satisfactorily in that event as long as the field-dependent optimized structure corresponds to the same minimum as the field-free optimized structure. This was the case in previous work on ammonia, 35 which has a double minimum potential. However, it is sometimes not the case for the endohedral fullerenes considered here, especially Li@C 60 . In fact, we were unable to determine the NR contribution in the x direction, i. e., perpendicular to the symmetry plane, for that molecule. It was possible to obtain α nr xx , based on the alternative analytical formulation, 26-28 utilizing field-free dipole (first) derivatives and the Hessian. The analytical polarizability components in the other two directions were, In addition to the situation just discussed, it was also found that the electric field can sometimes lead to a change of electronic state, as detected by a sudden jump in the computed polarizability. This further limited the range of applicable field strengths and, thus, the range of properties that could be computed with sufficient statistical confidence.
Our results for the static NR contributions to the diagonal (hyper) polarizability components are shown in Table 4. Most of the values were obtained at the (U)B3LYP/6-31G level. For comparison, a few calculations were also done for Li@C 60 at the UB3LYP/6-31+G level. As seen from the Table, diffuse basis functions have a non-negligible effect on the computed values, although the effect is smaller than on the electronic properties (cf. Table 1). Note in particular that there is no qualitative change of any vibrational property upon going from 6-31G to 6-31+G, in contrast to the electronic properties where 6-31G gives a negative value for γ zzzz , while it is positive for 6-31+G. Thus, we expect that the values of vibrational properties obtained with the 6-31G basis are qualitatively correct, although the accuracy becomes worse for γ nr than it is for α nr or β nr .
For α the vibrational contributions are quite small in comparison with their electronic counterparts (cf. Table 2). In the case of [Li@C 60 ] + this appears, at first glance, to contrast with what Whitehouse and Buckingham (WB) 4 have previously found. However, their values were obtained by classical averaging in the high temperature limit -in this case above 20 K -while ours are for 0 K. Another difference is that WB obtain the complete vibrational polarizability and second hyperpolarizability, albeit very approximately, whereas we have not included the so-called curvature contribution. 30 For the polarizability one would "normally" expect the latter to be substantially smaller than the NR term, but endohedral fullerenes are not "normal" molecules and that may not be the case here. Of the several approximations in the WB treatment, one of the most questionable is the spherical approximation for the field-free potential, which Zhang et al., have shown does not qualitatively reproduce the low energy vibrational spectrum of the neutral. This is not to mention that WB considered a rigid cage and motion along the C 5v symmetry axis, rather than the lower energy C 3v symmetry axis.
In contrast to α, the NR contributions to β and γ , are quite large. For the cation, in particular, two diagonal components of γ nr are larger than the corresponding electronic components. Such an increase in the relative magnitude of the NR (hyper) polarizabilities, as compared to the electronic values, as the order of nonlinearity increases is often observed in conjugated systems. 25,36 However, it is usually not to as large an extent as found here. The relationship with degree of nonlinearity may be connected with the fact that only dipole derivatives enter into the linear NR polarizabilities, whereas β nr and γ nr depend additionally on polarizability derivatives, as the perturbation expressions for these quantities show. 28 Because of the conjugation the polarizability derivatives tend to be large (small changes in bond length alternation cause large changes in α). Furthermore, the higher-order vibrational polarizabilities (as opposed to the linear vibrational polarizability) depend upon electrical and mechanical anharmonicity which, undoubtedly, is especially large for the systems we are considering. WB give an expression for γ , but no numerical values. In addition, they considered only linear terms in their field-dependent vibrational Hamiltonian. This omits contributions to the vibrational hyperpolarizability that are generally important as suggested previously.
In Table 5 we show the computed NR contribution to several components of the dc-Pockels effect for Li@C 60 , i.e., β nr ijy (−ω; ω, 0) and β nr ijz (−ω; ω, 0) where ij = xx, yy, yz, and zz. [see eq. (8)]. Comparison with the corresponding electronic values of Table 2 shows that the vibrational contribution is relatively small, compared to the corresponding static electronic component, but not negligible. The values obtained here may be compared with those of a typical donor-acceptor molecule, H 2 N-(CH=CH) 3 -NO 2 , for which the ratio β nr iii (−ω; ω, 0) ω→∞ /β e iii (0; 0, 0) along the dipole direction was found to be about 0.7 at the RHF level and about 0.2 at the MP2 level. 36 For LiC 60 this ratio is 0.13 at the UB3LYP level.

Conclusions
We have computed both electronic and NR vibrational contributions to the (hyper) polarizabilities of the prototype endohedral fullerene Li@C 60 and its cation. A number of these properties were obtained for the first time. In other cases, our results differ quite signicantly from those previously determined using more approximate approaches. The latter include the static electronic properties calculated by Campbell et al. 5,6 Although, for the cation, there is a large difference between our values of the static vibrational contribution to α and those reported by Whitehouse and Buckingham, 4 these results are not really comparable because their calculations include the effect of temperature. In addition, they applied several strong approximations, such as assuming a spherical field-free potential inside the cage. On the other hand, our calculations do not include higher-order vibrational contributions omitted in the NR treatment. It would be worthwhile to add temperature-dependence to the NR approach as we plan to do in the future. Whereas the NR contribution to the static α is quite small for both endohedral fullerenes, it becomes quite large for the static hyperpolarizabilities. This contribution is reduced for the dynamic Pockels effect, computed in the infinite optical frequency approximation, but is still not negligible. For [Li@C 60 ] + the calculated (hyper) polarizabilities are roughly comparable to those of the hypothetical noninteracting system obtained by charge transfer of the Li valence electron to the cage giving Li + + C 60 . The same is true for the linear polarizability of the neutral but the noninteracting charge transfer model completely breaks down for the hyperpolarizabilities.
We consider our work as a substantial step towards the final goal of a full computational characterization of the linear and nonlinear electric dipole properties of endohedral fullerenes. As far as vibrational contributions are concerned, in addition to the NR treatment of low-order perturbation terms, there is an established method for obtaining all remaining contributions through calculation of zero-point vibrationally averaged properties at the relaxed field-dependent geometry. 30 What is needed is a robust procedure for carrying out the geometry optimizations when the PES has multiple minima and/or other strongly anharmonic features. Work is in progress on a reduced dimensionality scheme that may be combined with quasidegenerate perturbation theory to treat such circumstances. 35