Chemical Potential of Integer Electron Systems

TL;DR

This paper introduces a density embedding method to analyze how finite-distance interactions between an atom and a metal surface affect the atom's electron count and chemical potential, revealing smoothed electron number steps similar to temperature effects.

Contribution

The paper presents a novel density embedding approach that accounts for finite-distance interactions, reducing the range of chemical potentials leading to integer electron counts.

Findings

Finite-distance interactions smooth out the N(μ) staircase.

Fractional occupations occur only at sharply-defined μ.

Method demonstrated on a model atom-metal surface system.

Abstract

A truly isolated atom always has an integer number of electrons. If placed in contact with a far-away metallic reservoir, a {\em range} of metallic chemical potentials $\mu$ will lead to an identical number of electrons, $N$, on the atom. We formulate a density embedding method in which the range of $\mu$ leading to integer $N$ decreases due to finite-distance interactions between the metal and the atom. The typical $N(\mu)$ staircase function is smoothed out due to these finite-distance interactions, resembling finite-temperature effects. Fractional occupations on the atom occur only for sharply-defined $\mu$'s. We illustrate the new method with the simplest model system designed to mimic an atom near a metal surface. Because calculating fractional charges is important in various fields, from electrolysis to catalysis, solar cells and organic electronics, we anticipate several potential uses of the proposed approach.

Equations

$$E(N)=(1-\omega )E(m)+\omega E(m+1),$$

$$-I<\mu <-A,$$

$$v(\mathbf{r})=v_{atom}(\mathbf{r})+v_{metal}(\mathbf{r}),$$

$${\mu }_{atom}={\mu }_{metal}=\mu ,$$

$$-I_{atom}<{\mu }_{metal}=\mu <-A_{atom},$$

$$v_{metal}(x)=\frac{-V_0}{1+e^{s(x-R)}},$$

$$v_{atom}(x)=-Z{\cosh }^{-2}(\gamma x),$$

$$n(x)=\frac{1}{2\pi i}{\int }_{C_{\mu }}{G}_E(x,x)dE.$$

$$n_{atom}(x)=\omega n_{p+1}(x)+(1-\omega )n_p(x),$$

$$\varepsilon_{i}^{(0)}=-\frac{\gamma^{2}}{8}\Big{\{}-[1+2(i-1)]+\sqrt{1+(8Z\gamma^{2})}\Big{\}}^{2}~{}~{},$$

$$\omega^{(0)}=\int dx\Big{[}n(x)-n_{metal}^{(0)}(x)-n_{core}^{(0)}(x)\Big{]}~{}~{},$$

$${\chi }_{\alpha }^{(i)}(x,x^{\prime })=\frac{d{n}_{\alpha }^{(i)}(x)}{d{v}_{eff,\alpha }^{(i)}(x^{\prime })},$$

$$N_{atom}^{(i+1)}=N_{atom}^{(i)}+\frac{{\mu }_{metal}^{(i)}-{\mu }_{atom}^{(i)}}{d{\mu }_{atom}^{(i)}/d{N}_{atom}^{(i)}},$$

$$\frac{d{\mu }_{atom}^{(i)}}{d{N}_{atom}^{(i)}}=\int dx{n}_{\mathrm{HOMO}}^{(i)}(x)\frac{d{v}_R^{(i)}(x)}{d{N}_{atom}^{(i)}}.$$

$$\begin{split}\frac{dv_{R}^{(i)}(x)}{dN_{atom}^{(i)}}&=\int dx^{\prime}n_{\rm HOMO}^{(i)}(x)\Big{[}\Big{(}\chi_{metal}^{(i)}(x,x^{\prime})\Big{)}^{-1}+\\ &\Big{(}\chi_{atom}^{(i)}(x,x^{\prime})\Big{)}^{-1}\Big{]}~{}~{}.\end{split}$$

$$\begin{array}{rl}{E}_{atom}& \equiv \sum _{i=1}^{p_{atom}}{\epsilon }_{atom}^i+w{\epsilon }_{atom}^{\mathrm{LUMO}}\\ & -\int dx{v}_R(x)n_{atom}(x).\end{array}$$

$$v_R(x)=v_{R,\mathrm{ext}}(x)+v_{R,\mathrm{kin}}(x),$$

$$\frac{n_{metal}(x)}{n(x)}{v}_{atom}(x)+\frac{n_{atom}(x)}{n(x)}{v}_{metal}(x),$$

$$\bar{n}=\frac{1}{e^{({\epsilon }_i^{(0)}-\mu )/kT}+1}$$

Full text

1 Introduction

