Energetics of charged metal clusters containing vacancies

TL;DR

This paper develops a theoretical approach combining Kohn-Sham calculations and analytical expansions to study the energetics of large charged metal clusters with vacancies, providing insights into their ionization, electron affinity, and vacancy formation energies.

Contribution

It introduces a new method that combines Kohn-Sham results with analytical expansions to analyze vacancy effects in metal clusters, including size dependence and asymptotic behaviors.

Findings

Derived expressions for ionization potential and electron affinity as corrections to work function.

Calculated vacancy formation energies and their size dependence for Na, Mg, Al clusters.

Analyzed the asymptotic behavior of vacancy formation mechanisms and their impact on cluster energetics.

Abstract

We study theoretically large metal clusters containing vacancies. We propose an approach, which combines the Kohn-Sham results for monovacancy in a bulk of metal and analytical expansions in small parameters $c_{v}$ (relative concentration of vacancies) and $R_{N,v}^{-1}$, $R_{N,v}$ being cluster radius. We obtain expressions of the ionization potential and electron affinity in the form of corrections to electron work function, which require only the characteristics of 3D defect-free metal. The Kohn-Sham method is used to calculate the electron profiles, ionization potential, electron affinity, electrical capacitance, dissociation, cohesion and monovacancy-formation energies of the small perfect clusters Na$_{N}$, Mg$_{N}$, Al$_{N}$ ($N \leq 270$) and the clusters containing a monovacancy ($N\geq 12$) in the stabilized-jellium model. The quantum-sized dependences for monovacancy-formation energies are calculated for the Schottky scenario and the "bubble blowing" scenario, and their asymptotic behavior is also determined. It is shown that the asymptotical behaviors of size dependences for these two mechanisms differ from each other and weakly depend on the number of atoms in the cluster. The contribution of monovacancy to energetics of charged clusters, the size dependences of their characteristics and asymptotics is discussed. It is shown that difference between the characteristics for the neutral and charged cluster is entirely determined by size dependences of ionization potential and electron affinity. Obtained analytical dependences may be useful for the analysis of the results of photoionization experiments and for the estimation of the size dependences of the vacancy concentration including the vicinity of the melting point.

Equations

$$R_N=N^{1/3}{r}_0,$$

$${\mathrm{IP}}_N=E_N^{+}(R_N)-E_N(R_N),$$

$${\mathrm{EA}}_N=E_N(R_N)-E_N^{-}(R_N),$$

$$\mu (R_N)={\mu }_0+\frac{{\mu }_1}{R_N}+\frac{{\mu }_2}{R_N^2}+O(R_N^{-3}),$$

$$\sigma (R_N)={\sigma }_0+\frac{{\sigma }_1}{R_N}+\frac{{\sigma }_2}{R_N^2}+O(R_N^{-3}),$$

$${\mu }_1=\frac{2{\sigma }_0}{\bar{n}},$$

$${\mu }_2={\mu }_1\left({\delta }_1-\frac{{\sigma }_0}{B_0}\right),$$

$$R_{N,v}=R_N{\left(1+c_v\right)}^{1/3},c_v=N_v/N.$$

$$R_v={\left(\frac{3}{4\pi n_a{c}_v}\right)}^{1/3}\gg r_0={\left(\frac{3}{4\pi n_a}\right)}^{1/3},$$

$${\mathrm{IP}}_{N,v}=W_{\mathrm{eff},v}-\frac{{\mu }_1}{R_{N,v}}-\frac{{\mu }_2}{R_{N,v}^2}+\frac{e^2}{2{\mathcal{C}}_{N,v}},$$

$${\mathrm{EA}}_{N,v}=W_{\mathrm{eff},v}-\frac{{\mu }_1}{R_{N,v}}-\frac{{\mu }_2}{R_{N,v}^2}-\frac{e^2}{2{\mathcal{C}}_{N,v}},$$

$${\mathcal{C}}_{N,v}=R_{N,v}.$$

$${\mathrm{IP}}_N=W_0+\frac{\alpha e^2}{N^{1/3}},\alpha =\frac{1}{2r_0}-\frac{{\tilde{\mu }}_1}{e^2},$$

$${\mathrm{EA}}_N=W_0-\frac{\beta e^2}{N^{1/3}},\beta =\frac{1}{2r_0}+\frac{{\tilde{\mu }}_1}{e^2},$$

$$-\frac{\hbar^{2}}{2m}\nabla^{2}\Psi({\bf r})+\Big{[}v_{\rm eff}(r)+\sum\limits_{i=1}^{N_{v}}\delta v_{{\rm eff},v}(\textbf{r}-\textbf{R}_{i})\Big{]}\Psi(\textbf{r})=\varepsilon_{b}\Psi({\bf r}),$$

$$\Psi (\mathbf{r})=\psi (\mathbf{r})u(\mathbf{r}),$$

$${\epsilon }_b={\epsilon }^{(0)}+{\epsilon }^{(1)},$$

$${\epsilon }^{(0)}=T_0+⟨\delta v_{\mathrm{eff},v}{⟩}_{R_v},$$

$$\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+\delta v_{\mathrm{eff},v}(r)-T_0\right]u_{\mathrm{WS}}(r)=0$$

$$\nabla u_{\mathrm{WS}}(r){|}_{r=R_v}=0,u_{\mathrm{WS}}(r){|}_{r=L_v}=0,$$

$$u_{\mathrm{WS}}(r)=\frac{A}{\sqrt{4\pi }}\frac{\sin [q_0(r-L_v)]}{q_{0}r}$$

$$\tan [q_0(R_v-L_v)]-q_0{R}_v=0,$$

$$T_0=\frac{{\hslash }^2{q}_0^2}{2m}.$$

$$T_0=\frac{3{\hslash }^2{L}_v}{2m{r}_0^3}{c}_v.$$

$$⟨\delta v_{\mathrm{eff},v}{⟩}_{R_v}=\frac{3c_v}{4\pi r_0^3}{\int }_0^{R_v}dr4\pi r^2\delta v_{\mathrm{eff},\mathrm{v}}(r),$$

$${\epsilon }^{(0)}=A_1{c}_v+O(c_v^2),$$

$$\left[-\frac{{\hslash }^2}{2m}{\nabla }^2-\frac{{\hslash }^2}{m}\sum _{i=1}^{N_v}\frac{\nabla u_{\mathrm{WS}}(\varrho )}{u_{\mathrm{WS}}(\varrho )}\nabla +v_{\mathrm{eff}}(r)-{\epsilon }^{(1)}\right]\psi (\mathbf{r})=0,$$

$${\epsilon }^{(1)}=\frac{{\hslash }^2{k}^2}{2m_{\mathrm{eff}}}.$$

$$W_{\mathrm{eff},v}=W+\delta W_v^{\mathrm{bulk}},$$

$$\delta W_v^{\mathrm{bulk}}=-{\epsilon }^{(0)}.$$

$$\mu (x)=e\varphi (x)+{⟨\delta v⟩}_{\mathrm{WS}}\theta (-x)+\frac{\delta G}{\delta n(x)}=\text{mbox}const,$$

$$\frac{d^2\varphi }{d{x}^2}=-4\pi e[n(x)-\rho (x)],\rho (x)=\bar{n}\theta (-x).$$

$$W=-\bar{\mu }.$$

Full text

The Kohn-Sham method is used to calculate the electron profiles, ionization potential, electron affinity, electrical capacitance; dissociation, cohesion and monovacancy-formation energies of the small perfect clusters NaN, MgN, AlN ($N\leq 270$) and the clusters containing a monovacancy ($N\geq 12$) in the stabilized-jellium model. The quantum-sized dependences for monovacancy-formation energies are calculated for the Schottky scenario and the ‘‘bubble blowing’’ scenario, and their asymptotic behavior is also determined. It is shown that the asymptotical behaviors of size dependences for these two mechanisms differ from each other and weakly depend on the number of atoms in the cluster. The contribution of monovacancy to energetics of charged clusters, the size dependences of their characteristics and asymptotics is discussed. It is shown that difference between the characteristics for the neutral and charged cluster is entirely determined by size dependences of ionization potential and electron affinity. Obtained analytical dependences may be useful for the analysis of the results of photoionization experiments and for the estimation of the size dependences of the vacancy concentration including the vicinity of the melting point.

I Introduction

