General static polarizability in spherical neutral metal clusters and fullerenes
within Thomas-Fermi theory

Clusters are by definition, mesoscopic systems formed from $3-10^{7}$ reinhard2008introduction atoms of various elements. The theoretical methods of investigation are going from molecular to bulk (solid-state) domain, both quantum and classical approaches.

As we have mentioned before, the clusters formed from metallic elements have the property that the electrons from the unclosed shell are loosely bounded and their behavior is close to HEG one. Also, we treat the problem in the Born-Oppenheimer approximation. This means that their dynamics and de-localization are high enough to consider that they $see$ the ionic background in a averaged manner. By $ionic$ $background$ we understand the system obtained from the coupling between charged nuclei and the $core$ electrons which determine a net positive charge.

Our working framework is the so called $jellium$ $model$ in which the ionic background is approximated to a homogeneous positive charge distribution while the free electrons will also have an almost constant density in the bulk region. In clusters the jellium model it is also applied and the ionic background determines a smooth Coulomb potential with the appropriate symmetries. Once the Coulomb potential is generated it may be included as an input in different theoretical approaches (Hartree-Fock, DFT, etc.) to derive the electron density. The self-consistent jellium model proved to be a very appropriate method predicting quantitative results in good agreement with the experimental data brack1993physics , de1993physics . The simplest geometry possible is that of a sphere, but can be somehow a troublemaker in numerical simulations since it has a discontinuity at the edge. A more refined model is the so called $soft$ $jellium$ which works with a Wood-Saxon radial profile coupled with any possible angular dependence, giving access to any possible geometry of the ionic part. The jellium density is written in its most general form as:

with $R(\theta,\phi)=R_{0}(1+\sum_{l,m}\alpha_{lm}Y_{lm}(\theta,\phi))$. The Weitzecker-Seitzs radius $r_{s}$ is a parameter of the bulk domain interpreted as the volume occupied by a single atom while $Z$ is the difference between the numbers of protons and the number of bound electrons and so $Z|e|$ is the charge of jellium. Therefore, the jellium density satisfies the condition: $\int\rho_{jel}dr^{3}=Z$. The limiting case of spherical sharp distribution is obtained when $\sigma_{jel}\to 0$.

Thomas-Fermi model (TF) was derived independently by L. Thomas and E. Fermi in 1927, soon after Schroedinger equation, 1926. The basic approximation of the model is to treat the electron density distribution in the atom using a local approximation for the kinetic term i.e. the HEG approximation: $E_{kin}=\frac{3\gamma}{5}\int_{\mathbb{R}^{3}}\rho^{5/3}(\vec{r})\mathrm{d}r^{3}$ , with $\gamma=(3\pi^{2})^{2/3}\hbar^{2}/2m$. If $V(\vec{r})$ is the external potential , the total energy density functional can be written:

The ground-state energy and the corresponding density distribution are obtained within a Ritz variational principle, searching for the minimum of this quantity:

The condition $\rho\in\mathfrak{L}^{5/3}(\mathbb{R}^{3})$ refers to the fact that density is a $5/3$ integrable function over $\mathbb{R}^{3}$ and so, is a condition for finite kinetic energy to emerge. The constrain associated with the condition $\int_{\mathbb{R}^{3}}\rho(\vec{r})\mathrm{d}r^{3}=N$ is introduced through Euler-Lagrange multiplier technique and solving the variational problem, the resulting Thomas-Fermi equation is be derived:

In this equation, $\Phi(\vec{r})$, $\Phi:\mathbb{R}^{3}\to\mathbb{R}$ is the energetic Coulomb potential while $\mu$ is the chemical potential:

In turn, the Coulomb potential is the generated by the charge distribution and is the solution of the Poisson’s equation $\Delta\Phi(\vec{r})=|e|\Delta V(\vec{r})+e^{2}\rho(\vec{r})/\varepsilon_{0}$. In the case of neutral electric systems, the chemical potential is null lieb1977thomas , but this feature is maintained only within the Thomas-Fermi method. The additional terms, as Dirac or Wietzacker contributions, break this property. In the jellium approximation $V(\vec{r})$ is induced by $\rho_{jel}$ and the Thomas-Fermi becomes in differential form:

In the absence of external interactions, for the ground-state, the TF equation leads to $\gamma\rho_{0}^{2/3}(r)=\Phi_{0}(r)$ with $\rho_{0}$ the ground state density of the free electrons. We shall study in the following the static linear response in TF approximation considering a time-independent potential of an arbitrary form. If the coupling strength to the free electrons is $\lambda$:

And the TF equation for the new stationary state, becomes:

Then the interaction energy density is $\rho v(\vec{r})$. This perturbation will induce a spatial change of the charge which can be treated in a power expansion in the coupling strength $\lambda$: $\rho(\vec{r})=\rho_{0}(r)+\lambda\rho^{1}(\vec{r})+\lambda^{2}\rho^{2}(\vec{r})+...$.

If we resume to first order term $(\lambda\ll e^{2}/(4\pi\varepsilon r_{0}))$, then $\rho^{1}$ should satisfy $\int_{\mathbb{R}^{3}}\rho^{1}(\vec{r})\mathrm{d}r^{3}=0$ since $\int_{\mathbb{R}^{3}}\rho_{0}(\vec{r})\mathrm{d}r^{3}=\int_{\mathbb{R}^{3}}\rho(\vec{r})\mathrm{d}r^{3}=Z$.

In the presence of the external potential the initial spherical symmetry is broken and the densities varify the properties:

The linearized kinetic density energy term and the potential $\Phi$ are:

Then the TF equation for the perturbed part becomes:

Working in spherical coordinates we can consider the following expansion of $\rho^{1}$

Using the expansions for potential $\eqref{potexpansion}$ and for density $\eqref{densexpansion}$ in equation $\eqref{eq:tfro1}$ we deduce that $u_{lm}(r)$ functions will satisfy the radial equation:

Our equation is an approximate version of the equation (23) deduced in serra1989static in the absence of Weizsacker, exchange and correlation terms. This fact can be seen by setting the $beta$ factor to $0$ and excluding exchange-correlation effects. Nonetheless, their derivation is done on the basis of variational method, while our is not. We stress in this paper the fact that the complexity of the equation proposed in serra1989static is much more involved being an integro-differential equation and the presence of those supplementary terms is not necessary for semi-quantitative results. In fact, the results have the same level of accuracy with ours and since the present equation is pure differential we can take advantage of its simplicity in order to study the main effects which contribute to polarizability in metal clusters. Moreover, while in the reference, the Weiszacker correction is used with different constants and is considered to be the essential reason for which the differential equation is derived, our derivation has no such condition.

Concerning the boundary conditions, in the origin, in order to have finite $\rho^{1}(0)$ we require that $\lim\limits_{r\to 0}u_{lm}(r)/r$ to be finite. Concerning the behavior at infinity we shall ask for the coefficient $u_{lm}(r)$ to follow the behavior of the perturbation i.e. $u_{lm}(r)=3/(2\gamma)v_{lm}$ when $r\to\infty$. From TF eq $\eqref{eq:tfro1}$, the asymptotic behavior of $u_{lm}$ is:

Here, the $q_{lm}$ term is the multipole moment associated with the induced charge:

From numerical point of view we solve the equation $\eqref{eq:tfro1}$ with the associated boundary conditions as it follows: first we guess the term $q_{lm}$ (considering the particular system to be studied, the magnitude can be easily guessed) and solve the equation in such a way that the solution satisfies the asymptotic behavior mentioned above for the selected value of $q_{lm}$. With the solution constructed in these way we find a new value of $q_{lm}$ and repeat the procedure until we reach convergence condition of the solution.

Even though can seem to be a long road to the convergence, in practice, this is reached within 10 iterations.

Again, in comparison with serra1989static we have different asymptotic boundary condition, simply from the form of our equation and the meaning of the unknown. The entire method of iteration for finding the polarizability as a parameter of asymptotic behavior is original, at the best of our knowledge and essential for the results. The above mentioned reference does not discuss such matters.

In this section we shall apply the method described above to the specific case of dipolar response. The applied field is $v(\vec{r})=rcos\theta$ and consequently $v_{10}=r^{2}$. The equation $\eqref{eq:ultimate}$ becomes:

Of course, this can and must be solved numerically, but there are some specific cases with analytic solutions which can be helpful for the boundary conditions.

For $\rho_{0}(r)=0$, which usually describes the regions with large $r$ where the density must be null we have the solution which respects the boundary condition (19) for $q_{10}$, $u_{10}(r)=Ar^{2}+\frac{B}{r}$. In the more general case of $\rho_{0}(r)=const$ we have a more elaborate solution ($k^{2}=6\pi\rho_{0}^{1/3}/\gamma$):

Sodium clusters represent a textbook metal cluster due to the nature of the element which has a single electron on the last shell, namely the $Na$ with the atomic number $Z=11$ and the electronic configuration $1s^{2}2s^{2}2p^{6}3s^{1}$. This element has been taken into account in this work, due to their close to spherical symmetry for medium sized clusters and due to the strong metallic character. The classical electromagnetism provides, in the frame of small metal sphere model, a polarizability connected with the radius by:

Experimental data reveals higher poralizabilities for all $Na$ clusters, only in the high radius limit, the classical value is reached.

In our calculations, different clusters were taken into account as having spherical symmetry and a constant density of atoms. The electrons on the first two atomic levels were considered as core electrons and so, the jellium model reduces to a sphere of a radius connected to the number of atoms through the Weitzecker-Seitzs radius. The positive jellium charge of $N|e|$ distributed by a Wood-Saxon profile like in $\eqref{jellium}$ (but no angular dependence), with a sharp fall of density around the radius of the cluster ($\sigma\simeq 0.8a_{0}$ usually used reinhard2008introduction ) and Seitz radius $r_{0}=3.93a_{0}$, atomic units of length, see Fig. 1a), for the case $N=40$. Other parametrization of smaller $\sigma$ have been explored, but due to sensitivity of the method far from the center of the cluster, this parametrizations give worse results.

Taking into account the spherical symmetry, the TF equation reduces to radial differential equation:

For the same number of atoms, in Fig  1a) (continuous line) is plotted the electron density $\rho$ as obtained from above equation $\eqref{TFNa}$. In practice, clusters with the number of atoms between $20$ and $100$ were investigated, but in Fig.  1 just a generic plot of the jellium density and the electron density is presented in units of $1/r_{0}^{3}$, with $r_{0}$ the radius of the cluster.

This distribution manifests a tail beyond the jellium volume associated with the quantum behavior of the electrons. The ground state electrostatic potential $\Phi_{0}$ is plotted in Fig.  1b), while the radial dependence of the induced charge density is represented in Fig  1c) (blue-filled). In figure Fig.  1c) we’re drawn for comparison the asymptotic function as described by $\eqref{eq:cond1}$ and actual solution.

The proposed equation $\eqref{eq:ultimate}$ has been used in the dipolar case with the analytic limits $\eqref{eq:cond1}$,$\eqref{eq:cond2}$ and the polarizability was obtain in a quite good agreement with the experimental resultstikhonov2001measurement . The results can be seen and compared with reference tikhonov2001measurement in tabel  1 and in  2 where the quantity $\alpha/N$ as a function of $N$ is plotted. The solid horizontal line is associated with the classic solution.

For low numbers of atoms, the equation fail to describe quantitatively the polarizability due to the fact that the jellium model and the spherical shape of the cluster are no longer realistic approximations while in the range $N=40$, $N=100$ we can see an error bellow $10\%$ (plotted in  2), the only deviations being a consequence of shell effects. Also, the tendency of decreasing to a constant (bulk) value with the number of atoms involved, explained as a classical limit of our semi-classical treatment of the electron system.

From numerical solution obtained with the above method, we have observed that it is possible an empirical parametrization of the approximative $u_{10}(r)$ solution in the general case as:

Which allows us to formulate a final empirical sum rule-like expression for static dipole polarizability. The potential of this observation is that links the linear response only to the ground state properties of the system as in the usual moments of the response function:

We make comparison with the well known sum rule for static dipole polarizability [de1993physics ] exhibited by spherical metal clusters in which the main contribution is given through the so called $spilled-out$ electron which are considered outside the jellium region and appear proportional to $\delta$ in the approximation : $\alpha\simeq r_{0}^{3}(1+\delta)$ snider1983density .