When the Hohenberg-Kohn theorem [1] was extended to fractional electron numbers by Perdew, Parr, Levy, and Balduz (PPLB, [2, 3, 4]), a result of far-reaching consequences was found: The ground-state energy of an $N$-electron system, $E(N)$, is a piecewise-continuous linear function of $N$:

[TABLE]

where $E(m)$ and $E(m+1)$ are the ground-state energies of the (integer) $m$ and $(m+1)$-electron systems, and $0\leq\omega\leq 1$. At strictly zero temperature, an atom or molecule that is in equilibrium with a far-away metal reservoir, will be neutral in the ground state for any chemical potential $\mu$ in the range

[TABLE]

where $I$ is the ionization potential and $A$ the (positive) electron affinity of the neutral atom. For chemical potentials lower than $-I$, the atom transfers one electron to the reservoir. For chemical potentials higher than $-A$, the atom receives one electron from the reservoir. The number of electrons in the atom is thus a staircase function of the chemical potential (black dot-dash line in Fig. 1), which is clearly only sharply defined for non-integer numbers. The range of $\mu$ that is consistent with the integer $m$ is the fundamental energy gap of the atom, $E_{g}=I-A$, which is thus given by the total discontinuity in the derivative of $E(N)$ with respect to $N$ at $N=m$. All properties of the system involving derivatives of the energy with respect to $N$ are similarly undefined at the integers at zero temperature.

A smoothening of the discontinuities at integer numbers of electrons and a range of $\mu$ that is narrower than $I-A$ can be found by applying techniques of the grand-canonical ensemble at finite temperature [5]. The main result of our work is that sharper values of $\mu$ can be found even at zero temperature by considering finite distances from the metal reservoir. To show this, an unambiguous definition is needed for the charge of the atom when it is located at an interacting distance from the metal. We provide such definition by requiring that the chemical potential of the two fragments (metal and atom) be equal while satisfying the standard constraint of density-embedding methods, i.e. that the sum of the two fragment densities be equal to the total electronic density. With this definition of fragments, the regions of strictly integer numbers of electrons on the atom are narrower than $I-A$ when the atom is at an interacting distance from the reservoir (red and blue lines in Fig. 1). Outside of the shaded regions in Fig. 1, the atom acquires a fractional number of electrons. At large separations between the atom and the metal, our model recovers the PPLB results. At shorter separations, the regions of integer occupations shrink but do not collapse to a single point. Due to the finite-distance interactions, the effective values of $I$ and $A$ are different from those of the isolated atom. As a result, the narrowing of the integer windows is not symmetric with respect to ${(I+A)}/{2}$ and is markedly different near different integer occupations.

The method is described in Sec. 2 and illustrated through explicit numerical computation in Sec. 3. We end with a brief summary and outlook in Sec. 4.

2 Chemical-potential constrained Partition-DFT

Consider a system of electrons in an external potential $v({\bf r})$ that can be written as:

[TABLE]

where $v_{metal}({\bf r})$ describes a background periodic or semi-periodic metallic potential supporting a continuum of electronic levels occupied up to a Fermi energy, $\epsilon_{F}$, and $v_{atom}({\bf r})$ is a localized potential such as the Coulomb or screened-Coulomb potential of an atom. The partition of Eq.(3) is useful when one wants to describe an atomic defect in a solid or an atom adsorbed on a metal surface.

The task of finding the number of electrons on the atom, $N_{atom}$, is non-trivial unless $v_{atom}({\bf r})$ is far from all regions where $v_{metal}({\bf r})$ is non-zero, in which case one recovers the black staircase function of Fig.1 with $\mu=\epsilon_{F}$. The total density $n({\bf r})$ for the combined system of atom and metal can be partitioned as $n_{atom}({\bf r})+n_{metal}({\bf r})$ in many different ways. Partition Density Functional Theory (P-DFT, [6, 7, 8]) provides an elegant, unambiguous method for performing such a partition when the number of electrons is finite and the external potential for each fragment vanishes in all directions as $|{\bf r}|\to\infty$. Fragments in P-DFT are isolated from each other and are in contact with a far-away electronic reservoir through which they can exchange electrons. The interaction energy between the fragments is recovered by means of a unique global embedding potential, referred to here as the reactivity potential, $v_{R}({\bf r})$. The prescription to determine $N_{atom}$ becomes simple: Minimize the sum of the fragment energies (i.e. atom and metal) subject to the constraint that the fragment densities sum to to the correct total density, and then calculate the number of electrons in the atom as $N_{atom}=\int n_{atom}({\bf r})d{\bf r}$. This number is in general not an integer because each P-DFT fragment energy is given by the ensemble expression of Eq.(1), where the non-integer $\omega$ is one of the parameters to be optimized during the energy minimization.