Frenkel theory of melting assumes an abrupt increase of the vacancies concentration at the triple point, as well as the decrease in monovacancy-formation energies with increasing their concentration 179 ; Ubbelohde ; Berry-2013 . The equilibrium vacancy concentration is estimated from thermodynamic considerations based on the monovacancy-formation energy, the magnitude of which, in turn, can be extracted from positron annihilation spectroscopy BPR-2015 . At the melting point, the relative concentration of vacancies in metals is a fraction of a percent. Despite such small concentrations, vacancies significantly influence properties of solids. In the case of radiation damage in 3D metals or metal islands, the vacancy concentration can be even tens of percent.

Small metal nanoclusters at low temperatures can be in superconducting state, which results in a strong modification of the energy spectrum. It is known that superconducting correlations depend crucially on the density of states near the Fermi energy. Certain shapes of nanoclusters support highly enhanced density of states near electronic shell closings, see, e.g., recent works KresinOvchin-2006 ; HalderKresin-2015 . For ideal nanoclusters, the highest energy occupied electronic levels become strongly degenerate at spherical shell closings (‘‘magic’’ numbers). The presence of vacancies, as shown below, will lead to a change in the magic numbers of atoms and might result in smearing of the effect of shell structure on superconducting correlations including gap in excitation spectrum.

Initially, in the experiments it was established that the temperature melting of clusters on a substrate and free clusters decreases with a decrease in their size Sambles ; Borel . Interpretations of this mesoscopic phenomena were developed in a number of papers Yang-2007 ; Hendy ; Guisbiers ; Safaei ; Luo ; Chandra ; qi-2016 . A popular point of view is based on thermodynamical considerations: near the temperature melting – the smaller cluster size, the lower a monovacancy-formation energy  while the vacancy concentration is not size- dependent Yang-2007 ; Guisbiers .

Modern mass spectroscopic and calorimetric methods, allowing to study in detail the process of premelting and postmelting in metal clusters consisting of a countable number of atoms Jarrold-2005 ; Haberland-2009 ; Jarrold-2010 ; Zamith , have shown that the melting temperature is characterized by an oscillatory size dependence, and also has dimensional anomalies (for example, for Al), poorly described by simple models Jarrold-2005 ; Jarrold-2010 . Also, in the melting process the diffusion of surface vacancies into bulk is more favorable for clusters with unfilled electronic shells than for clusters with the magic number of atoms Haberland-2009 . These facts stimulate growing interest in the description of phase transition from the solid to the liquid state, as the configuration excitation of voids-vacancies in the clusters. For example, the question about the size of the monovacancy-formation energies, the vacancy concentration and a relation to the melting process remains open.

For the first time, mass-spectrometric measurements of dissociation energy cluster ions Na${}_{N}^{+}$ were reported in Brechignac , and for Al${}_{N}^{+}$ in Martin . Traditionally, according to such data, and also based on the measured ionization potentials, the cohesion energy of neutral clusters is calculated.

Recently, the ionization potential of Al32-95 clusters and its temperature dependence in the range 65 – 230 K were measured in Kresin-2015 . With increasing temperature, a decrease of the ionization potential is not significant ($\sim 10$ meV). The melting point of such clusters is in the range 600 – 700 K Jarrold-2005 ; Jarrold-2010 .

The energy characteristics of solid clusters have been intensively studied within different models and approaches, including triaxially deformed ordinary jellium and ab initio simulations (see 38 ; 7 ; Landman-95 ; Vieira ; Akola-2000 ; Kanhere ; Aguado-2009 and references therein). We use model of spherical stabilized jellium 25 that does not contain adjustable parameters. It is also convenient and transparent for the analysis of the role played by vacancies in cluster melting or other excitation processes.

However, with the exception of the works Alonso-88 ; Itoh , self-consistent calculations for monovacancy-formation energy in clusters and the impact on it of electronic spectrum quantization has not yet been addressed. Therefore, one of the actual problems that can be formulated in connection with melting of small-sized aggregates, is the study of the size of its electron affinity and ionization potential in the case of clusters contained vacancies. In our model, the monovacancy is represented in the form of spherical void of atomic size in a homogeneous positively charged background due to the ions Alonso-88 ; Ziesche .

An electron binding energy for small perfect clusters may be calculated numerically only. In the opposite case of large clusters the binding energy is close to the binding energy for the 3D metal with the first size correction due to surface curvature. In Landman-95 ; Seidl-96 ; Seidl-1998 an exhaustive survey of the semiclassical description of IP and EA in the form of the first dimensional correction to the work function is given. In our work the second size corrections are also taken into account and the developed procedure is suitable for calculations of the ionization potential and electron affinity of metal clusters containing point defects and impurities. As an example, vacancies are considered. This approach can be easily adapted for calculation of the positron attachment energy for such clusters.

We propose a method, which combines several approaches. The first one is based on a density-functional solution in the stabilized jellium model for metal monovacancy in the bulk ignoring an external surface. This problem was addressed by us earlier. We found a shift in the energy of the ground state of the electrons BVP-2014-1 and positrons BVP-2014-2 due to the presence of a ‘‘subsystem’’ of vacancies in 3D metal.

In the first approximation it is assumed that vacancies are noninteracting because of their small concentration. However, there are experiments showing that in some metallic systems the partial ordering of vacancies is observed Amelinckx , i.e. the energy of vacancies formation has to be concentration dependent.

The second one uses a solution for a defect-free metal in presence of an external flat surface, but with lower atomic density. The lower density of the atoms is due to the existence of superlattice of vacancies of relative concentration $c_{v}$ in the defected metal. Using $c_{v}$ as a small parameter, all metal characteristics are evaluated as series expansions. In these expansions the zero-order terms are the characteristics of defect-free metal, and the first-order terms are expressed through them. This approach allows to obtain the vacancy dependence of electron work function for 3D metal. These results are then used to evaluate an asymptotic size dependence of the ionization potential of a spherical metal cluster and the electron affinity.

The consistent procedure for the calculation of a size dependence of ionization potential and electron affinity of a large spherical metal cluster containing vacancies is presented. In the framework of effective medium approach, the perturbation theory over the small parameters $c_{v}$, $R_{v}/R_{N,v}$ and $L_{v}/R_{v}$ is proposed for the ionization potential and electron affinity ($R_{v}$ is the average distance between vacancies and $L_{v}$ is the electron-vacancy scattering length).

For small clusters Rb, K, Na, Li, Mg and Al containing monovacancy we performed the Kohn-Sham calculations of the electron and effective potential profiles; the ionization potential and an electron affinity; the dissociation, the cohesion, the vacancy formation energies. For demonstration and analysis of results, we limited ourselves to metals Na, Mg and Al.

The size dependences for monovacancy-formation energies are calculated for the Schottky scenario and the ‘‘bubble blowing’’ scenario, and their asymptotic behavior is also determined. It is shown that the asymptotical behaviors of size dependences for these two mechanisms differ from each other and weakly depend on the number of atoms in the cluster. The contribution of monovacancy to energetics of charged clusters, the size dependences of their characteristics and asymptotics is discussed. It is shown that difference between the characteristics for the neutral and charged cluster is entirely determined by size dependences of ionization potential and electron affinity.

II Ionization Potential and Electron Affinity

II.1 General relation

By definition, the ionization potential and the electron affinity of a metal cluster with the $N$ atoms and radius

[TABLE]

where $r_{0}$ is the radius of Wigner-Seitz cell per unit atom, have the form

[TABLE]

where $E_{N}^{+}\equiv E_{N}^{N_{e}-1}$, $E_{N}^{-}\equiv E_{N}^{N_{e}+1}$ and $E_{N}$ are total energies of charged and neutral spheres; $N_{e}=ZN$ is the total number of electrons in a neutral cluster, $Z$ is a valency of metal.

The basis of the semiclassical approximation is the expansion of the electron chemical potential $\mu$ of valence electrons and the surface energy per unit area $\sigma$ of the neutral cluster in powers of the inverse radius $R_{N}^{-1}$ (liquid drop model) 1985 ; pppppp :

[TABLE]

where $\mu_{0}$ and $\sigma_{0}$ correspond to flat surface ($R_{N}\rightarrow\infty$).

From the condition of mechanical equilibrium of the cluster, the sum rules were obtained in Refs. 177 ; iakPog-1995 ; P-1995 ; keiPog-1996 . In particular