The Buckminster fullerene represents one of the most studied molecule in the last decades due to its high symmetry, special features, high stability, etc. Consequently, the polarizability has been studied gueorguiev2004quantum in many models, the best theoretical results being obtained in the frame of DFT-LCAO while other theories obtained large errors in respect with the experimental data.

Essentially, fullerene is a carbon molecule with the atoms placed on a structure similar with that of a soccer ball. Due to the fact that carbon is not a genuine metal, one could argue that to study it along with true metallic clusters like sodium, it is a bad mistake. Nonetheless, even if from electronic structure and band gap point of view our approach can not be justified, it is a known fact that the optical spectra from fullerene exhibit a large, well localized plasmon. Further, this plasmon it is explained as being a surface plasmon scully2005photoexcitation and so, it can be reasonably concluded that in the dynamic regime, the electrons from fullerene behave close to the ones from a metal. For this reason we have applied Thomas-Fermi equation for the semi-delocalized electrons and obtain a good description for polarizability, Mie’s plasmon centroid and the density distribution for the ground-state.

We have used as jellium model, a gaussian distribution centered on the radius on which the carbon nuclei are placed but with a small width, described by the equation $\rho_{jel}(r)\propto exp(-\sigma(r-r_{0})^{2})$. Regarding the charge, our jellium model contain the core electrons from the $1s^{2}2s^{2}$ while the other 240 electrons from $2p^{2}$ are considered quasi-free and taken into account in the TF calculations. The results are quite sensitive to the width ($\Delta=$ full width at half maximum) of the jellium gaussian distribution and for that reason we have performed our calculations with different values for this quantity between $0.01\mathring{A}$ and $0.6\mathring{A}$. This impediment is hard to be avoided since the physical meaning of this width is the radius of the core electrons in which they can be accounted as part of jellium, but in an approximative physical way, our values cover the usual atomic values for this feature and provide a set of very close values to the experimental data.

Also, another approach to the jellium model was considered the spherical homogeneous shell (discussed in almost all papers using jellium model in $C_{60}$ investigations rudel2002imaging puska1993photoabsorption weaver1991electronic ) centered on the mean radius of $r_{0}=3.54\mathring{A}$ and with a width of $1.5\mathring{A}$rudel2002imaging .

The results in the electron density are physical and in good agreement with the experimental values for the inner, $\simeq 1.8\mathring{A}$ and outer radius of the fullerene $\simeq 5.1\mathring{A}$, see  5a),c). From the calculations of dipolar polarizability, we have obtained the expected volumic shift in density on the direction of the potential gradient ($Oz$ axes) shown in Fig  5b),d) while the value for polarizability is sweepings an interval between $80\mathring{A}^{3}$ and $85\mathring{A}^{3}$, depending on the chosen full width at half maximum of the gaussian jellium.

In Fig  6 we have plotted the results for $C_{60}$ polarizability for different parametrizations of the jellium model. The results are unexpectedly close to the experimental value of $78A^{3}$. While the global aspect of ground state electron density has no essential dependence on $\Delta$, the values of the density far from center or the cluster influence the value of the polarizability, fact which explains the spectrum of obtained values.

In the case of spherical homogeneous shell the obtained polarizability was $\alpha=92\mathring{A}^{3}$. The fact that the results from the gaussian parametrization of the jellium were closer to the experimental value raise the question of weather this aspect is a mathematical property of the equations involved, or simply the gaussian case is more realistic as discussed due to asymptotic tail of the core electrons from carbon atoms which must be taken into account in the geometry of jellium.

In order to test further the power of our equation and the validity of the approximations involved we calculate the polarizability of fullerenes with $180$ and $240$ carbon atoms also in a spherical symmetry and with a gaussian profile. The obtained results are around $260\mathring{A}^{3}$ for $C_{180}$ comparable with the RPA result of $300\mathring{A}^{3}$ zope2008static and $340\mathring{A}^{3}$ for $C_{240}$ compared with $432\mathring{A}^{3}$ from RPA zope2008static .

As the size of fullerene increases, the method starts to fail, one of the reasons being the fact that the spherical symmetry begins to be broken. Nonetheless, the results are still comparable with those from more involved methods gueorguiev2004quantum from computational point of view.