In the case of the potential of Eq.(3), however, $v_{metal}({\bf r})$ does not vanish as $|{\bf r}|\to\infty$ in all directions, and one of the fragment energies is infinite. The approach of P-DFT is thus not directly applicable.

In lieu of an energy minimization, we propose here to impose a chemical-potential equalization constraint, shown to be equivalent to energy-minimization for the case of finite systems [9]. The prescription is just as simple: Find the fragment densities that equalize the chemical potentials of the fragments and the chemical potential of the combined system:

[TABLE]

while adding to the correct total density. The resulting density of the atom is an ensemble ground-state density of $v_{atom}$ that is modified by the addition of $v_{R}({\bf r})$, which is identical for both atomic and metallic fragments.

When $N_{atom}$ is an integer, $\mu_{atom}$ is defined only within a range, so Eq. 4 is applicable only for non-integer values of $N_{atom}$. For integer occupations, the condition of Eq. 4 is modified taking into account Eq. 2:

[TABLE]

where $I_{atom}$ and $A_{atom}$ are computed in the presence of $v_{R}({\bf r})$. We consider our method converged if either $N_{atom}$ is non-integer and condition 4 is satisfied or if $N_{atom}$ is integer and condition 5 is satisfied. In the following section, we successfully apply this method to a model system that mimics an atom-metal interface in 1-D; however, the rigorous derivation of the conditions for the existence of a unique reactivity potential for systems with semi-infinite fragments is still not established.

3 Simple Illustration

We choose the simplest non-trivial system that exhibits the features we need: One semi-infinite fragment (the ‘metal’) and one finite fragment with a small number of bound states (the ‘atom’). The total number of electrons is infinite, but the electrons are non-interacting and restricted to move in only one dimension.

3.1 Model System

The metal is represented by a potential that goes to a negative constant $-V_{0}$ as $x\to-\infty$:

[TABLE]

and is populated with non-interacting ‘spinless electrons’ up to the Fermi level $\epsilon_{F}$, with $-V_{0}<\epsilon_{F}<0$. In Eq. 6, $R$ is the separation between the metal surface and the center of the atomic potential, and $s$ is a parameter that determines the steepness of the step. The form of the potential allows it to be smooth enough to be used with finite-difference methods on a spatial grid while preserving a steep step-like feature. The atom is represented by a finite potential with a finite number of bound states:

[TABLE]

where $Z$ and $\gamma$ are parameters that control the depth and width of the well. We use $V_{0}=3.5$, $Z=2$, and $\gamma=0.5$ throughout the paper. The total external potential is then just the sum of $v_{metal}$ and $v_{atom}$ according to Eq.(3), as shown in Fig. 2.

The full system of metal plus atom produces a continuum of states. The chemical potential of this system is its Fermi energy. The reactivity potential $v_{R}(x)$ ensures fragment densities sum to the total density of the system. The densities of the total system and of the metal fragment are calculated as the integral [10, 11]:

[TABLE]

Here, $G_{E}(x,x^{\prime})$ is the corresponding Green’s function that can be found numerically exactly. The integral is evaluated over $C_{\mu}$, a contour in the complex energy plane containing all possible occupied states. The total system density $n(x)$ for a large separation $R=15$ is shown in Fig. 3. We can see that as the chemical potential of the system increases through the energy levels of the isolated atomic potential, the density near the atom increases in large jumps every time the chemical potential reaches a bound state. The atomic densities have the ensemble form [2]:

[TABLE]

where $p$ is the lower bounding integer of $N_{atom}$, $0\leq\omega<1$, and $N_{atom}=p+\omega$. Calculations of the atomic densities at integer occupations are trivial.

The eigenvalues of $v_{atom}(x)$, when isolated, are known for any chosen $Z$ and $\gamma$ [12]:

[TABLE]

where $i$ runs from 0 to the maximum number of bound states. The superscript ‘$(0)$’ indicates that the atom does not interact with the metal. These eigenvalues broaden into resonances when the atom couples to the metal and $|\epsilon_{F}|>|{\varepsilon_{i}}|$.

3.2 Search for chemical-potential Equalization