[TABLE]

where $\delta_{1}=\sigma_{1}/\sigma_{0}$. Next, it is convenient to use the values $\widetilde{\mu}_{1}\equiv\mu_{1}/r_{0}$, $\widetilde{\mu}_{2}\equiv\mu_{2}/r_{0}^{2}$ having a dimension of energy.

Consider a metallic cluster consisting of $N$ atoms and containing $N_{v}$ vacancies. Then radius of the cluster is

[TABLE]

Clusters of atoms have a structural periodicity, which is not translational, but rather has a property of ‘‘spherical periodicity’’ Martin-1996 . This periodicity is due to spherical layers of atoms (atomic shells or coordination spheres). By this principle, one can also consider the vacancy ‘‘subsystem’’ in a cluster with distribution over spherical layers. We divide conditionally the cluster into $N_{v}$ Wigner-Seitz supercells of radius

[TABLE]

where $n_{a}$ is the concentration of atoms. Respectively the supercellular ‘‘muffin-tin’’ potential is replaced by a spherical symmetric one.

For a large cluster with dilute subsystem of vacancies, we will use the expressions

[TABLE]

where $W_{{\rm eff},v}$ is the effective electron work function of 3D metal containing vacancies, $-e$ is the electron charge. The last term (7) is the charging energy of metallic sphere with electric capacitance

[TABLE]

The analysis of experimental data for small metallic clusters is usually carried out according to the formulae

[TABLE]

which will be used by us as reference expressions.

II.2 Volume term for vacancy shift of the grand energy of electrons

The translational symmetry of the lattice of a solid is at the basis of calculation of the electron work functions. Therefore, to take into account the contribution of vacancies to work function, we have to assume their periodic location in the form of a ‘‘superlattice’’ in a metal. In this case the vacancy shift of the bottom of the electron conductivity band can be introduced.

The problem of description of our system within simple models is that the bulk of the metal is not homogeneous due to vacancies, so the density functional method in the jellium model must be solved as a three-dimensional problem. Within the one-dimensional problem, it is impossible to describe a set of spherically symmetric vacancies and a planar external metal surface.

The wave function of the ground state of the electron is determined by the Schrödinger equation

[TABLE]

in which the vacancies are centered at the points $\textbf{R}_{i}$ ($R_{i}<R_{N,v}$). One-electron effective potentials are represented in such manner that the spherically symmetric potential $v_{\rm eff}(r)$ forms the bottom of the conduction band and the surface barrier, and $\delta v_{{\rm ef},v}(\textbf{r}-\textbf{R}_{i})$ forms the ith vacancy (Fig. 1). In (10) the energy $\varepsilon_{b}$ is a shift in the bottom of the conduction band due to the presence of vacancies.

Suppose that in 3D metal the potential field of the supercell has translational periodicity. Then the electron wave function of the ground state can be written in the form,

[TABLE]

where the function $\psi({\bf r})$ (for 3D metal it is $e^{{\rm i}\textbf{k}\textbf{r}}$) is modulated on the vacancy scale by the function $u({\bf r})$. The function $u({\bf r})$ is the superposition of the usual Wigner-Seitz solutions $u_{\rm WS}(|\textbf{r}-\textbf{R}_{i}|)$ inside cells, centered at $\textbf{R}_{i}$.

Next, $\varepsilon_{b}$ can be represented as the sum

[TABLE]

[TABLE]

where $T_{0}$ is the cellular energy eigenvalue; $\langle\delta v_{{\rm eff},v}\rangle_{R_{v}}$ is averaged over the volume of the supercell the contribution of potential energy from the electron-vacancy potential $\delta v_{{\rm eff},v}(r)=v_{{\rm eff},v}(r)-\bar{v}_{\rm eff}$, where $\bar{v}_{\rm eff}$ is the position of the bottom of the conduction band in absence of vacancies (Fig. 1).

The value $T_{0}$ is determined from equation for supercell

[TABLE]

with the boundary conditions

[TABLE]

where $L_{v}$ is the scattering length of electron on the vacancy (Table I). For an electron, the vacancy represents a potential hillock, therefore $L_{v}>0$.

Table I. The results of the Kohn-Sham calculation of $L_{v}$ for the electron scattering on monovacancy BVP-2014-1 ; Bab .

[TABLE]

Substitution of the expression for the electron wave function

[TABLE]

in the second boundary condition gives the equation

[TABLE]

from which we obtain

[TABLE]

This approach was first used in Spr .

The interaction of an electron with vacancies can also be described using the Fermi optical approximation

[TABLE]

Here it should be noted that the smaller $c_{v}$ the better agreement between results of calculations using (16) (18) and (19). Solving Eqs. (17) and (19), for example, for Al with $c_{v}=0.01$ leads to the values $T_{0}=0.0388$ eV and 0.0296 eV, respectively, which evidences that an accuracy of the optical approximation is lower.

The value of $T_{0}$ is determined by the $s$-phase of scattering ($\rightarrow-L_{v}k$ as $k\rightarrow 0$) for $k=q_{0}$ and can be tested by the value of residual resistivity of vacancies. All phases of electron scattering for a monovacancy give contributions to the resistivity, but the main contribution is produced by the $s$-phase of scattering for $k=k_{\rm F}$ (see Table 1 in BVP-2014-1 ). The residual resistivity of vacancies estimated in BVP-2014-1 turned out to be 2-3 times lower than the experimental values.

The value $\langle\delta v_{{\rm eff},v}\rangle_{R_{v}}$ in (13) is determined in the mean-field approximation by the expression

[TABLE]

in which it is convenient to substitute $R_{v}=r_{0}c_{v}^{-1/3}$.

Using the numerical solution of the problem of electron scattering on the vacancy potential BVP-2014-1 we get

[TABLE]

where $A_{1}=4.10$ eV and 13.3 eV for Na and Al, respectively.

We substitute (11) and (14) into the equation (10), which after simple transformations can be rewritten in the form

[TABLE]

where $\varrho\equiv|\textbf{r}-\textbf{R}_{i}|$.

The Eq. (22) contains the component of effective potential in the form of a cross-term. It can be treated as a perturbation. Earlier, analogous procedures were carried out in 3D perfect metals for calculations of ground state energy $\varepsilon^{(0)}$ and effective masses $m_{\rm eff}$ of electrons by Bardeen 244 , Cohen and Ham 245 and for positrons by Stott and Kubica 246 ,

[TABLE]

In Refs. 244 and 245 the set of cellular wave functions was used. In Ref. 246 the function $\psi({\bf r})$ was expanded in the set of plane waves in a crystal.

It turned out that the influence of vacancies on effective masses of electrons BVP-2014-1 and positrons BVP-2014-2 in melting point ($c_{v}\approx 10^{-3}$) is insignificant. In order to estimate the effective mass, we used in BVP-2014-1 the Bardeen approach 244 with $s$- and $p$-phase of scattering for $k=q_{0}$. For $c_{v}=10^{-2}$ the calculated values of the electron effective mass $m_{\rm eff}=1.006$, 1.01, 1.012 and 1.014 for Na, Cu, Al and Pb, respectively. The small excess of $m_{\rm eff}$ over $m$ indicates that the repulsion prevails over the attraction in the scattering of electrons by the supercell potential with the radius $R_{v}$.

II.3 $c_{v}$-expansions

We propose a method which combines two approaches. The first approach is based on a self-consisten solution for the monovacancy in 3D metal neglecting its surface. The second approach is based on a stabilized jellium model for uniform 3D metal with flat boundary, but the atom density of a metal is decreased due to the presence of vacancies superlattice of relative concentration $c_{v}$. In this case, the effective electron work function $W_{{\rm eff},v}$ can be presented as a sum

[TABLE]

where $W$ is traditionally calculated by the density functional method a characteristic consisting of a volume component and a surface dipole barrier, and

[TABLE]

For a semi-infinite metal [the $x$ axis is perpendicular to the interface metal ($x\leq 0$) - vacuum ($x>0$)] containing vacancies the equilibrium profile $n(x)$ satisfies the the Euler-Lagrange equation

[TABLE]

where the electrostatic potential $\phi(x)$ is determined by the integration of the Poisson equation

[TABLE]

Here $\mu$ is the electron chemical potential; $\left\langle\delta v\right\rangle_{\rm WS}$ is the stabilization potential, $\theta(-x)$ is the unit step function, $G[n]$ is the universal functional corresponding to ordinary jellium model with energy per electron $\varepsilon_{\rm J}=\varepsilon_{\rm t}+\varepsilon_{\rm xc}$ (for detail see Ref. 25 ); distribution of the positive charge $\rho({\bf r})$ is homogeneous inside the metal and zero outside it.

The condition $\mu(x)$ = const follows from the equivalence of the choice of the coordinate in (25). It is convenient to take $x=-\infty$, where the gradient terms (g. t.) vanish. Then the work function of the electrons can be found as

[TABLE]

Let us make a comparison between a defect-free half-infinity metal and a metal containing vacancies. The scenario for the formation of vacancies is not important in this context. It may be the Schottky mechanism or the ‘‘bubble blowing’’ mechanism Pog-94 . It is important that the number of atoms in the sample is so large, that in the spaces between vacancies a density of atoms is the same as in perfect sample.

We average the charge density over the vacancy supercell. As a result, the electroneutrality condition in the bulk of a ‘‘fictitious, defect-free’’ metal is

[TABLE]

where $Z$ is the valency of metal. The electron density satisfies the condition $\frac{4}{3}\pi r_{s}^{3}\bar{n}=1$; $r_{s}=r_{0}/Z^{1/3}$.

Carrying out the series expansion with respect to the small dimensionless parameter $c_{v}$ in condition (28) and limiting ourselves to linear terms, we have

[TABLE]

where $\bar{n}^{0}$ is the density of homogeneous electron gas for $c_{v}=0$.

By analogy with the Refs. iakPog-1995 ; P-1995 ; keiPog-1996 , in which a method is developed for determination of size corrections of energy characteristics of spherical clusters in the stabilized jellium ($R^{-1}$ is a small dimension parameter), we represent characteristics of a metal containing vacancies as

[TABLE]

This allows us to expand the equations (25) and (26) into series. Restricting ourselves to zero and first-order terms of expansion (superscripts ‘‘0’’ and ‘‘1’’, respectively), we have

[TABLE]

[TABLE]

[TABLE]

In (32) $g^{0}\equiv n^{0}(x)\varepsilon_{\rm J}\left(n^{0}(x)\right)$ is the energy density in LDA.

From condition $x=-\infty$, we have sum-rules

[TABLE]

[TABLE]

where the energy per electron in the stabilized jellium $\varepsilon_{\rm SJ}=\varepsilon_{\rm J}+\tilde{\varepsilon}$, $\tilde{\varepsilon}$ is the electrostatic self-energy of the uniform negative in the Wigner-Seitz cell 25 .

Note that $\bar{\mu}_{0}$ in Eq. (3) and $\bar{\mu}^{0}$ in Eq. (30) are the same, and $\bar{\mu}_{1}$ and $\bar{\mu}^{1}$ have different dimensions and are expressed in terms of characteristics of a defect-free metal with a flat surface.

Integration in Eq. (33) gives

[TABLE]

Due to the fact that the vacancy shift of the electrostatic potential in the bulk $\bar{\phi}^{1}$ is only due to a decrease in average electron and ions densities in the metal, by analogy with self-compressed clusters P-1995 ; keiPog-1996 , where the effect is inversed, we have instead of (36)

[TABLE]

Here the primes denote derivatives with respect to $\bar{n}^{0}$.

Using expression for the bulk modulus, $B^{0}=(\bar{n}^{0})^{3}(\bar{\varepsilon}_{\rm SJ}^{0})^{\prime\prime}$, the formulas (29) and (37), we finally have for first-order term

[TABLE]

As a result of this approach

[TABLE]

and vacancy contribution

[TABLE]

as a whole is expressed only through the characteristics of a defect-free metal. For defect-free Al and Na $W^{0}=4.30$ eV and 2.93 eV, respectively.

In the above exact formulas, the only $\bar{\phi}^{0}$ and $(\bar{\phi}^{0})^{\prime}$ terms require self-consistent calculations. Table II shows the values of the components formula (38). As we see, the contribution from the surface barrier $\sim(\bar{\phi}^{0})^{\prime}$ is very significant, competing with the bulk contribution in (38) (a sum of terms in parentheses). In general, $W^{1}$ are negative in sign. In almost all cases, the major contribution of the effect is in the magnitude of $\delta W_{v}^{\rm bulk}$ in (40).

For Al at the melting point $c_{v}\approx 10^{-3}$. Consequently, the contribution of equilibrium bulk vacancies to the electron work function is approximately $-$0.2 eV. According to the paper Gul , the concentration $c_{v}$ should be much higher, then, respectively, the effect of vacancies will increase by many times.

Table II. The results (in eV) of the Kohn-Sham calculations for componets of formula (38) and vacancy contribution $W_{\rm eff}^{1}$ (40) to the electron work function.

[TABLE]

Small independent parameters $R^{-1}$ and $c_{v}$ appear, when applying these approaches to clusters.

II.4 Size vacancy contribution

Suppose that in large spherical clusters the potential field vacancy suppercells has a ‘‘spherical periodicity’’ and function $\psi({\bf r})$ in (11), which varies on the scale of the whole cluster, is modulated by the function $u({\bf r})$. Then the boundary conditions for $\psi({\bf r})$ on boundary lead to a discrete energy spectrum of electrons.

In the simplest case, assuming infinitely deep square well, the energy of the ground state of electrons in the potential well $v_{\rm eff}(r)$ of the whole cluster was found in PR-2017 . This approach uses the boundary condition

[TABLE]

for Eq. (22) and allows to derive an analytical expression for vacancy quantum correction. The solution was obtained by expanding the function $\psi({\bf r})$ in terms of the complete orthonormal set of eigenfunctions corresponding to the deep potential well. In this case, the term $\varepsilon^{(1)}$ in Eq. (12) can be written as

[TABLE]

The operator $\delta V(r)$ in (22) is used as a perturbation,

[TABLE]

The diagonal matrix element of the perturbation is the potential field $\delta V(r)$,

[TABLE]

averaged over the ground state with the quantum numbers $n_{r}=1$, $l=0$

[TABLE]

As a result of the integration in Eq. (44) (see Appendix A), using (4), (7), (42) and (80), the final expression of the ionization potential becomes

[TABLE]

where $D_{1}\approx 2.7$. Here we neglected the size dependence in $\delta W_{v}^{\rm bulk}$ and the vacancy dependence in $\widetilde{\mu}_{1}$ and $\widetilde{\mu}_{2}$. In our work PR-2017 we used $\delta W_{v}^{\rm bulk}$ instead of $W_{{\rm eff},v}$.

The corresponding expression for electron affinity ${\rm EA}_{N,v}$ is obtained from (46) by replacing $\alpha\rightarrow-\beta$.

In deriving $\varepsilon^{(1)}$, we used the boundary condition (41). It does not take into account electron spillover. The simplest way to account for the spillover is to use perturbation theory Kubo-66 , taking into account the finite depth of the spherical well $\approx\bar{v}_{\rm eff}$, $\bar{v}_{\rm eff}=-6$ and $-16$ eV for Na and Al, respectively (Fig.1). There is no analytical solution for a spherical well, and for a cluster-cube the result is known PogKurbVas-2005 . This means that that result $\sim\hbar^{2}$ in Eq. (46) for hard wall-cube must be multiplied by

[TABLE]

where $k_{0}=\sqrt{-2m\bar{v}_{\rm eff}}$. The account of the spillover will result in the appearance the term $-A/N^{1/3}$ in the last parenthesis in the (46), where $A\approx 1$ for Na and Al. Thus it is obvious that in the derivation of the term we neglected the terms proportional $1/N$.

Therefore, taking into account that $c_{v}\ll 1$, condition (41) limits the value of $c_{v}$ from below. As a result of the electron spillover neglection, the applicability of formula (46) is possible under the condition

[TABLE]

The lower size limit where the concepts of supercells work is two supercells shells ($N_{v}\geq 55$). The use of this technique for one supercells shells ($N_{v}=13$) is already an extrapolation.