To obtain a single point on the $N_{atom}$ versus $\mu$ plot in Fig. 1, we perform a numerical algorithm for a set value of $\mu$. This algorithm consists of the ‘inner’ inversion that computes the reactivity potential at the current guess of $N_{atom}$ and the ‘outer’ loop that updates $N_{atom}$ until one of the chemical-potential equalization conditions, Eq. 4 or Eq. 5, is satisfied. Our inversion method requires the precomputed total density $n(x)$ for each $\mu$. We set $\mu_{metal}$ equal to $\mu$ and do not vary it throughout the inversion procedure.

We choose $v_{R,Guess}^{(0)}(x)=0$ as our initial guess for $v_{R}(x)$. To calculate the initial guess for $N_{atom}$, we start by calculating the isolated atomic density, $n_{atom}^{(0)}(x)$. We separate $n_{atom}^{(0)}(x)$ into contributions from the density of the highest occupied atomic orbital (HOMO) and the density due to the core electrons, $n_{atom}^{(0)}(x)=n_{core}^{(0)}(x)+n_{\rm HOMO}^{(0)}(x)$. The number of states included in the core region, $N_{core}$, is equal to the number of eigenvalues of the isolated atom which are below $\mu$. The initial guess for the atomic occupation is $N_{atom}^{(0)}=N_{core}+\omega^{(0)}$, where $\omega^{(0)}$ is found as:

[TABLE]

and $n_{metal}^{(0)}(x)$ is the density of the isolated metal fragment.

At each iteration $i\geq 0$ of the ‘outer’ loop, we use the current values $N_{atom}^{(i)}$ and $v_{R,Guess}^{(i)}(x)$ to compute the $v_{R}^{(i)}(x)$ that minimizes the difference between $n(x)$ and $n_{f}^{(i)}(x)=n_{metal}^{(i)}(x)+n_{atom}^{(i)}(x)$ to numerical precision (i.e. $v_{R,Guess}^{(i)}(x)$ is used as an initial guess to find $v_{R}^{(i)}(x)$ at fixed $N_{atom}^{(i)}$ and $\mu_{metal}$). The resulting fragment densities are used to calculate the fragment responses $\chi_{\alpha}^{(i)}(x,x^{\prime})$ that will be used to update $v_{R,Guess}^{(i)}(x)$:

[TABLE]

where $v_{eff,\alpha}^{(i)}(x)=v_{\alpha}(x)+v_{R}^{(i)}(x)$ and $\alpha$ is either ‘metal’ or ‘atom.’

If $N_{atom}^{(i)}$ is not an integer, then $\mu_{atom}^{(i)}$ equals the HOMO energy $\varepsilon_{\rm HOMO}^{(i)}$ in the presence of $v_{R}^{(i)}(x)$, and we use Eq. 4 to check if the algorithm has converged. On the other hand, if $N_{atom}^{(i)}$ is an integer, we use the convergence criteria of Eq. 5, with $-I_{atom}=\varepsilon_{\rm HOMO}^{(i)}$ and $-A_{atom}=\varepsilon_{\rm LUMO}^{(i)}$, where $\varepsilon_{\rm LUMO}^{(i)}$ is the energy of lowest unoccupied atomic orbital in the presence of $v_{R}^{(i)}(x)$.

If neither of the conditions is met, we continue by calculating $N_{atom}^{(i+1)}$ as:

[TABLE]

where:

[TABLE]

In Eq. 14, the derivative on the right hand side is given by:

[TABLE]

Eq. 15 is also used to update the guess $v_{R,Guess}^{(i+1)}(x)=v_{R}^{(i)}(x)+dv_{R}^{(i)}(x)$.

3.3 Energies, densities, and reactivity potentials

The origin of the discontinuities of the chemical potential can be understood in terms of the atomic orbitals $\varepsilon_{i}$ (in the presence of $v_{R}(x)$). Near integer occupations, the energy of the HOMO shifts up from the left and the energy of the LUMO shifts down from the right, as we see in Fig. 4. Even for separations as small as $R=3$, levels do not equalize.

The effect of the finite-distance interactions on the energy of the atom can be seen in Fig. 5. The energy of the atom is defined as the sum of occupied orbitals minus the energy contribution from the reactivity potential:

[TABLE]

In Fig. 5, the dashed line shows the energy at large separation, $R=15$. It consists of straight line segments [2, 3, 4]. At short distances (e.g. $R=3$, solid red in Fig. 5) the line segments have a slight curvature. As shown in the inset plot of Fig. 5, the curvature is more noticeable for $N_{atom}$ in the range of $2$ to $3$, where $E_{atom}$ values are more evenly spaced. This curvature is the consequence of the inter-fragment interactions, but it does not smoothen the cusps at integer occupations.