The straightforward application of inequality (6) and the boundary condition (41) in the paper PR-2017 leads to the erroneous application of formula for ${\rm IP}_{N,v}$ to the case of a monovacancy in a cluster (see Fig. 3 in PR-2017 ). The perturbation theory made it possible to introduce more realistic description of the cluster. Because of this, without changing the form of (80), the limits of applicability (47) of the theory for a dilute subsystem of vacancies in a cluster are formulated in the case of a discrete electron spectrum.

For $c_{v}=0.01$ and the first two coordination spheres of supercells, the total number of atoms in the cluster is $N=1300$ and 5500, which corresponds to $R_{N,v}\approx 1.8$ and 2.9 nm for Al, respectively.

However, the choice of the applicability criteria of (46) is more consistent from the condition of expansion of the electron chemical potential in power series of the inverse radius ($c_{v}=0$)

[TABLE]

In its turn, this leads to the values $N\geq 9.68\cdot 10^{3}$ $(R\approx 4.5$ nm) and $N\geq 5.45\cdot 10^{4}$ ($R\approx 6$ nm) for Na and Al, respectively.

The large clusters are detectable in the experiment. For example, Hoffmann reported on a photoelectron spectroscopy and determination of the steps of the Coulomb staircase $(E_{C}=e^{2}/R_{N})$ of ionized clusters Al${}_{N\leq 32000}^{-}$ with an accuracy of approximately $10^{-2}$ eV.

The value of $\delta_{1}$ in (4) was determined repeatedly Pog-94 . We use calculated within a stabilized jellium model values $\delta_{1}=\sigma_{1}/\sigma_{0}=0.32r_{0},\,0.57r_{0}$ for Na and Al, respectively [$\sigma_{1}=\gamma/2$ in Table VI of work Ziesche ]. The experimental values of $W_{0}=2.75$ eV, IP${}_{1}=5.14$ eV, $\sigma_{0}/B_{0}=0.70$ $a_{0}$ for Na and $W_{0}=4.28$ eV, IP${}_{1}=5.99$ eV, $\sigma_{0}/B_{0}=0.40$ $a_{0}$ for Al are taken from 86 .

In Fig. 2 the dependences (9) and (46) are shown for Na and Al. Due to the fact that the formula (9) gives the value IP1 only 10% higher than the experimental one, and the experimental values of $W_{0}$ for a number of materials are strongly dependent on measurement techniques, it was assumed from measurements IPN≤100 for magic numbers (see Fig. 28 in 38 ) using (9) to determine the true values $W_{0}$. On the other hand, the the last term in (46) that is conditioned by quantization lead to non-linear behavior of IP$(N^{-1/3})$ and EA$(N^{-1/3})$ in the region of small $N$.

The presented analytical approach seems to be promising for experimental estimation of the concentrations of point defects or impurities in metal clusters. To this end, it is first needed to calculate the scattering length of electrons on the corresponding defect in 3D metal. In particular, the problem of vacancies concentration in the cluster at the melting point can be solved. Using Fig. 2 as grid of values, if the experimental value of IPN,v resides at one of the curves, the value $c_{v}(N)$ is fixed for a given temperature.

II.5 Small clusters containing monovacancy

For simplicity, the density of the positive charge background of defected and perfect (free-defect) clusters is chosen in the form

[TABLE]

and

[TABLE]

respectively.

Wave functions and eigenvalues of energies $\varepsilon_{j,v}$ are the solution of the Kohn-Sham equations with effective single-electron potential, including electrostatic and exchange-correlation potential in the LDA. The energy is measured from the vacuum level, that is, on the energy of an electron with zero kinetic energy, located far from the sample ($r\gg R_{N,v}$).

The electrostatic potential is a solution of the Poisson equation for a fixed condition

[TABLE]

where $Q$ is the total charge of the cluster.

Fig. 3 shows the profiles of the electron distribution and effective potential for a defect-free cluster and cluster with vacancy, containing the same number of atoms $N=12$. The vacancy radii for Na and Al clusters are shown in the figure. Their radii differ according to the definition (5). The insets show the tails of electron profiles and the potentials far from clusters. Despite the fact that electronic distributions are rapidly decreasing functions of radial distance $r$, the potential tails extend far (the calculation was carried out to $r\approx R_{N}+900\,a_{0}$). For charged clusters, the electrostatic potential decreases asymptotically beyond $\sim Q/r$. It is interesting to compare the profiles of electrons and potentials in Fig. 3 for Na12 and Al12 with analogous profiles for 3D metals in Fig. 1. For large clusters, the spatial profile becomes similar to the profile near the surface of a semi-infinite metal, containing a large number of the Friedel oscillations. Obtained profiles allow us to calculate the total energy of the cluster $E^{\pm}_{N,v}$ (and $E^{\pm}_{N}$).

Fig. 4 shows the results of direct calculations of IP and EA by the formula (2), as well as the asymptotics (9). A difference between perfect clusters and defective clusters can be understood from this figure as well as from Tables III and IV. For Na with the increase of $N$, starting from 12, this difference can be 0.1 – 0.5 eV (for Mg and Al about twice as much). Maximum difference is observed at the transition from a completely filled shell to an empty one. As $N$ increases, this difference is leveled. For clarity, the results of the calculations are given in coordinates $N^{-1/3}$. For clusters with a monovacancy $c_{v}=1/N$, therefore there is a correspondence $N^{-1/3}=c_{v}^{1/3}$. In our case $c_{v}\rightarrow 0$ for $N\rightarrow\infty$. If the vacancy is not single, but their concentration is low (vacancies do not interact with each other), according to our results, it is possible to track the dependence of energy characteristics on the concentration of vacancies.

The ionization potential and the electron affinity show a strong oscillatory behavior due to the spherical shell structure. They tend to $W_{0}$ asymptotically slow enough, due to high orbital degeneracy and large angular quantum numbers $l$. According to the results of experiments, the oscillations should be much weaker. Going beyond LDA and using local spin density approximation (LSDA) allows to reduce oscillations Ekardt-84 .

Using the Koopmans’ theorem, the formulas (2) can be rewritten Perdew-88 ; Sabin in form

[TABLE]

where $\varepsilon_{N,v}^{\rm HO}$ / $\varepsilon_{N,v}^{\rm LU}$ and $\mathcal{C}_{N,v}^{\pm}$ are the energy of the upper occupied / lower unoccupied orbital and electrical capacitances, respectively.

The spectra for perfect Na and Al clusters, shown in Fig. 5 demonstrate the filling of electronic shells as the number of electrons is increased. For partially filled shells $\varepsilon_{N,v}^{\rm HO}=\varepsilon_{N,v}^{\rm LU}\approx\mu(R_{N,v})$. The maximum values of $\varepsilon_{N,v}^{\rm HO}$ correspond to completely filled shells, and magic numbers of atoms $N^{*}$ for spherical perfect clusters and clusters with vacancy does not coincide in all cases. For Na the obtained values are $N^{*}=$2, 8, 18, 20, 34, 40, 58, 68, 90, 92, 106, 132, 138, 168, 186, 196, 198, (230), 232, (252), 254; for Mg $N^{*}=1$, 4, 9, 10, 17, 20, 29, 34, 45, 46, 53, 66, 69, 78, 93, 98, (99), {115}, (116), {126}, 127, 134, 153, 156, 169, 178, 199, 204, 219; for Al $N^{*}=$6, (30), 44, 46, 52, 62, {66}, (84), (102), {104}, {136}, 146, {154}, (180), (202), 204. In round brackets, the obtained values are given for defective clusters which do not coincide with the corresponding values for perfect clusters, and in braces – on the contrary.

As $R_{N,v}$ increases, the values $\varepsilon_{N,v}^{\rm HO}$ and $\varepsilon_{N,v}^{\rm LU}$, oscillate and tend to $\mu(R_{N,v})$ for $R\rightarrow\infty$. Amplitude of oscillations decreases approximately as $R_{N,v}^{-3}$.

Let us denote the difference

[TABLE]

and return to Fig. 4. At first sight, the sign $\Delta({\rm IP}_{N})>0$ is unexpected (circles are placed above the points for the same $N$). Exceptions are clusters with such $N$ that maximum contribution is given by levels with low $l$ ($s,p$- and partially $d$-orbital). In Fig. 4 these narrow areas are enclosed between vertical dashed lines.

Table III. Computed quantities ${\rm IP}$ and ${\rm EA}$ (in eV) for Na clusters. Experimental values of ${\rm IP}_{\rm exp}$ (Na17-94 – 38 , Na132-140 – Kostko ) and ${\rm EA}_{\rm exp}$ (Na17-196 – Kostko ).

[TABLE]

From an analysis of the asymptotic behavior of ${\rm IP}_{N}$ and ${\rm IP}_{N,v}$ the basic vacancy dependence is contained in the work function $W_{{\rm eff},v}<W_{0}$ (see (39) and Table II).

In the case of small defected clusters, the perturbation of the vacancy becomes essential, provided concentration $c_{v}\sim R_{N}^{-3}$. It follows from Fig. 3 that the behavior of $v_{{\rm eff},v}(r)$ is such that the electrons are squeezed out by a vacancy from the center of the cluster to its surface and are grouped, mainly, in the spherical layer $r_{0}<r<R_{N}$. And when integrating in spherical coordinates, the expression for the total energy gives the main contribution to the energy. This is confirmed by the spectral values of the energies, corresponding to the points (circles) in Fig. 5. As an example, we give the values $\varepsilon_{n_{r},l}$ ($n_{r}$ and $l$ are the radial and orbital quantum numbers) for perfect and defective (in parentheses) clusters: $\varepsilon_{0,0}=-4.925$ ($-$4.577), $\varepsilon_{0,1}=-3.871$ ($-$3.831), $\varepsilon_{0,2}^{\rm HO,LU}=-2.595$ ($-$2.708) eV for Na12; and $\varepsilon_{0,0}=-5.073$ ($-$4.755), $\varepsilon_{0,1}=-4.177$ ($-$4.135), $\varepsilon_{0,2}^{\rm HO}=-3.119$ ($-$3.189), $\varepsilon_{1,0}^{\rm LU}=-2.787$ ($-$2.048) eV for Na18.

Table IV. Computed quantities ${\rm IP}$ and ${\rm EA}$ (in eV) for Al clusters. Experimental values of ${\rm IP}_{\rm exp}$ (Al20-63 – 38 , Al89,90 – Kresin-2015 ) and ${\rm EA}_{\rm exp}$ (Al20-23 – 38 ).

[TABLE]

With increasing $N$, the contribution from the cluster bulk becomes more and more important and for large enough $N$ the point and the circles are interchanged, that is, the difference $\Delta({\rm IP}_{N\rightarrow\infty})$ becomes negative by the sign.

Self-consistent values of IP, EA, $\varepsilon^{\rm HO}$ and $\varepsilon^{\rm LU}$, which are found according to the general formulas (2), allow us to evaluate capacitances of charged and neutral clusters using the expressions (52)

[TABLE]

[TABLE]

Similar formulas for $\mathcal{C}_{N}$ correspond to a defect-free clusters.

In classical electrostatics, the capacitance of conducting sphere is determined by its radii $R_{N,v}$. Surface roughness on an atomic scale (atoms have a finite volume) does not allow to determine exactly the boundary Pog-94 ; P-1994 . In the jellium, the boundary of the ion core always corresponds to the coordinate $r=R_{N,v}$. However, the electronic cloud is more and more ‘‘splashes out’’ beyond the boundary of the core as its radius $R_{N,v}$ decreases. Moreover, such a ‘‘splashing’’ depends on the sign of the excess cluster charge (Fig. 3). As a result, the quantities $\mathcal{C}_{N,v}^{+}$, $\mathcal{C}_{N,v}^{-}$ and $\mathcal{C}_{N,v}$ are equal to each other only in the limit $N\rightarrow\infty$.

Fig. 6 shows the results of calculations for capacitances $\mathcal{C}_{N}$ and $\mathcal{C}_{N,v}$ normalized by their radius $R_{N}$ and $R_{N,v}$ (atomic units), respectively. Sign-alternating difference

[TABLE]

Is observed for certain intervals of $N$, which correspond to filling of $s-$ and $p-$ electronic shells. The difference $\Delta\tilde{\mathcal{C}}_{N}$ is determined mainly by the ratio of the quantities $\varepsilon^{\rm HO}$ and $\varepsilon^{\rm LU}$ for different $l$ in perfect and defective clusters, which can vary depending on the radial quantum number. The capacitance of defective clusters at filling shells with small $l$ is larger than for perfect ones, and the opposite relation is observed for large $l$.

Using the experimental values of ${\rm IP}_{1}$ and ${\rm EA}_{1}$ for the Na atom ($R_{1}=r_{0}$), and also the condition $\varepsilon^{\rm HO}_{1}=\varepsilon^{\rm LU}_{1}$ for unfilled shells, as a test, we get the value $\mathcal{C}_{1}/r_{0}=1.8$. It is in a satisfactory agreement with the calculated values for the smallest clusters. For non-closed electronic shells, the cluster may have a lower symmetry (spheroidal or triaxially deformed droplet Landman-95 ).

The charging effect expressed in the electric capacitance of the cluster anions and cations depends on the sign of the excess charge. Normalization of the capacitance allows to give a simple interpretation of the results of calculations: an excess negative charge leads to an effective increase in the electron cloud (effective radius) of the cluster, and excess positive charge results in a decrease of both effective radius and capacitance. This is qualitatively confirmed by the behavior of electronic profiles in Fig. 3.

In Jarrold-2009 , in addition to measurements of thermal capacity of cluster anions and cations Al35-70, the ionization potentials and electron affinity were calculated. Calculations are carried out by the DFT under the condition of a global minimum of the total energy for various configurations of atoms. The results (in eV) presented in Fig. 9 in Jarrold-2009 are approximated by us in the form

[TABLE]

[TABLE]

It follows from (9) that ${\rm IP}_{N}$ and ${\rm EA}_{N}$ must approach the same quantity $W_{0}$ with increase of $N$. Such a tendency does not agree with (54) and (55). Also the comparison of (54) and (55) from one side and (9) from another side allows to extract values of $\widetilde{\mu}_{1}$. The results are ambiguous: it follows from (54) that $\widetilde{\mu}_{1}\approx+0.58$ eV, while $\widetilde{\mu}_{1}\approx-1.51$ eV from (55). Let us remind that $\widetilde{\mu}_{1}\approx+0.68$ eV for Na Seidl-1998 ; pppppp .

III Dissociation and cohesion energies

The dissociation energy of a neutral metallic (Me) cluster according to the reaction Me${}_{N}\rightarrow$Me${}_{N-1}+$ Meat is determined by the difference in total energies

[TABLE]

In the stabilized jellium, the energy of the atom $E_{\rm at}$ is the total energy of a metal sphere of radius $r_{0}$.

By definition, the cohesive energy $\varepsilon^{\rm coh}_{N}$ is the binding energy (atoms in the cluster) per atom. It is determined by the difference in the total energy of the $N$ free atoms and the energy of a cluster consisting of $N$ atoms

[TABLE]

For $N\rightarrow\infty$, $\varepsilon^{\rm coh}_{N}\rightarrow\varepsilon^{\rm coh}_{\infty}\equiv\varepsilon^{\rm coh}(r_{0})$. The calculated values of $\varepsilon^{\rm coh}(r_{0})=$ 1.16, 1.17 and 3.97 eV, respectively for Na, Mg and Al, are in satisfactory agreement with the experimental values of $\varepsilon^{\rm coh}_{\infty}=$ 1.11, 1.51 and 3.39 eV (see Ref. Ziesche and references therein).

The binding equation has the form

[TABLE]

Asymptotics of the size dependence of the cohesion energy (57) is well known 179

[TABLE]

where the last term can be written as $-Z\mu_{1}/R_{N}$.

It should be noted that even in the works of Frenkel and Langmuir it is noted that for some substances at low temperatures a universal ratio is observed

[TABLE]

which is constructed from the observed values: the average distance between atoms $r_{0}$, the surface energy $\sigma$, and the heat evaporation $q=\varepsilon^{\rm coh}(r_{0})$ (see Table 2 in Pog-94 ). Using this relation, the asymptotics (59) can be rewritten in a form convenient for estimations

[TABLE]

Next, using Eqs. (59) and (56), we find a coincidence of the asymptotics $\varepsilon^{\rm coh}_{N}$ (59) and $\varepsilon^{\rm dis}_{N}$.

It is of interest to determine the effect of charging on dissociation and cohesion energies of clusters. The previous results in electronic model and molecular dynamics simulations show that fragmentation consists mainly in emission of single atoms Bennemann . Due to the fact that the ion work function is greater than the heat of evaporation of the neutral atom PogKurbVas-2005 , using Eqs. (56), (57) and expressions