The atomic fragment density at large values of $R$ jumps abruptly when going through integer occupations, as can been seen in the top ($N_{atom}=1$) and middle ($N_{atom}=2$) panels of Fig. 6. For each value of $N_{atom}$, increasing the Fermi energy of the system changes almost exclusively $n_{metal}(x)$. As these changes occur, we observe an increase in the value of the metal density accompanied by a decrease in the period of density (Friedel) oscillations. The bottom panel of Fig. 6 shows the representative behavior of fragment densities at small separations. Densities corresponding to non-integer values of $N_{atom}$ begin to appear. We note that, in this regime, the density of the metal fragment appears unchanged for different values of $\mu_{metal}$. The density response of the system to infinitesimal changes of $\mu$ is thus largely localized to either atom or metal fragments.

Finally, we point out a connection between the features of $v_{R}(x)$ and the known features of exact Kohn-Sham (KS) potentials, $v_{s}(x)$, at interfaces [1, 3, 13, 14, 15, 16, 17]. Since we work with non-interacting electrons, $v_{R}(x)$ has contributions only from the non-additive external potential ($v_{R,\mathrm{ext}}$) and from the non-additive kinetic potential ($v_{R,\mathrm{kin}}$):

[TABLE]

where $v_{R,\mathrm{ext}}(x)$ is given by:

[TABLE]

and $v_{R,\mathrm{kin}}(x)=v_{R}(x)-v_{R,\mathrm{ext}}(x)$. In Fig. 7, we plot $v_{R}(x)$ and its components at large inter-fragment separation when $N_{atom}=2$. We observe that $v_{R,\mathrm{ext}}(x)$ has a well in the low density region. In contrast, $v_{R,\mathrm{kin}}(x)$ has a step-like feature analogous to the feature known to be present in $v_{s}(x)$ when two inequivalent fragments are separated (note from Eqs. 17-18 that $v_{R,\mathrm{kin}}(x)=v_{s}[n_{metal}](x){n_{metal}(x)}/{n(x)}+v_{s}[n_{atom}](x){n_{atom}(x)}/{n(x)}-v_{s}[n](x)$). The magnitude of the feature in $v_{R,\mathrm{kin}}(x)$ highlights the importance of non-additive non-interacting kinetic energy functional approximations for practical embedding calculations. [18, 19]

3.4 Finite Distance as Finite Temperature

The smoothening of the $N$ vs. $\mu$ staircase in Fig.1 suggests a possible analogy between finite distances and finite temperatures. In Fig. 8, we compare our calculated $N_{atom}$ to the average number of particles $\bar{n}$ from a Fermi-Dirac (FD) distribution:

[TABLE]

where $k$ is Boltzmann’s constant and $T$ is the temperature. It is apparent from the figures that the analogy is not exact. The FD distribution at specified (unphysical) temperatures can capture some of the behavior of $N_{atom}(\mu)$ around the step between integer numbers or the upper region of the curve as it flattens near the integer. It cannot capture both at once, or correctly follow the behavior of the lower region as it rises from the lower integer.

4 Conclusions and Outlook

We have shown that the chemical potential of an integer-electron system can be smaller than $I-A$ when the system (here, an atom) is at interacting distances from a metallic reservoir of electrons. A continuous change in a global molecular property, $\mu$, distorts the density of one fragment (either metal or atom) markedly more than the density of the other fragment. The typical $N_{atom}$ vs. $\mu$ staircase function is smoothed-out as a result of the finite-distance interactions between the atom and the metal surface. Our method is useful for calculations on semi-infinite systems and allows treatment of different fragments with different computational techniques. For example, an atomic or a molecular fragment can be treated with an accurate wave-function method and the semi-infinite metal fragment can be treated with a more innate Green’s function method. The method provides a convenient way to account for the finite-distance interactions near the metal surface.

In the extension of Frozen-density embedding [20] to fragments with non-integer particle numbers [21], the total energy is minimized under the constraint that each fragment density integrates to a pre-established fractional value. In this method, each fragment has a different chemical potential along with a different embedding potential, and the fractional charges on the fragments are not an output but an input for the calculation. As an alternative, we have proposed chemical-potential equalization as the main criterion for determining fractional charges in density embedding. Because calculating fractional charges is important in various fields, from electrolysis [22, 23] to catalysis [24], solar cells and organic electronics [25, 26], we anticipate several potential uses of the proposed approach.

Next steps include the investigation of exchange-correlation effects, the application of chemical-potential constrained P-DFT to realistic systems of atoms and molecules adsorbed on metal surfaces, and the calculation of surface resonance lifetimes through complex scaling [27].