[TABLE]

the energy differences can be expressed as

[TABLE]

Figs. 7 and 8 show the dissociation and cohesion energies of perfect clusters and clusters containing monovacancy. Size dependence of the dissociation energy in Fig. 7 is represented by quantum oscillations around the asymptotical dependence. The values $\varepsilon^{\rm dis}_{N,v}$ for defective cluster, with large $l$, are larger than for a perfect one, and for small $l$, they interchange. For perfect and defective clusters, in addition to a change of order filling of electronic levels, a significant difference in the behavior of dissociation energy with the same number of atoms $N$ is also seen: for small $l$, dissociation energy of perfect clusters decreases with increasing $N$, while it increases for defective clusters. Comparison of the data in Fig. 8 and 9 confirms the accuracy of formula (58), and also explains the difference in the position of local maxima. The calculated values of $\varepsilon^{\rm coh}_{N,v}$ are closer to the experimental values, obtained for $T=150$ K Martin , than $\varepsilon^{\rm coh}_{N}$. We note that near the phase transition the quantities $\varepsilon^{\rm dis}_{N}$, extracted from measurements of melting temperature and latent heat of the transition can be negative Jarrold-2010 .

Thus, we can conclude that the most stable defect-free clusters are those, which have latest filled levels with a small $l$, and for defective one - on the contrary. In the experiments, the size oscillations $\varepsilon^{\rm dis}_{N}$, are apparently suppressed by temperature effects (see Fig. 9 in Brechignac ).

As shown in Fig. 8 cohesion energy of cluster anions and cations is different from neutral clusters. Excess positive charge leads to a decrease in the energy of cohesion due to the increase of the forces of electrostatic repulsion, while excess negative charge leads to the opposite effect. A behavior of dissociation and cohesion energies of charged clusters is described by the Eq. (61).

IV Monovacancy-formation energy

A considerable number of papers deal with first principles (ab initio) calculations of the monovacancy-formation energies in metals Freysoldt-2014 . In the stabilized jellium and the liquid drop models the energy of cohesion of an atom and the monovacancy-formation energies are investigated in form of the Padé expansion Ziesche (see also Pog-94 ). Within the notation of the work Pog-94 , the results of Ziesche can be expressed as

[TABLE]

where $\tilde{\delta}_{1}\equiv\delta_{1}/r_{0}$ and $\tilde{\delta}_{2}\equiv\delta_{2}/r_{0}^{2}$.

Calculated in BVP-2014-2 $\varepsilon_{\infty}^{\rm vac}\equiv\varepsilon^{\rm vac}(r_{0})=$ 0.33, 0.72 and 1.00 eV agree with the experimental values of 0.335, 0.84 and 0.73 eV for Na, Mg and Al, respectively Kraftmakher .

Using values $\varepsilon^{\rm coh}(r_{0})$ and $\varepsilon^{\rm vac}(r_{0})$, and $\tilde{\delta}_{2}=-0.13$ (Na), $-0.015$ (Mg) and +0.22 (Al) from Ziesche , we find $\tilde{\delta}_{1}=0.32$ (Na), 0.54 (Mg) and 0.57 (Al). The values $\delta_{1}$ and $\delta_{2}$ are necessary for us in order to construct the asymptotics of the monovacancy-formation energies.

For clusters, self-consistent calculations of $\varepsilon^{\rm vac}_{N,v}$ have not been performed due to the need for detailed description of the vacancy formation process. Therefore, it is of interest to elucidate, which of the two scenarios of the vacancy formation is favorable.

Within the Schottky scenario, an atom is extracted from the surface of a perfect sphere, and in the final state the vacancy is in the center of the sphere. In this case, we have

[TABLE]

where $E_{N-1,\,v}$ is the total energy of sphere with vacancy (the spherical layer between $r=r_{0}$ and $r=R_{N-1,\,v}$ contains $N-1$ atoms).

According to another scenario Pog-94 , the number of atoms does not change, but the ‘‘bubble’’ of radius $r_{0}$ is blown in the center of the system. In this case, we have

[TABLE]

The comparison of Eqs. (63) and (64) demonstrates the advantage of the second mechanism by relation

[TABLE]

We construct the asymptotics of the monovacancy-formation energy. Its size dependence is determined by the difference in the total energies of the spheres calculated by the formulas (63) and (64) in the limit $N\rightarrow\infty$, and reduces to the difference of total surface energies.

For the vacancy blowing mechanism, using (5), we obtain

[TABLE]

or the Schottky mechanism, in accordance with (65) and $R_{N-1,v}=R_{N}$, the asymptotics are determined by the sum of the expressions (66) and (59). Asymptotic dependence $\varepsilon_{N,v}^{\rm vac,Sh}$ weakly depends on $N$, and the dependence (66) demonstrates a decrease in vacancy formation energy with an increase of $N$, which agrees with the conclusions of the paper Sinder ; Delavari , but it contradicts the conclusions of Refs. Yang-2007 ; Hendy ; Guisbiers .

The asymptotic size behavior of the vacancy formation energy (66) can be qualitatively compared with the results of Sinder in which the energies of the cluster with or without vacancy have been calculated on the basis of the tight-binding approximation. Representing the expression (66) in the form

[TABLE]

we obtained that the values $C=1.21,1.50$ and 1.03 for Na, Mg and Al, respectively, are qualitatively consistent to $C=1.33$ and 1.46 for Cu and $\beta-$Ti from Sinder .

Fig. 9 shows the results of calculations of the monovacancy-formation energy by two mechanisms. These calculations confirm the formula (65), namely, the advantage of blowing a vacancy. All dependences experience strong fluctuations. For some $N$, especially for Al, the values $\varepsilon_{N,v}^{\rm vac,blow}$ become negative in narrow ranges of $N$. Such areas are shown in Fig. 4 and above, the paper already contains comments on hierarchy of electronic states in such clusters.

Table V. Monovacancy-formation energies (in eV) for neutral NaN.

[TABLE]

The paper Itoh reports on the results of ab initio calculations on NaN=55,147,309 clusters that show icosahedral growth. The monovacancy-formation energy in cluster depends on the site at which the vacancy is created. A vacancy at the center or the first atom shell is found to cost much higher energy compared to other sites in the clusters.

In Table V, we compared the vacancy formation energy according to the Shottky mechanism (the atom is removed from the cluster center to outside) and according to the blowing bubble mechanism (the removed atom is placed at the flat icosahedra surface Itoh ). On one hand, the comparison of $\varepsilon_{N,v}^{\rm vac,Sh}$ (is denoted by $E_{r}$ in Ref. Itoh ) reveals a difference of about 0.4 – 0.5 eV, but, on the other hand, it shows a similar size dependence. Values of $\varepsilon_{N,v}^{\rm vac,blow}$ (is denoted by $E_{v}$ in Ref. Itoh ) agree much better. In Table 5, in order to demonstrate the dependence of $\varepsilon_{N,v}^{\rm vac,blow}$ on the location of vacancy formation, two values are presented, which are separated by the slash. They correspond to the atom displacement from the cluster center and from the first icosahedral atomic shell to the surface. The major aim of our simple treatment is not on the absolute values of clusters characteristics, but on their change upon vacancy formation.

The difference between the energies of formation of vacancies in a charged and neutral cluster in accordance with two scenarios ($E_{N}^{\pm}\rightarrow E^{\pm}_{N-1,\,v}+E_{\rm at}$ and $E_{N}^{\pm}\rightarrow E^{\pm}_{N,\,v}$) by analogy with (61) can be represented in the form of relations

[TABLE]

The character of the size dependence $\varepsilon_{N,v}^{\rm vac,\pm}$ (67) is fully supported by direct calculations and is determined by the behavior of IP and EA in Fig. 4.

The above calculations in LDA correspond to zero temperature. Perhaps, for a density of atoms corresponding to finite temperature, reduced symmetry of clusters, as well as the use of LSDA for exchange-correlation energy, strong oscillatory behavior of the energy characteristics will be suppressed. Anyway, the size behavior of the results of direct calculations agrees with its asymptotics.

In quasi-thermodynamics, the probability of the appearance of vacancies in a cluster at finite temperature $T$, can be estimated from the condition of free energy variation $\Delta F_{N,v}^{\rm vac,blow}$,

[TABLE]

As a result of the fact that when the vacancy is blown, the number of ions in the cluster does not change, the entropy contribution is provided only by the degenerate electron gas A-78 . The corresponding expression is

[TABLE]

For calculations in (69), equilibrium profiles of electron distributions in the stabilized jellium will be needed for given $N$ and $T$. At zero temperature and $N=12$ these profiles are presented on Fig. 3.

V Summary

In this paper, we propose a consistent formal procedure for finding ionization potential of a large metal cluster containing the vacancies that is based on a previously solved problem on the scattering of electrons by a monovacancy in bulk metal by the Kohn-Sham method in the stabilized jellium. Self-consistent profiles are used to determine the vacancy energy shift of the ground state in a metal cluster sphere as a series of size corrections. Spherical periodicity in the arrangement vacancies was assumed. The limits of applicability of this expansion to powers of the inverse radius are $R>4.5$ nm and $R>6$ nm for Na and Al, respectively.

In our work, the effect of bulk vacancies on the properties of a cluster is estimated. Only by introducing periodic arrangement of internal vacancies it is possible to consider the defected sample as a cellular media and to estimate the shift of its ionization potential. The presented approach seems to be promising for experimental estimation of the concentrations of bulk point defects or impurities in metal clusters. To this end, it is first needed to calculate the scattering length of electrons on the corresponding defect in 3D metal. The obtained analytical expressions are convenient for analysis of the results of photoionization experiments. In particular, the concentration of internal vacancies in a cluster formally can be estimated near the melting point using these expressions. The jellium and ab initio calculations Landman-95 ; Akola-2000 ; Kanhere ; Itoh reveal a rather subtle interplay between geometric and electronic shell effects, and evidences that the quantum mechanical description of the metallic bonding is crucial for understanding quantitatively the variation in melting temperatures observed experimentally for free clusters.

There is no doubt that the site of vacancies in the bulk, near or on the surface Itoh will affect the characteristics of the cluster differently. And without a doubt it is difficult to determine their contribution separately. However, due to the fact that the concentration of equilibrium internal vacancies in 3D metals exponentially depends on the inverse temperature, the surface energy and work function will contain, in addition to the linear temperature dependences associated with thermal expansion, also the exponential dependence which is easier to observe near the melting of 3D metals. Due to the fact that the energy of the vacancy formation depends on its site in the cluster (on the surface it is approximately half the bulk), perhaps we should expect the weak exponential temperature dependences of their concentrations and, correspondingly, exponentional temperature-dependent contribution to both IP and EA. The ability to distinguish and separate these dependences for large clusters or island films depends of course on the accuracy of the experiment.

The Fowler-Nordheim theory is applicable for planar field electron emission. As a modification, we can propose the use of the temperature and vacancy dependence of the electron work function. For special electrode geometry, adaptation of the theory is also possible in an analytical form (see, for example, Yuasa ).

The self-consistent calculations of the electronic profiles for perfect clusters and clusters with a vacancy allowed us to determine the total energy of the neutral and charged defective cluster and then calculate the dissociation, cohesion, vacancy-formation energies, the ionization potential, and electron affinity as well as the electrical capacitance.

The results of calculations for Na and Al are compared with the asymptotics and results for defect-free clusters. The ionization potential for the smallest cluster with a vacancy is greater than for a perfect cluster (approximately 0.1 eV for Na and 0.5 eV for Al). A maximum difference is observed at the transition from a completely filled shell to the empty one. As $N$ increases, this difference disappears.

Magic numbers of atoms for perfect clusters and clusters with monovacancy are different, especially for Al. Normalized electrical capacities of clusters always exceed unity and contain quantum size fluctuations. In this case, for defective clusters with partially filled electronic shells capacities are significantly larger than for perfect clusters.

The size dependence of the cohesion energy contains local maxima. Clusters corresponding to them are more stable, that is, have large-scale binding, dissociation, and vacancy-formation energies than their neighbors. For small clusters, such maxima appear at the complete filling of the next electron shell. The positions of the maxima for defect and defect-free clusters are different, which is due not only to the difference in their sizes, but also due to the character of the behavior of the electron wave functions.

The energy of cohesion of charged cluster anions and cations is different from the cohesive energy of neutral clusters. Excessive positive charge leads to a decrease in energy through increased forces of electrostatic repulsion, and excess negative charge leads to the opposite effect.

The quantum-size dependences of the vacancy-formation energies in the Schottky and the ‘‘bubble blowing’’ scenarios, and their asymptotic tendencies were determined. Strong size fluctuations in the entire cluster size range were found. Size asymptotics for these two mechanisms are different from each other and are weakly dependent on the number of atoms in the cluster. The nature of the size dependence of vacancy-formation energies from excess charge in a cluster is determined by the behavior of the ionization potential of the cluster and the electron affinity.

With the increase in $N$ the dissociation energy is either increases, or has a local minimum (in the areas between the maxima), while the vacancy-formation energies decrease monotonically.

Figs. 7 – 9 indicate that all characteristics of charged and neutral clusters differ from each other. It is reasonable to assume that charging might control the melting temperature of the clusters (see experimental size dependent melting temperature of anions and cations Al35-70 at Fig. 4 in Jarrold-2009 ).

In this model, the relaxation of the cluster volume was not taken into account. In the limit of large clusters, the effect of self-compression on the ionization potential is analytically described in iakPog-1995 ; keiPog-1996 , and for small clusters is numerically investigated in a spherically averaged pseudopotential model Vieira . Relaxation of ionic distribution in the cluster will lead to a decrease in the total energy. More coherent, but also more labor-intensive, are ab initio methods with the selection of the coordinates of ions under the minimum condition of the total energy of the cluster. Such a procedure is implemented in Jarrold-2009 only for Al30-70 clusters, i.e., clusters with a small atoms number.

Acknowledgements.

We are grateful to Walter V. Pogosov for reading the manuscript.

Appendix A Vacancy quantum correction

The straightforward estimation for $\left\langle\delta V\right\rangle_{R_{N,v}}$ in (44) gives the surprising result:

[TABLE]

where $\langle...\rangle$ denotes the integration over the cluster volume. It seems as if $\left\langle\delta V\right\rangle_{R_{N,v}}$ is proportional to $(R_{v}R_{N,v})^{-1}$. In this case the hierarchy of terms in our expansion (42) would be broken, because the previous term is proportional to $R_{N,v}^{-2}$. However, below we shall demonstrate the emergence of an extra factor to $\left\langle\delta V\right\rangle_{R_{N,v}}$ which is proportional to $\xi\left(R_{v}/R_{N,v}\right)$,

[TABLE]

It is a result of the integration in Eq. (44) over angles PR-2017 .

Due to the fact that the perturbation $\delta V(r)$ occurs on a scale of supercell, in (44), it will be reasonable to proceed to integration over supercell and use the Green formula:

[TABLE]

where $\Omega_{i}$ is a volume of space under condition $L_{v}<\varrho<R_{v}$.

Let us introduce the gradient expansion of the squared wave function $\psi^{2}(r)$ near the point $\textbf{R}_{i}$,

[TABLE]

Using expansions of derivatives over the small parameter $\varrho/R_{i}\equiv R_{v}/R_{i}$ we can take the derivatives over r at the centers of cells $\textbf{R}_{i}$,

[TABLE]

and

[TABLE]

Then

[TABLE]

Using expression (16), we obtain

[TABLE]

The terms under the sum in Eq. (76) have the opposite signs.

The integral over the surface of ith cell is evaluated exactly (see Fig. 10),

[TABLE]

Hence, the integration over angles at $R_{i}\gg R_{v}$ gives the additional power of $R_{v}/R_{i}$. Now, taking account for the factor $R_{v}/R_{i}$ from Eq. (78), we replace in Eq. (76) summation by integration and obtain

[TABLE]

[TABLE]

Next, after the change of the summation over cell numbers to the integration over the cluster volume it is possible to reduce the second term in (76) to a surface integral. The integral vanishes due to the boundary condition (41). Finally, using Eqs. (78) and (79), we have first nonzero term

[TABLE]

where $D_{1}=12D_{0}/\pi$. It is worth noting that the excluded volumes inside the supercells contribute to neglected terms $\sim\xi^{3}$.