900

# Tools for Predicting the Properties of Nanomaterials

Authored by: James R. Chelikowsky

# 1 Handbook of Nanophysics

Print publication date:  September  2010
Online publication date:  September  2010

Print ISBN: 9781420075403
eBook ISBN: 9781420075410

10.1201/9781420075410-5

#### Abstract

Materials at the nanoscale have been the subject of intensive study owing to their unusual electrical, magnetic, and optical properties [19]. The combination of new synthesis techniques offers unprecedented opportunities to tailor systems without resort to changing their chemical makeup. In particular, the physical properties of the material can be modified by confinement. If we confine a material in three, two, or one dimensions, we create a nanocrystal, a nanowire, or a nanofilm, respectively. Such confined systems often possess dramatically different properties than their macroscopic counterparts. As an example, consider the optical properties of the semiconductor cadmium selenide. By creating nanocrystals of different sizes, typically ∼5 nm in diameter, the entire optical spectrum can be spanned. Likewise, silicon can be changed from an optically inactive material at the macroscopic scale to an optically active nanocrystal.

#### 3.1  Introduction

Materials at the nanoscale have been the subject of intensive study owing to their unusual electrical, magnetic, and optical properties [19]. The combination of new synthesis techniques offers unprecedented opportunities to tailor systems without resort to changing their chemical makeup. In particular, the physical properties of the material can be modified by confinement. If we confine a material in three, two, or one dimensions, we create a nanocrystal, a nanowire, or a nanofilm, respectively. Such confined systems often possess dramatically different properties than their macroscopic counterparts. As an example, consider the optical properties of the semiconductor cadmium selenide. By creating nanocrystals of different sizes, typically ∼5 nm in diameter, the entire optical spectrum can be spanned. Likewise, silicon can be changed from an optically inactive material at the macroscopic scale to an optically active nanocrystal.

The physical confinement of a material changes its properties when the confinement length is comparable to the quantum length scale. This can easily be seen by considering the uncertainty principle. At best, the uncertainty of a particle’s momentum, Δp, and its position, Δx, must be such that ΔpΔxћ, where h is Planck’s constant divided by 2π. For a particle in a box of length scale, Δx, the kinetic energy of the particle will scale as ∼1/Δx 2 and rapidly increase for small values of Δx. If the electron confinement energy becomes comparable to the total energy of the electron, then its physical properties will clearly change, i.e., once the confining dimension approaches the delocalization length of an electron in a solid, “quantum confinement” occurs.

We can restate this notion in a similar context by estimating a confinement length, which can vary from solid to solid, and by considering the wave properties of the electron. For example, the de Broglie wavelength of an electron is given by λ = h/p. If λ is in the order of Δx, then we expect confinement and quantum effects to be important. Of course, this is essentially the same criterion from the uncertainty principle as λ ∼ Δx would give ΔpΔxh with p ∼ Δp. For a simple metal, we can estimate the momentum from a free electron gas: p = ћk f. k f is the wave vector such that the Fermi energy, $E f = ℏ 2 k f 2 / 2 m$

. This yields k f = (3π2 n)1/3 where n is the electron density. The value of k f for a typical simple metal such as Na is ∼1010 m-1. This would give a value of λ = h/p = 2π/k f or roughly ∼1 nm. We would expect quantum confinement to occur on a length scale of a nanometer for this example.

To observe the role of quantum confinement in real materials, we need to be able to construct materials routinely on the nanometer scale. This is the case for many materials. Nanocrystals of many materials can be made, although sometimes it is difficult to determine whether crystallinity is preserved. Such nanostruc-tures provide a unique opportunity to study the properties at nanometer scales and to reveal the underlying physics occurring at reduced dimensionality. From a practical point of view, nano-structures are promising building blocks in nanotechnologies, e.g., the smallest nominal length scale in a modern CPU chip is in the order of 100 nm or less. At length scales less than this, the “band structure” of a material may no longer appear to be quasi-continuous. Rather the electronic energy levels may be described by discrete quantum energy levels, which change with the size of the system.

An understanding of the physics of confinement is necessary to provide the fundamental science for the development of nano optical, magnetic, and electronic device applications. This understanding can best be obtained by utilizing ab initio approaches [10]. These approaches can provide valuable insights into nanoscale phenomena without empirical parameters or adjustments extrapolated from bulk properties [11,12]. As recently as 15 years ago, it was declared that ab initio approaches would not be useful to systems with more than a hundred atoms or so [11,12]. Of course, hardware advances have occurred since the mid-1990s, but more significant advances have occurred in the area of algorithms and new ideas. These ideas have allowed one to progress at a much faster rate than suggested by Moore’s law. We will review some of these advances and illustrate their application to nanocrystals. We will focus on two examples: a silicon nanocrystal and an iron nanocrystal. These two examples will illustrate the behavior of quantum confinement on the optical gap in a semiconductor and on the magnetic moment of a ferromagnetic metal.

#### 3.2  The Quantum Problem

The spatial and energetic distributions of electrons with the quantum theory of materials can be described by a solution of the Kohn–Sham eigenvalue equation [13], which can be justified using density functional theory [13,14]:

3.1 $[ − ℏ 2 ∇ 2 2 m + V ext ( r → ) + V H ( r → ) + V xc ( r → ) ] Ψ n ( r → ) = E n Ψ n ( r → )$

where m is the electron mass. The eigenvalues correspond to energy levels, En , and the eigenfunctions or wave functions are given by Ψ n ; a solution of the Kohn–Sham equation gives the energetic and spatial distribution of the electrons. The external potential, V ext, is a potential that does not depend on the electronic solution. The external potential can be taken as a linear superposition of atomic potentials corresponding to the Coulomb potential produced by the nuclear charge. In the case of an isolated hydrogen atom, V ext = –e2/r. The potential arising from the electron–electron interactions can be divided into two parts. One part represents the “classical” electrostatic terms and is called the “Hartree” or “Coulomb” potential:

3.2 $∇ 2 V H ( r → ) = − 4 π e ρ ( r → )$

where ρ is the electron charge density; it is obtained by summing up the square of the occupied eigenfunctions and gives the probability of finding an electron at the point $r →$

3.3 $ρ ( r → ) = e ∑ n , occup | Ψ n ( r → ) | 2$

The second part of the screening potential, the “exchange-correlation” part of the potential, V xc, is quantum mechanical in nature and effectively contains the physics of the Pauli exclusion principle. A common approximation for this part of the potential arises from the local density approximation, i.e., the potential depends only on the charge density at the point of interest, $V xc r → = V xc [ ρ ( r → ) ]$

.

Figure 3.1   Self-consistent field loop. The loop is repeated until the “input” and “output” charge densities are equal to within some specified tolerance.

In principle, the density functional theory is exact, provided one is given an exact functional for V xc. This is an outstanding research problem. It is commonly assumed that the functional extracted for a homogeneous electron gas [15] is “universal” and can be approximated by resort to the inhomogeneous gas problem. The procedure for generating a self-consistent field (SCF) potential is given in Figure 3.1. The SCF cycle is initiated with a potential constructed by a superposition of atomic densities for a nanostructure of interest. (Charge densities are easy to obtain for an atom. Under the assumption of a spherically symmetric atom, the Kohn-Sham equation becomes one dimensional and can be solved by doing a radial integration.) The atomic densities are used to solve a Poisson equation for the Hartree potential, and a density functional is used to obtain the exchange-correlation potential. A screening potential composed of the Hartree and exchange-correlation potentials is then added to the fixed external potential, after which the Kohn-Sham equation is solved. The resulting wave functions from this solution are then employed to construct a new potential and the cycle is repeated. In practice, the “output” and “input” potentials are mixed using a scheme that accounts for the history of the previous iterations [16,17].

This procedure is difficult because the eigenvalues can span a large range of energies and the corresponding eigenfunctions span disparate length scales. Consider a heavy element such as Pb. Electrons in the 1s state of Pb possess relativistic energies and are strongly localized around the nucleus. In contrast, the Pb 6s electrons are loosely bound and delocalized. Attempting to describe the energies and wave functions for these states is not trivial and cannot easily be accomplished using simple basis functions such as plane waves. Moreover, the tightly bound core electrons in atoms are not chemically active and can be removed from the Kohn-Sham equation without significant loss of accuracy by using the pseudopotential model of materials.

Figure 3.2   Pseudopotential model of a solid.

The pseudopotential model is quite general and reflects the physical content of the periodic table. In Figure 3.2, the pseudopotential model is illustrated for a crystal. In the pseudopo-tential model of a material, the electron states are decomposed into core states and valence states, e.g., in silicon the 1s 22s 22p 6 states represent the core states, and the 3s23p2 states represent the valence states. The pseudopotential represents the potential arising from a combination of core states and the nuclear charge: the so-called ion-core pseudopotential. The ion-core pseudopotential is assumed to be completely transferable from the atom to a cluster or to a nanostructure.

By replacing the external potential in the Kohn-Sham equation with an ion-core pseudopotential, we can avoid considering the core states altogether. The solution of the Kohn-Sham equation using pseudopotentials will yield only the valence states. The energy and length scales are then set by the valence states; it becomes no more difficult to solve for the electronic states of a heavy element such as Pb when compared to a light element such as C.

#### 3.2.1  Constructing Pseudopotentials from Density Functional Theory

Here we will focus on recipes for creating ion-core pseudopoten-tials within the density functional theory, although pseudopo-tentials can also be constructed from experimental data [18]. The construction of ion-core pseudopotentials has become an active area of electronic structure theory. Methods for constructing such potentials have centered on ab initio or “first-principles” pseudopo-tentials; i.e., the informational content on which the pseudopoten-tial is based does not involve any experimental input.

The first step in the construction process is to consider an electronic structure calculation for a free atom. For example, in the case of a silicon atom the Kohn-Sham equation [13] can be solved for the eigenvalues and wave functions. Knowing the valence wave functions, i.e., 3s 2 and 3p 2, states and corresponding eigenvalues, the pseudo wave functions can be constructed. Solving the Kohn-Sham problem for an atom is an easy numerical calculation as the atomic densities are assumed to possess spherical symmetry and the problem reduces to a one-dimensional radial integration. Once we know the solution for an “all-electron” potential, we can invert the Kohn-Sham equation and find the total pseudopotential. We can “unscreen” the total potential and extract the ion-core pseudopotential. This ion-core potential, which arises from tightly bound core electrons and the nuclear charge, is not expected to change from one environment to another. The issue of this “transferability” is one that must be addressed according to the system of interest. The immediate issue is how to define pseudo-wave functions that can be used to define the corresponding pseudopotential.

Suppose we insist that the pseudo-wave function be identical to the all-electron wave function outside of the core region. For example, let us consider the 3s state for a silicon atom. We want the pseudo-wave function to be identical to the all-electron state outside the core region:

3.4

where

• $ϕ 3 s p$ is a pseudo-wave function for the 3s state
• r c defines the core size

This assignment will guarantee that the pseudo-wave function will possess properties identical to the all-electron wave function, ψ3s , in the region away from the ion core.

For r < r c, we alter the all-electron wave function. We are free to do this as we do not expect the valence wave function within the core region to alter the chemical properties of the system. We choose to make the pseudo-wave function smooth and nodeless in the core region. This initiative will provide rapid convergence with simple basis functions. One other criterion is mandated. Namely, the integral of the pseudocharge density within the core should be equal to the integral of the all-electron charge density. Without this condition, the pseudo-wave function differs by a scaling factor from the all-electron wave function. Pseudopotentials constructed with this constraint are called “norm conserving” [19]. Since we expect the bonding in a solid to be highly dependent on the tails of the valence wave functions, it is imperative that the normalized pseudo-wave function be identical to the all-electron wave functions.

There are many ways of constructing “norm-conserving” pseudopotentials as within the core the pseudo-wave function is not unique. One of the most straightforward construction procedures is from Kerker [20] and was later extended by Troullier and Martins [21].

3.5 $ϕ l p ( r ) = { r l exp ( p ( r ) ) r ≤ r c ψ l ( r ) r > r c$

p(r) is taken to be a polynomial of the form

3.6 $p ( r ) = c 0 + ∑ n = 1 6 c 2 n r 2 n$

This form assures us that the pseudo-wave function is nodeless and by taking even powers there is no cusp associated with the pseudo-wave function. The parameters, c 2n , are fixed by the following criteria: (a) The all-electron and pseudo-wave functions have the same valence eigenvalue. (b) The pseudo-wave function is nodeless and be identical to the all-electron wave function for r > r c. (c) The pseudo-wave function must be continuous as well as the first four derivatives at r c. (d) The pseudopotential has zero curvature at the origin. This construction is easy to implement and extend to include other constraints. An example of an atomic pseudo-wave function for Si is given in Figure 3.3 where it is compared to an all-electron wave function. Unlike the 3s all-electron wave function, the pseudo-wave function is node-less. The pseudo-wave function is much easier to express as a Fourier transform or a combination of Gaussian orbitals than the all-electron wave function.

Once the pseudo-wave function is constructed, the Kohn-Sham equation can be inverted to arrive at the ion-core pseudopotential

3.7 $V ion , l p ( r ) = ℏ 2 ∇ 2 ϕ l p 2 m ϕ l p − E n , l − V H ( r ) − V xc [ ρ ( r ) ]$

The ion-core pseudopotential is well behaved as $ϕ l p$

has no nodes; however, the resulting ion-core pseudopotential is both state dependent and energy dependent. The energy dependence is usually weak. For example, the 4s state in silicon computed by the pseudopotential constructed from the 3s state is usually accurate.

Figure 3.3   An all-electron and a pseudo-wave function for the silicon 3s radial wave function.

Physically this happens because the 4s state is extended and experiences the potential in a region where the ion-core potential has assumed a simple −Zve 2/r behavior where Zv is the number of valence electrons. However, the state dependence through l is an issue, the difference between a potential generated via a 3s state and a 3p can be an issue. In particular, for first-row elements such as C or O, the nonlocality is quite large as there are no p states within the core region. For the first-row transition elements such as Fe or Cu, this is also an issue as again there are no d-states within the core. This state dependence complicates the use of pseudopotentials. The state dependence or the nonlocal character of the ion-core pseudopotential for an atom can be expressed as

3.8 $V ion p ( r → ) = ∑ l = 0 ∞ P l † V l , ion p ( r ) P l = P s † V s , ion p ( r ) P s + P p † V p , ion p ( r ) P p + P d † V d , ion p ( r ) P d + ⋯$

$P l$

is an operator that projects out the lth-component. The ion-core pseudopotential is often termed “semi-local” as the potential is radial local, but possesses an angular dependence that is not local.

An additional advantage of the norm-conserving potential concerns the logarithmic derivative of the pseudo-wave function [22]. An identity exists:

3.9

The energy derivative of the logarithmic derivative of the pseudo-wave function is fixed by the amount of charge within a radius, R. The radial derivative of the wave function, ϕ, is related to the scattering phase shift from elementary quantum mechanics. For a norm-conserving pseudopotential, the scattering phase shift at R = r c and at the eigenvalue of interest is identical to the all-electron case as Q all elect(r c) = Q pseudo(r c). The scattering properties of the pseudopotential and the all-electron potential have the same energy variation to first order when transferred to other systems.

There is some flexibility in constructing pseudopotentials; the pseudo-wave functions are not unique. This aspect of the pseudo-wave function was recognized early in its inception, i.e., there are a variety of ways to construct the wave functions in the core. The non-uniqueness of the pseudo-wave function and the pseudopotential can be exploited to optimize the convergence of the pseudopotentials for the basis of interest. Much effort has been made to construct “soft” pseudopotentials. By “soft,” one means a “rapidly” convergent calculation using plane waves as a basis. Typically, soft potentials are characterized by a “large” core size, i.e., a larger value for rc. However, as the core becomes larger, the “goodness” of the pseudo-wave function can be compromised as the transferability of the pseudopotential becomes more limited.

A schematic illustration of the difference between an “all-electron” potential and a pseudopotential is given in Figure 3.4.

Figure 3.4   Schematic all-electron potential and pseudopotential. Outside of the core radius, r c, the potentials and wave functions are identical.

#### 3.2.2  Algorithms for Solving the Kohn-Sham Equation

The Kohn-Sham equation as cast in Equation 3.1 can be solved using a variety of techniques. Often the wave functions can be expanded in a basis such as plane waves or Gaussians and the resulting secular equations can be solved using standard diago-nalization packages such as those found in the widely used code: vasp [23]. vasp is a particularly robust code, with a wide following, but it was not constructed with a parallel computing environment in mind.

Here we focus on a different approach that is particularly targeted at highly parallel computing platforms. We solve the Kohn-Sham equation without resort to an explicit basis [2432]. We solve for the wave functions on a uniform grid within a fixed domain. The wave functions outside of the domain are required to vanish for confined systems or we can assume periodic boundary conditions for systems with trans-lational symmetry [3335]. In contrast to methods employing an explicit basis, such boundary conditions are easily incorporated. In particular, real space methods do not require the use of supercells for localized systems. As such, charged systems can easily be examined without considering any electrostatic divergences. The problem is typically solved on a uniform grid as indicated in Figure 3.5.

Within a “real space” approach, one can solve the eigenvalue problem using a finite element or finite difference approach [2432]. We use a higher order finite difference approach owing to its simplicity in implementation. The Laplacian operator can be expressed using

Figure 3.5   Real space geometry for a confined system. The grid is uniform and the wave function is taken to vanish outside the domain of interest.

3.10 $( ∂ 2 ψ ∂ x 2 ) x 0 ≈ ∑ n = − N N A n ψ ( x 0 + n h , y , z ) ,$

where

• h is the grid spacing
• N is the number of nearest grid points
• An are the coefficients for evaluating the required derivatives [36]

The error scales as O(h 2N +2). Typically, N ≈ 6 − 8 is used and there is a trade off between using a higher value for N and a coarser grid, or a smaller value for N and a finer grid. Because pseudopotentials are used in this implementation, the wave function structure is “smooth” and a higher finite difference expression for the kinetic energy converges quickly for a fine grid.

In real space, we can easily incorporate the nonlocal nature of the ion-core pseudopotential [24]. The Kleinman-Bylander form [37] can be expressed in real space as

3.11

where $u l m p$

are the reference atomic pseudo-wave functions. The nonlocal nature of the pseudopotential is apparent from the definition of G lm , the value of these coefficients are dependent on the pseudo-wave function, $ϕ l p$ , acted on by the operator ΔVl . This is very similar in spirit to the pseudopotential defined by Phillips and Kleinman [38].

Once the secular equation is created, the eigenvalue problem can be solved using iterative methods [32,39,40]. Typically, a method such as a preconditioned Davidson method can be used [32]. This is a robust and efficient method, which never requires one to store explicitly the Hamiltonian matrix.

Recent work avoids an explicit diagonalization and instead improves the wave functions by filtering approximate wave functions using a damped Chebyshev polynomial filtered sub-space iteration [32]. In this approach, only the initial iteration necessitates solving an eigenvalue problem, which can be handled by means of any efficient eigensolver. This step is used to provide a good initial subspace (or good initial approximation to the wave functions). Because the subspace dimension is slightly larger than the number of wanted eigenvalues, the method does not utilize as much memory as standard restarted eigensolvers such as ARPACK and TRLan (Thick—Restart, Lanczos) [41,42]. Moreover, the cost of orthogonalization is much reduced as the filtering approach only requires a subspace with dimension slightly larger than the number of occupied states and orthogo-nalization is performed only once per SCF iteration. In contrast, standard eigensolvers using restart usually require a subspace at least twice as large and the orthogonalization and other costs related to updating the eigenvectors are much higher.

The essential idea of the filtering method is to start with an approximate initial eigenbasis, {ψ n }, corresponding to occupied states of the initial Hamiltonian, and then to improve adap-tively the subspace by polynomial filtering. That is, at a given self-consistent step, a polynomial filter, Pm (t), of order m is constructed for the current Hamiltonian $H$

. As the eigen-basis is updated, the polynomial will be different at each SCF step since $H$ will change. The goal of the filter is to make the subspace spanned by ${ ψ ̂ n } = P m ( H ) { ψ n }$ approximate the eigen subspace corresponding to the occupied states of $H$ . There is no need to make the new subspace, ${ ψ ̂ n }$ , approximate the wanted eigen subspace of $H$ to high accuracy at intermediate steps. Instead, the filtering is designed so that the new subspace obtained at each self-consistent iteration step will progressively approximate the wanted eigen space of the final Hamiltonian when self-consistency is reached.

This can be efficiently achieved by exploiting the Chebyshev polynomials, Cm , for the polynomials Pm. In principle, any set of polynomials would work where the value of the polynomial is large over the interval of interest and damped elsewhere. Specifically, we wish to exploit the fast growth property of Chebyshev polynomials outside of the [-1, 1] interval. All that is required to obtain a good filter at a given SCF step, is to provide a lower bound and an upper bound of an interval of the spectrum of the current Hamiltonian $H$

. The lower bound can be readily obtained from the Ritz values computed from the previous step, and the upper bound can be inexpensively obtained by a very small number of (e.g., 4 or 5) Lanczos steps [32]. The main cost of the filtering at each iteration is in performing the products of the polynomial of the Hamiltonian by the basis vectors; this operation can be simplified by utilizing recursion relations. To construct a “damped” Chebyshev polynomial on the interval [a, b] to the interval [-1, 1], one can use an affine mapping such that

3.12 $l ( t ) = t − ( a + b ) / 2 ( b − a ) / 2 .$

Figure 3.6   A damped Chebyshev polynomial, C 6. Th e shaded area corresponds to eigenvalue spectrum regime that will be enhanced by the filtering operation (see text).

The interval is chosen to encompass the energy interval containing the eigen space to be filtered. The filtering operation can then be expressed as

3.13 ${ ψ ̂ n } = C m ( l ( H ) ) { ψ n } .$

This computation is accomplished by exploiting the convenient three-term recursion property of Chebyshev polynomials:

3.14

An example of a damped Chebyshev polynomial as defined by Equations 3.12 and 3.14 is given in Figure 3.6 where we have taken the lower bound as a = 0.2 and the upper bound as b = 2. In this schematic example, the filtering would enhance the eigenvalue components in the shaded region.

The filtering procedure for the self-consistent cycle is illustrated in Figure 3.7. Unlike traditional methods, the cycle only requires one explicit diagonalization step. Instead of repeating the diagonalization step within the self-consistent loop, a filtering operation is used to create a new basis in which the desired eigen subspace is enhanced. After the new basis, ${ ψ ̂ n }$

, is formed, the basis is orthogonalized. The orthogonalization step scales as the cube of the number of occupied states and as such this method is not an “order-n” method. However, the prefactor is sufficiently small that the method is much faster than previous implementations of real space methods [32]. The cycle is repeated until the “input” and “output” density is unchanged within some specified tolerance, e.g., the eigenvalues must not change by ∼0.001 eV, or the like.

In Table 3.1, we compare the timings using the Chebyshev filtering method along with explicit diagonalization solvers using the TRLan and the ARPACK. These timings are for a modest-sized nanocrystal: Si525H276. The Hamiltonian size is 292,584 × 292,584 and 1194 eigenvalues were determined. The numerical runs were performed on the SGI Altix 3700 cluster at the Minnesota Supercomputing Institute. The CPU type is a 1.3 GHz Intel Madison processor. Although the number of matrix-vector products and SCF iterations is similar, the total time with filtering is over an order of magnitude faster compared to ARPACK and a factor of better than four versus the TRLan. The scaling of the algorithm with the number of processors is shown in Figure 3.8. Such improved timings are not limited to this particular example. Our focus here is on silicon nanocrystals, our method is not limited to semiconducting or insulating systems, we have also used this method to examine metallic systems such as liquid lead [43] and magnetic systems such as iron clusters [8].

Figure 3.7   Self-consistent cycle using damped Chebyshev filtering. Atomic units (e = ћ = m) are used here.

### Table 3.1   Comparison of Computational Timings for Various Methods for a Nanocrystal: Si525H276

Method

SCF Its

CPU(s)

Filtering

11

5,947

ARPACK

10

62,026

TRLan

10

26,853

Note: While the number of SCF iterations is comparable for all three methods, the total time with filtering methods can be dramatically reduced.

Figure 3.8   Examples of performance scaling for the real space pseudopotential code.

#### 3.3.1.1  Intrinsic Properties

If we wish to examine the intrinsic properties of nanocrystals, we need to deal with the crystal surface, which at the nanoscale becomes increasingly important. Experimentally, this issue is often handled by passivating the intrinsic surface with surfactants or hydrogenating the surface. We chose to use the latter approach by capping surface dangling bonds with hydrogen atoms [44]. The largest nanocrystal we examined contained over ten thousand atoms: Si9041H1860, which is approximately 7 nm in diameter [32,45]. A ball and stick model of a typical nanocrystal is shown in Figure 3.9.

Figure 3.9   The ball and stick model of a hydrogenated silicon quantum dot. The interior consists of a diamond fragment. The surface of the fragment is capped with hydrogen atoms.

A solution of the Kohn-Sham equations yields the distribution of eigenvalues. For sufficiently large nanocrystals, one expects the distribution to approach that of crystalline silicon, i.e., the distribution of states should approach that of the crystalline density of states. The range of sizes for these nanocrystals allows us to make comparisons with the bulk crystal. This comparison is made in Figure 3.10. The “density of states” (DOS) for the eigenvalue spectrum for the nanocrystal shares similar structure as in the bulk crystal. For example, the sharp peak at about 4 eV below the top of the valence band arises from an M 1 critical point along the [110] direction in the bulk [18]. This feature is clearly present in the nanocrystal eigenvalue spectra. The only notable differences between the crystal and the nanocrystal occur because of the presence of the Si-H bonds and the lack of completely evolved bands.

We can also examine the evolution of the ionization potentials (I) and the electron affinities (A) for the nanocrystal.

3.15 $I = E ( N − 1 ) − E ( N ) A = E ( N ) − E ( N + 1 )$

In principle, these are ground state properties and, if the correct functional were known, these quantities would be accurately predicted by density functional theory. For atoms and molecules, one typically extracts accurate values for I and A, e.g., for the first-row atoms, the error is typically less than ∼5%. The difference between the ionization potential and the electron affinity can be associated with the quasi-particle gap: E qp = I - A. If the exciton (electron-hole) interaction is small, this gap can be compared to the optical gap. However, for silicon nanocrys-tals the exciton energy is believed to be on the order of ∼1 eV for nanocrystals of less than ∼1 nm [46].

Figure 3.10   The eigenvalue density of states for the Si9041H1860 nanocrystal (top panel) and the electronic density of states for crystalline silicon (bottom panel). Th e highest occupied state in both panels is taken to be the zero energy reference.

Figure 3.11   The evolution of the ionization potential (IP) and electron affinity (EA) with quantum dot size. Also shown are the eigenvalue levels for the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO).

We can examine the scaling of the ionization potential and affinity by assuming a simple scaling and fitting to the calculated values (shown in Figure 3.11):

3.16 $I ( D ) = I ∞ + A / D α A ( D ) = A ∞ + B / D β$

where D is the dot diameter. A fit of these quantities results in I = 4.5 eV, A = 3.9 eV, α = 1.1 and β = 1.08. The fit gives a quasi-particle gap of E qp(D → ∞) = I - A = 0.6 eV in the limit of an infinitely large dot. This value is in good agreement with the gap found for crystalline silicon using the local density approximation [47,48]. The gap is not in good agreement with experiment owing to the failure of the local density approximation to describe band gaps of bulk semiconductors in general.

A key aspect of our study is that we can examine the scaling of the ionization potential and electron affinity for nanocrystals ranging from silane (SiH4) to systems containing thousands of atoms. We not only verify the limiting value of the quasi-particle gap, we can ascertain how this limit is reached, i.e., how the ionization potential and electron affinity scale with the size of the dot and what the relationship is between these quantities and the highest occupied and lowest empty energy levels. Since values of I and A from LDA are reasonably accurate for atoms and molecules, one can ask how the size of the nanocrystal affects the accuracy of LDA for predicting I and A. Unfortunately, experimental values for the ionization potentials and electron affinities are not known for hydrogenated silicon clusters and nanocrystals, a notable exception being silane (where the electron affinity is negative) [49]. However, there is some theoretical evidence that the issue involves errors in the ionization energies as opposed to the electron affinities. In Figure 3.12, we display the ionization potentials and electron affinities from the GW approximation [4951] and compare to calculations using the LDA approximation. Within GW, the electron self-energy is obtained by calculating the lowest order diagram in the dynamically screened Coulomb interaction [49]. It is expected that the values for (IP,EA) are more accurately replicated by GW than the LDA work. The LDA and GW values for the EA are nearly identical. However, GW yields slightly larger IP values than LDA. As a function of size, the difference between the IP and GW values is small for the smallest nanocrystals and appears to saturate for the larger ones.

Figure 3.12   First ionization potential and electron affinity of passivated silicon clusters, calculated within the GW approximation [49] solid lines and ΔSCF dotted lines. Experimental data from Ref. [104] SCF results include spin-polarization eff ects.

Given the energy as a function of position, we can also examine structural and vibrational properties of nanocrystals. Vibrational properties are easier to describe as they converge more rapidly to the bulk values than do electronic states; however, because they involve small changes in energy with position, the wave functions need to be better converged.

Vibrational modes in nanocrystals can play an important role within the context of the photovoltaic applications. In particular, vibrational properties are directly related to the phonon-assisted optical transitions. This is an important consideration for silicon: in the bulk limit, the lowest optical transitions are indirect and must be phonon assisted. It is expected that the vibrational dynamics in the presence of a free carrier will be different from the case of an intrinsic nanocrystal. Strains induced by a free uncompensated carrier might break the symmetry of the vibrational modes, thereby increasing their participation probability in the optical transitions.

Owing to the localized nature of nanocrystals, it is feasible to predict vibrational modes calculations by the direct force-constant method [52]. The dynamical matrix of the system is constructed by displacing all atoms one by one from their equilibrium positions along the Cartesian directions and finding the forces induced on the other atoms of the nanocrystal. We determine the forces using the Hellmann-Feynman theorem in real space [33] and employ a symmetrized form of the dynamical matrix expression [53]. The elements of the dynamical matrix, $D i j αβ$

, are given by

3.17 $D i j αβ = − 1 2 [ F i α ( { R } + d j β ) − F i α ( { R } − d j β ) 2 d j β ] 000000000 + F j β ( { R } + d i α ) − F j β ( { R } − d i α ) 2 d i α$

where

• $F j α$ is the force on atom α in the direction i
• ${ R } + d j β$ is the atomic configuration where only the atom β is displaced along j from its equilibrium position

The value of displacement was chosen to be 0.015 Å. The equilibrium structure was relaxed so that the maximum residual forces were less than 10-3 eV/Å. For this accuracy, the grid spacing h is reduced to 0.4 a.u. (1 a.u. = 0.529 Å) as compared to a value of about 0.7 a.u., typically used for electronic properties. The Chebyshev filtering algorithm is especially well suited for this procedure as the initial diagonalization need not be repeated when the geometry changes are small.

The vibrational modes frequencies and corresponding eigenvectors can be obtained from the dynamical equation

3.18 $∑ β , k [ ω 2 δ αβ δ i k − D i k αβ M α M β ] A k β = 0$

where M a is the mass of atom labeled by α.

We apply this procedure to a Si nanocrystal: Si29H36. The surface of this nanocrystal is passivated by hydrogen atoms in order to electronically passivate the nanocrystal, i.e., remove any dangling bond states from the gap [54]. (Si clusters without hydrogen passivation has been examined previously, assuming notable surface reconstructions [55]). The choice of the number of silicon atoms is dictated so that the outermost Si atoms are passivated by no more than two hydrogen atoms. Our nanocrystal has a bulk-like geometry with the central atom having the Td symmetry. The bonding distance to the four equivalent nearest neighbors is 2.31 Å.

We also consider cation nanocrystals by removing an electron from the system: (Si29H36)+. A relaxation of this nanocrystal leads to a distortion owing to a Jahn-Teller effect structure, i.e., there is a partial filling of the highest occupied state. This symmetry breaking leads to the central atom of the charged nanocrystal being bonded to four atoms with different bond lengths. Two bond lengths become slightly shorter (2.30 Å) and two others lengthen (2.35 Å). This distortion propagates throughout the nanocrystal.

Figure 3.13   Vibrational density of states of Si29H36.

Vibrational densities of states of the two nanocrystals are shown in Figure 3.13. The DOS for the nanocrystal can be divided into four distinct regions. (Vibration spectra of various silicon clusters were studied in a number of earlier papers [56,57] using model potential calculations [58]). The lowest energy vibrations involve only the Si atoms. We will assign them to “core modes.” In the core modes all Si atoms are involved, while the passivat-ing hydrogens remain static. The lowest frequency modes in this region (below ∼250 cm-1) correspond to modes with no bond stretching, i.e., bond bending modes dominate. For modes above ∼250 cm-1 stretching dynamics become more important. In the Si-Si “stretching-bending” part of the spectrum, the important part is the highest peak right below 500 cm-1. The nature of this nanocrystal peak can be ascribed to modes that are characteristic transverse optical (TO) mode of crystalline silicon. The TO mode is extensively discussed in the literature as it is Raman active. The mode is sensitive to the nanocrystal size [59,60]. The next region of the vibration modes, just above 500 cm-1, is related to the surface vibrations of Si and H atoms. The fourth distinct region is located around 1000 cm-1 and is related to the H atoms scissor-bending modes; the Si atoms do not move in this case. The highest energy vibrations, above 2000 cm-1, are the Si-H stretching modes.

The cation nanocrystal, (Si29H36)+, has a similar vibrational spectrum. However, in the charged system we see that hydrogen actively participates in the vibrations well below 500 cm-1 (Figure 3.14). The peak at 500 cm-1 is red shifted in the case of the changed system, which is associated with a much larger contribution from the surface atoms. The vibrational DOS spectrum of the hydrogen, shown in Figure 3.14, demonstrates that hydrogen atoms vibrate at frequencies all the way down to 400 cm-1. This is rather unusual, and occurs owing to an enhanced charge transfer on the surface associated with the hole charge. As noted in Ref. [55], the distribution of the hole density affects both lattice relaxations and the corresponding vibrational spectra. In order to quantify the location of the hole charge in the (Si29H36)+ nanocrystal, we calculated the “electron localization function” (ELF) for both neutral and charged structures [61]. Figure 3.15 presents two orientations of the ELF plots. On the left-hand side plots (a) and (c), the neutral nanocrystals is shown; on the right-hand side, plots (b) and (d), the charged nanocrystal is shown. Flat cuts of the ELF at the surface of the charged nanocrystal show where the hole is located.

Figure 3.14   Vibrational density of states for the hydrogen contributions in Si29H36.

Figure 3.15   Electron localization function: (a,c) the neutral Si29H36 nanocrystal shown in two orientations; (b,d) the Si29H36 1+ cation, also in two orientations.

#### 3.3.1.2  Extrinsic Properties

Doping a small percentage of foreign atoms in bulk semiconductors can profoundly change their electronic properties and makes possible the creation of modern electronic devices [62]. Phosphorus-doped crystalline Si introduces defect energy states close to the conduction band. For such shallow donors, electrons can be easily thermally excited, greatly enhancing the conductivity of the original pure semiconductor by orders of magnitude at room temperature.

The evolution of the semiconductor industry requires continued miniaturization. The industry is maintaining exponential gains in the performance of electronic circuits by designing devices ever smaller than the previous generation. This device miniaturization is now approaching the nanometer-scale. As a consequence, it is of the utmost importance to understand how doping operates at this length scale as quantum confinement is expected to alter the electronic properties of doped Si nanocrys-tals [63]. Also, doped Si nanowires have been synthesized and it has been demonstrated experimentally that they can be used as interconnects in electronic circuits or building blocks for semiconductor nanodevices [64,65]. Important questions arise as to whether the defect energy levels are shallow or not, e.g., at what length scale will device construction based on macroscopic laws be altered by quantum confinement?

Phosphorus-doped silicon nanocrystals represent the prototypical system for studying impurities in quantum dots. Recent experiments, designed to study this system, have utilized photo-luminescence [66,67] and electron spin resonance measurements [6870]. Electron spin resonance experiments probe the defect energy levels through hyperfine interaction. Hyperfine splitting (HFS) arises from the interaction between the electron spin of the defect level and the spin of the nucleus, which is directly related to the probability of finding a dopant electron localized on the impurity site [71]. A HFS much higher than the bulk value of 42 G has been observed for P-doped Si nanocrystals with radii of 10 nm [68]. A size dependence of the HFS of P atoms was also observed in Si nanocrystals [69,70].

Unfortunately, theoretical studies of shallow impurities in quantum dots are computationally challenging. Owing to the large number of atoms and to the low symmetry of the system involved, most total energy calculations have been limited to studying nanocrystals that are much smaller than the size synthesized in experiment [7275]. While empirical studies have been performed for impurities in large quantum dots, they often utilize parameters that are ad hoc extrapolations of bulk-like values [7678].

The same methods used to examine intrinsic silicon can be used for extrinsic silicon [45]. It is fairly routine to examine P-doped Si nanocrystals up to a diameter of 6 nm, which spans the entire range of experimental measurements [70]. The HFS size dependence is a consequence of strong quantum confinement, which also leads to the higher binding energy of the dopant electron. Hence, P is not a shallow donor in Si if the nanocrystals are less than 20 nm in diameter. In addition, we find that there is “critical” nanocrystal size below which the P donor is not stable against migration to the surface.

As for intrinsic properties, our calculations are based on density functional theory in the local density approximation [13,14]. However, the grid spacing is chosen to be 0.55 a.u., as a finer grid must be employed to converge the system owing to the presence of the P dopant [10,45]. The geometry of Si nanocrystal is taken to be bulk-like and roughly spherical in shape, in accord with experimental observation [79]. Again, the dangling bonds on the surface of the nanocrystal are passivated by H. The experimentally synthesized Si nanocrystals are usually embedded in an amorphous silicon dioxide matrix. The Si/SiO2 interface is in general not the same as H passivation. Nevertheless, both serve the role of satisfying the dangling bonds on the surface. The Si nanocrystals are then doped with one P atom, which substitutes a Si atom in the nanocrystal.

In Figure 3.16a, the defect state charge density along the [100] direction is illustrated. As the nanocrystal size increases, the defect wave function becomes more delocalized. This role of quantum confinement is observed in both experiments [70] and theoretical calculations [72,80]. The maxima of the charge density at around 0.2 nm correspond to the bond length between the P at the origin and its first Si neighbors. We can smooth out these atomic details by spherically averaging the defect wave function as shown in Figure 3.16b. We find that the defect wave function decays exponentially from the origin. This corresponds very well to the conventional understanding of defects in semiconductors: the defect ion and the defect electron form a hydrogen-like system with the wave function described as $ψ ∼ exp ( − r / a B eff )$

.

Figure 3.16   (a) Charge density for the dopant electron along the [100] direction for three P-doped Si nanocrystals with different radius. x is the coordinate along that direction. (b) The corresponding spherically averaged charge densities. (c) The effective Bohr radius a B eff corresponding to the dopant electron as a function of nanocrystal radius. (d) The calculated HFS of P-doped Si nanocrystals as a function of nanocrystal radius together with experimental data (▴) from Ref. [70]. Theoretical values for both the unrelaxed bulk geometries (•) and the fully relaxed structures (♦) are shown.

From the decay of the defect wave function, we can obtain an effective Bohr radius $a B eff$

and its dependence on nanocrystal size as plotted in Figure 3.16c. We find that the effective Bohr radius varies nearly linearly with nanocrystal radius R up to 3 nm where R is approximately five times $a B eff$ . Nonlinearity can be observed for very small Si nanocrystals, and $a B eff$ appears to converge to ∼0.2 nm as R → 0. The limit trends to the size of a phosphorus atom while retaining its sp 3 hybridization. In the bulk limit, $a B eff$ is ∼2.3 nm assuming a dielectric constant of 11.4 and effective electron mass of 0.26 m e. The dependence on R should trend to this bulk limit when the diameter of the nano-crystal is sufficiently large.

From the defect wave functions, we can evaluate the isotro-pic HFS as well [81,82]. Our calculated HFS is plotted in Figure 3.16d for a P atom located at the center of Si nanocrystal. There is very good agreement between experimental data [70] and our theoretical calculations. For nanocrystals with radius between 2-3 nm, the HFS is around 100 G, which considerably exceeds the bulk value of 42 G. Moreover, the HFS continues to increase as the nanocrystal size decreases. This is a consequence of the quantum confinement as illustrated in Figure 3.16a and b. The defect wave function becomes more localized at the P site as the radius decreases, leading to higher amplitude of the wave function at the P core. Only Si nanocrystals containing less than 1500 Si atoms are structurally optimized to a local energy minimum. For Si nanocrystals with radius larger than 1.5 nm, the effect of relaxation diminishes and the difference between unrelaxed and relaxed results is within ∼10%.

The experimental methodology used to measure the HFS ensures that only nanocrystals with exactly one impurity atom are probed. However, experiment has no control over the location of the impurity within the nanocrystal. The spatial distribution of impurities cannot be inferred from experimental data alone. Therefore, we considered the energetics of the doped nanocrystal by varying the P position along the [100] direction as illustrated in Figure 3.17. We avoid substituting the Si atoms on the surface of the nanocrystal by P. Our results for five of the small nanocrystals after relaxation to local energy minimum are shown in Figure 3.17. For Si nanocrystals with a diameter smaller than -2 nm, P tends to substitute Si near the surface. Otherwise, there is a bistable behavior in which both the center and the surface of the nanocrystal are energetically stable positions. This suggests that a “critical size” exists for nanocrystals. Below this size, P atoms will always be energetically expelled toward the surface.

We also calculated the binding energy E B of the defect electron as the P position changes. The binding energy is a measure of how strongly the defect electron interacts with the P atom, and is calculated by the energy required to ionize a P-doped Si nanocrystal by removing an electron (Jd) minus the energy gained by adding the electron to a pure Si nanocrystal (Ap). Figure 3.18a illustrates the typical situation for Si nanocrystals larger than 2 nm in diameter: the binding energy tends to decrease as the P moves toward the surface. This decrease occurs because the defect wave function becomes more distorted and less localized around P, leading to a loss in the Coulomb energy between the P ion and the defect electron. The change in the defect wave function explains why the center of the nanocrystal is energetically favorable. However, since the doped nanocrystal can relieve its stress by expelling the P atom toward the surface where there is more room for relaxation, positions close to the surface are always locally stable as depicted in Figure 3.17. The binding energy and P-induced stress compete with each other in determining the defect position within the Si nanocrystal.

Figure 3.17   Difference in energy and HFS as the P atom moves away from the center of the Si nanocrystal. Th e energies are with respect to the energy of the Si nanocrystal with P at the center. The x-axis measures the distance of P atom away from the origin in the unit of Si bond length as illustrated in the perspective view of a Si nanocrystal.

For Si nanocrystals less than 2 nm in diameter, the binding energy is higher close to the surface of the Si nanocrystal as shown in Figure 3.18b and c. A comparison of the binding energy between relaxed and unrelaxed structures suggests that the reversal in the trend is caused by relaxation. From Figure 3.18e, the P atom is found to relax toward the center of the nano-crystal, leading to a more localized defect wave function as in Figure 3.18d with better confinement inside the nanocrystal, and therefore higher binding energy. The relaxation of P atom toward the center of the nanocrystal causes an expansion of the Si nano-crystal in the perpendicular direction creating strain throughout the nanocrystal. This trade off between the binding energy and stress is only feasible for small Si nanocrystals as depicted in Figure 3.18f. As the diameter of the Si nanocrystal increases, the relaxation trends to a bulk-like geometry. This interplay between the binding energy and stress for small Si nanocrystals stabilizes the P atom close to the surface of the nanocrystal.

The HFS evaluated for different P positions inside the Si nanocrystal is also shown in Figure 3.17. There is a general trend for the HFS to drop drastically as the P atom approaches the surface of the nanocrystal. Smaller variations in the HFS are found for P positions away from the surface. However, for a large Si nanocrys-tal, the HFS varies within ∼10% of the value with P at the center position. From studies of the energetics and comparison to the experimental results in Figure 3.16d, it is possible for the synthesized Si nanocrystals to have P located close to the center of the nanocrystals.

An analysis of the defect wave function in Figure 3.16 can be fitted to an effective mass model [45]. Motivated by the defect wave functions having an approximate form of exp(-r/a B), we consider a hydrogen-like atom in a “dielectric box.” The potential that the electron experiences in atomic units is

Figure 3.18   (a-c) Changes in binding energy as the P atom moves away from the center of three Si nanocrystals with different sizes. The same x-axis is used as in Figure 3.17. (d) The charge density of the dopant electron along the [100] direction (x direction) for the Si86H76P cluster with the P atom located close to the surface as illustrated in (e). (f) A plot of the angle subtended by the P atom located close to the surface for five different Si nanocrystals. Only the number of Si atoms is used to label the x-axis for clarity. The solid line represents the results after relaxation, while the dashed line corresponds to unrelaxed bulk-like geometries in the figure.

3.19

where

• R is the radius of the well
• V0 the well depth

The dielectric constant ε(R) depends on nanocrystal size [45] and is assumed to follow Penn’s model ε(R) = 1 + ((11.4 - 1) / (1+(α/R) n )) [83]. The dielectric constant converges to 11.4 in the bulk limit. a and n will be used as fitting parameters. An approximate solution is obtained to the Schrodinger equation for the defect electron with an effective mass m* = 0.26 me under this potential by using a trial wave function $ψ = π / a B 3 exp ( − r / a B )$

. By applying the variational principle, the energy can be minimized with respect to the effective Bohr radius a B. Hence, a B can be found as a function of the well depth V 0 and radius R. Alternatively, V 0 can be inferred since we know a B from our first-principles results. In fact, the energy calculated from this model corresponds to the binding energy E B of the dopant electron. Therefore, we can fit the α and n parameters in Penn’s model by calculating the binding energy. Our calculated size dependence of the binding energy (E B = I d - A p as defined above) is plotted in Figure 3.19a and b. A detailed explanation of the trend can be found in Refs. [45,72].

Figure 3.19   (a) The ionization energy of P-doped Si nanocrystal (▾) and the electron affinity of pure Si nanocrystal (▴) plotted as a function of nanocrystal radius. (b) The calculated binding energy E B using eff ective mass theory plotted together with our ab initio results. From effective mass theory, the binding energy has contributions from the kinetic energy (KE), the Coulomb interaction between the dopant ion and its electron (PECoulomb), and the potential energy (PEQW) due to the quantum well. (c) Th e potential well depth V 0 as a function of nanocrystal radius from effective mass theory based on a hydrogen atom in a box.

Figure 3.19 illustrates the results from the fitting to the effective mass model. By using α = 5.4 nm and n = 1.6, we find that our model can reproduce almost exactly the binding energies from the first-principles calculations. The binding energy E B scales as R -11 where R the nanocrystal radius. For nanocrystals up to 6 nm in diameter, the binding energy is significantly larger than k B T at room temperature. An extrapolation of our results shows that a nanocrystal diameter of at least 20 nm is needed in order for P to be a shallow donor. Interestingly, the depth of the potential well V 0 depends on the nanocrystal radius R rather than being infinite or a constant. The finite V 0 explains why the dependence of HFS on R scales with an exponent smaller than three, which is a consequence of an infinitely deep quantum well [84,85]. The well represents the effect of quantum confinement on the wave function, which for very large nanocrystals diminishes, corresponding to a vanishing quantum well in the bulk limit. As the nanocrystal size decreases, the quantum well becomes deeper such that it can confine the defect electron more effectively as its kinetic energy increases.

#### 3.3.2  Iron Nanocrystals

The existence of spontaneous magnetization in metallic systems is an intriguing problem because of the extensive technological applications of magnetic phenomena and an incomplete theory of its fundamental mechanisms. In this scenario, clusters of metallic atoms serve as a bridge between the atomic limit and the bulk, and can form a basis for understanding the emergence of magnetization as a function of size. Several phenomena such as ferromagnetism, metallic behavior, and ferroelectricity have been intensely explored in bulk metals, but the way they manifest themselves in clusters is an open topic of debate. At the atomic level, ferromagnetism is associated with partially filled 3d orbitals. In solids, ferromagnetism may be understood in terms of the itinerant electron model [86], which assumes a partial delocalization of the 3d orbitals. In clusters of iron atoms, delo-calization is weaker owing to the presence of a surface, whose shape affects the magnetic properties of the cluster. Because of their small size, iron clusters containing a few tens to hundreds of atoms are superparamagnetic: The entire cluster serves as a single magnetic domain, with no internal grain boundaries [87]. Consequently, these clusters have strong magnetic moments, but exhibit no hysteresis.

The magnetic moment of nano-sized clusters has been measured as a function of temperature and size [8890], and several aspects of the experiment have not been fully clarified, despite the intense work on the subject [9196]. One intriguing experimental observation is that the specific heat of such clusters is lower than the Dulong-Petit value, which may be due to a magnetic phase transition [89]. In addition, the magnetic moment per atom does not decay monotonically as a function of the number of atoms and for fixed temperature. Possible explanations for this behavior are structural phase transitions, a strong dependence of magnetization with the shape of the cluster, or coupling with vibrational modes [89]. One difficulty is that the structure of such clusters is not well known. First-principles and model calculations have shown that clusters with up to 10 or 20 atoms assume a variety of exotic shapes in their lowest-energy configuration [97,98]. For larger clusters, there is indication for a stable body-centered cubic (BCC) structure, which is identical to ferromagnetic bulk iron [91].

The evolution of magnetic moment as a function of cluster size has attracted considerable attention from researchers in the field [8898]. A key question to be resolved is: What drives the suppression of magnetic moment as clusters grow in size? In the iron atom, the permanent magnetic moment arises from exchange splitting: the 3d orbitals (majority spin) are low in energy and completely occupied with five electrons, while the 3d orbitals (minority spin) are partially occupied with one electron, resulting in a magnetic moment of 4 μB, μB being the Bohr magneton. When atoms are assembled in a crystal, atomic orbitals hybridize and form energy bands: 4s orbitals create a wide band that remains partially filled, in contrast with the completely filled 4s orbital in the atom; while the 3d and 3d orbitals create narrower bands. Orbital hybridization together with the different bandwidths of the various 3d and 4s bands result in weaker magnetization, equivalent to 2.2 μB/atom in bulk iron.

In atomic clusters, orbital hybridization is not as strong because atoms on the surface of the cluster have fewer neighbors. The strength of hybridization can be quantified by the effective coordination number. A theoretical analysis of magnetization in clusters and thin slabs indicates that the dependence of the magnetic moment with the effective coordination number is approximately linear [93,95,96]. But the suppression of magnetic moment from orbital hybridization is not isotropic [99]. If we consider a layer of atoms for instance, the 3d orbitals oriented in the plane of atoms will hybridize more effectively than orbitals oriented normal to the plane. As a consequence, clusters with faceted surfaces are expected to have magnetic properties different from clusters with irregular surfaces, even if they have the same effective coordination number [99]. This effect is likely responsible for a nonmonotonic suppression of magnetic moment as a function of cluster size. In order to analyze the role of surface faceting more deeply, we have performed first-principles calculations of the magnetic moment of iron clusters with various geometries and with sizes ranging from 20 to 400 atoms.

The Kohn-Sham equation can be applied to this problem using a spin-density functional. We used the generalized gradient approximation (GGA) [100] and the computational details are as outlined elsewhere [8,24,31,101]. Obtaining an accurate description of the electronic and magnetic structures of iron clusters is more difficult than for simple metal clusters. Of course, the existence of a magnetic moment means an additional degree of freedom enters the problem. In principle, we could consider non-collinear magnetism and associate a magnetic vector at every point in space. Here we assume a collinear description owing to the high symmetry of the clusters considered. In either case, we need to consider a much larger configuration space for the electronic degrees of freedom. Another issue is the relatively localized nature of the 3d electronic states. For a real space approach, to obtain a fully converged solution, we need to employ a much finer grid spacing than for simple metals, typically 0.3 a.u. In contrast, for silicon one might use a spacing of 0.7 a.u. This finer grid required for iron results in a much larger Hamiltonian matrix and a corresponding increase in the computational load. As a consequence, while we can consider nano-crystals of silicon with over 10,000 atoms, nanocrystals of iron of this size are problematic.

The geometry of the iron clusters introduces a number of degrees of freedom. It is not currently possible to determine the definitive ground state for systems with dozens of atoms as myriads of clusters can exist with nearly degenerate energies. However, in most cases, it is not necessary to know the ground state. We are more interested in determining what structures are “reasonable” and representative of the observed ensemble, i.e., if two structures are within a few meV, these structures are not distinguishable.

We considered topologically distinct clusters, e.g., clusters of both icosahedral and BCC symmetry were explored in our work. In order to investigate the role of surface faceting, we constructed clusters with faceted and non-faceted surfaces. Faceted clusters are constructed by adding successive atomic layers around a nucleation point. Small faceted icosahedral clusters exist with sizes 13, 55, 147, and 309. Faceted BCC clusters are constructed with BCC local coordination and, differently from icosahedral ones, they do not need to be centered on an atom site.

We consider two families of cubic clusters: atom-centered or bridge-centered, respectively for clusters with nucleation point at an atom site or on the bridge between two neighboring atoms. The lattice parameter is equal to the bulk value, 2.87 A. Non-faceted clusters are built by adding shells of atoms around a nucleation point so that their distance to the nucleation point is less than a specified value. As a result, non-faceted clusters usually have narrow steps over otherwise planar surfaces and the overall shape is almost spherical. By construction, non-faceted clusters have well-defined point-group symmetries: Ih or Th for the icosahedral family, Oh for the atom-centered family, and D 4h for the bridge-centered family. Clusters constructed in that manner show low tension on the surface, making surface reconstruction less likely.

Figure 3.20   Density of states in the clusters Fe388 (a) and Fe393 (b), majority spin (upper panel), and minority spin(lower panel). For reference, the density of states in bulk iron is shown in dashed lines. The Fermi energy is chosen as energy reference.

As clusters grow in size, their properties approach the properties of bulk iron. Figure 3.20a shows the DOS for Fe388, with local BCC coordination. At this size range, the DOS assumes a shape typical of bulk iron, with a three-fold partition of the 3d bands. In addition, the cohesive energy of this cluster is only 77 meV lower than in bulk. This evidence suggests that interesting size effects will be predominantly observed in clusters smaller than Fe388. Figure 3.20b shows the DOS for Fe393, which belongs to the icosahedral family. This cluster has a very smooth DOS, with not much structure compared to Fe388 and bulk BCC iron. This is due to the icosahedral-like arrangement of atoms in Fe393. The overall dispersion of the 3d peak (4 eV for 3d and 6 eV for 3d ) is nevertheless similar in all the calculated DOS.

The magnetic moment is calculated as the expectation value of the total angular momentum:

3.20 $M = μ B ℏ [ g s 〈 S z 〉 + 〈 L z 〉 ] = μ B [ g s 2 ( n ↑ − n ↓ ) + 1 ℏ 〈 L z 〉 ]$

where gs = 2 is the electron gyromagnetic ratio. Figure 3.21 illustrates the approximately linear dependence between the magnetic moment and the spin moment, <S z>, throughout the whole size range. This results in an effective gyromagnetic ratio geff = 2.04μB/ħ, which is somewhat smaller than the gyromagnetic ratio in bulk BCC iron, 2.09μB/ħ. This is probably due to an underestimation in the orbital contribution, <L z>. In the absence of an external magnetic field, orbital magnetization arises from the spin-orbit interaction, which is included in the theory as a model potential,

3.21

where ξ = 80 meV/ħ 2 [93].

Figure 3.21   Magnetic moment versus spin moment calculated for the atom-centered BCC (“plus” signs) and bridge-centered BCC (crosses) iron clusters. Th e approximate ratio is M/<S z> = geff = 2.04 μB/ħ.

Figure 3.22 shows the magnetic moment of several clusters belonging to the three families studied: atom-centered BCC (top panel), bridge-centered BCC (middle panel), and icosahe-dral (bottom panel). Experimental data obtained by Billas and collaborators [88] are also shown. The suppression of magnetic moment as a function of size is readily observed. Also, clusters with faceted surfaces are predicted to have magnetic moments lower than other clusters with similar sizes. This is attributed to the more effective hybridization of d orbitals along the plane of the facets.

The non-monotonic behavior of the measured magnetic moment with size cluster can be attributed to the shape of the surface. Under this assumption, islands of low magnetic moment (observed at sizes 45, 85, and 188) are associated to clusters with faceted surfaces. In the icosahedral family, the islands of low magnetic moment are located around faceted clusters containing 55, 147, and 309 atoms. The first island is displaced by 10 units from the measured location. For the atom-centered and bridge-centered families, we found islands at (65, 175) and (92, 173) respectively, as indicated in Figure 3.22a and b. The first two islands are also close to the measured islands at 85 and 188. Clearly, there is no exact superposition in the location of calculated islands and measured islands. The magnetic moment was measured in clusters at 120 K [88,90]. At that temperature, vibrational modes or the occurrence of metastable configurations can shift the islands of low magnetic moment or make them more diffuse. Assuming that the non-monotonic decay of magnetic moment is dictated by the cluster shape, we also conclude that clusters with local structures different from the ones we discuss here (such as cobalt clusters with hexagonal-close packed coordination, or nickel clusters with face-centered cubic coordination) should have islands of low magnetic moment located at different “magic numbers,” according to the local atomic coordination.

Figure 3.22   Calculated magnetic moments for clusters in the atom-centered (“plus” signs, a), bridge-centered (crosses, b), and icosahedral (triangles, c) families. Experimental data [88] is shown in black diamonds with error bars. Some of the faceted and non-faceted clusters are depicted next to their corresponding data points. The dashed lines indicate the value of magnetic moment per atom in bulk iron.

#### 3.4  Conclusions

Here we reviewed tools for describing the electronic and vibra-tional properties of nanocrystals. The algorithm illustrated in this chapter replaces explicit eigenvalue calculations by an approximation of the wanted invariant subspace, obtained with the use of Chebyshev polynomial filters [32]. In this approach, only the initial self-consistent-field iteration requires solving an eigenvalue problem in order to provide a good initial subspace. In the remaining iterations, no iterative eigensolvers are involved. Instead, Chebyshev polynomials are used to refine the subspace. The subspace iteration at each step is easily one or more orders of magnitude faster than solving a corresponding eigenproblem by the most efficient eigen-algorithms. Moreover, the subspace iteration reaches self-consistency within roughly the same number of steps as an eigensolver-based approach.

We illustrated this algorithm by applying it to hydrogenated silicon nanocrystals for both electronic and vibrational modes. The largest dot we examined contained over 10,000 atoms and was ∼7 nm (Si9041H1860) in diameter. We examined the evolution of the electronic properties in these nanocrystals, which we found to assume a bulk-like configuration for dots larger than ∼5 nm. In addition, we obtained scaling relations for the ionization potential, the electron affinity, and the quasi-particle gap over the size regime of interest. We found the quasi-particle gap to approach the known bulk limit within density functional theory and suggested the remaining errors in the bulk limit occurred in the ionization potential. We did a similar calculation for vibrational modes, although the size of the nanocrystals is not as large, owing to computational issues, i.e., the need for more accurate forces, and that the vibrational modes converge to the bulk values more rapidly than the electronic states. We also examined the role of a residual charge on the vibrational modes.

We also examined the role of doping and the behavior of defect wave functions in P-doped Si nanocrystals with a diameter up to 6 nm by first-principles calculations. Our calculated HFS has very good agreement with experimental data. We found that the defect wave function has a functional form similar to the hydrogen 1s orbital. A model calculation of a hydrogen atom in a quantum well can be used to describe the defect electron. In addition, our study on the energetics of P location in the Si nanocrystals indicates that the P atom will be expelled toward the surface of the nanocrystal with diameter below a critical value of ∼2 nm.

The computational tools we outlined in this chapter are not restricted to insulating materials as is often done in computational methods targeted at large systems [102]. The filtering method we employ can equally well be applied to metallic systems and we have done so for liquid Pb [43] and for iron nanocrystals [8]. We reviewed our results for the nanocrystals of iron by examining the evolution of the magnetic moment in iron clusters containing 20 to 400 atoms using our real space pseudopotential method with damped Chebyshev filtering. Three families of clusters were studied. They were characterized by the arrangement of atoms: icosahedral, BCC centered on an atom site, and BCC centered on the bridge between two neighboring atoms. We found an overall decrease of magnetic moment as the clusters grow in size toward the bulk limit. Clusters with faceted surfaces are predicted to have magnetic moment lower than other clusters with similar size. As a result, the magnetic moments is observed to decrease as a function of size in a nonmonotonic manner, which explains measurements performed at low temperature.

The utility of this numerical approach should be widely applied to a variety of problems at the nanoscale. The method is sufficiently powerful that it can be applied to systems sufficiently large that the entire nano-regime can be examined from an isolated atom to a bulk crystal. Moreover, the method has recently been extended to include systems with partial periodicity, e.g., nanowires where the system is periodic along the axis of the nanowire [35,103].

#### Acknowledgments

This work was supported in part by the National Science Foundation under DMR-0551195 and by the U.S. Department of Energy under DE-FG02-06ER46286 and DE-FG02-06ER15760. Computational support is acknowledged from the Texas Advanced Computing Center (TACC) and the DOE National Energy Research Scientific Computing Center (NERSC).

#### References

L. Brus , J. Phys. Chem. 98, 3575 (1994).
A. P. Alivisatos , Science 271, 933 (1996).
A. P. Alivisatos , J. Phys. Chem. 100, 13226 (1996).
X. G. Peng , L. Manna , W. D. Yang , J. Wickham , E. Scher , A. Kadavanich , and A. P. Alivisatos , Nature 404, 59 (2000).
M. Law , J. Goldberger , and P. Yang , Ann. Rev. Mater. Res. 34, 83 (2004).
C. Burda , X. B. Chen , R. Narayanan , and M. A. El-Sayed , Chem. Rev. 105, 1025 (2005).
R. Walters , G. Bourianoff , and H. Atwater , Nat. Mater. 4, 143 (2005).
M. Tiago , Y. Zhou , M. M. G. Alemany , Y. Saad , and J. R. Chelikowsky , Phys. Rev. Lett. 97, 147201 (2006).
G. Rollmann , M. E. Gruner , A. Hucht , P. Entel , M. L. Tiago , and J. R. Chelikowsky , Phys. Rev. Lett. 99, 083402 (2007).
J. R. Chelikowsky , J. Phys. D: Appl. Phys. 33, R33 (2000).
L. W. Wang and A. Zunger , J. Chem. Phys. 100, 2394 (1994).
L. W. Wang and A. Zunger , J. Phys. Chem. 98, 2158 (1994).
W. Kohn and L. J. Sham , Phys. Rev. 140, A1133 (1965).
P. Hohenberg and W. Kohn , Phys. Rev. 136, B864 (1964).
D. M. Ceperley and B. J. Alder , Phys. Rev. Lett. 45, 566 (1980).
C. G. Broyden , Math. Comp. 19, 577 (1965), ISSN 00255718, URL http://www.jstor.org/stable/2003941.
J. R. Chelikowsky and M. L. Cohen , in Handbook of Semiconductors, edited by T. S. Moss and P. T. Landsberg (Elsevier, Amsterdam, the Netherlands, 1992), p. 59.
M. L. Cohen and J. R. Chelikowsky , Electronic Structure and Optical Properties of Semiconductors, 2nd edn. (Springer-Verlag, Berlin, Germany, 1989).
D. R. Hamann , M. Schluter , and C. Chiang , Phys. Rev. Lett. 43, 1494 (1979).
G. P. Kerker , J. Phys. C 13, L189 (1980).
N. Troullier and J. Martins , Phys. Rev. B 43, 1993 (1991).
G. Bachelet , D. R. Hamann , and M. Schluter , Phys. Rev. B 26, 4199 (1982).
G. Kresse and J. Furthmuller , Phys. Rev. B 54, 11169 (1996).
J. R. Chelikowsky , N. Troullier , and Y. Saad , Phys. Rev. Lett. 72, 1240 (1994).
E. L. Briggs , D. J. Sullivan , and J. Bernholc , Phys. Rev. B 52, R5471 (1995).
G. Zumbach , N. A. Modine , and E. Kaxiras , Solid State Commun. 99, 57 (1996).
J.-L. Fattebert and J. Bernholc , Phys. Rev. B 62, 1713 (2000).
T. L. Beck , Rev. Mod. Phys. 74, 1041 (2000).
M. Heikanen , T. Torsti , M. J. Puska , and R. M. Nieminen , Phys. Rev. B 63, 245106 (2001).
T. Torsti , T. Eirola , J. Enkovaara , T. Hakala , P. Havu , V. Havu , T. Hoynalanmaa , J. Ignatius , M. Lyly , I. Makkonen et al., Physica Status Solidi (b) 243, 1016 (2006).
L. Kronik , A. Makmal , M. L. Tiago , M. M. G. Alemany , M. Jain , X. Huang , Y. Saad , and J. R. Chelikowsky , Physica Status Solidi (b) 243, 1063 (2006).
Y. Zhou , Y. Saad , M. L. Tiago , and J. R. Chelikowsky , Phys. Rev. E 74, 066704 (2006).
M. M. G. Alemany , M. Jain , J. R. Chelikowsky , and L. Kronik , Phys. Rev. B 69, 075101 (2004).
M. Alemany , M. Jain , M. L. Tiago , Y. Zhou , Y. Saad , and J. R. Chelikowsky , Comp. Phys. Commun. 177, 339 (2007).
A. Natan , A. Mor , D. Naveh , L. Kronik , M. L. Tiago , S. P. Beckman , and J. R. Chelikowsky , Phys. Rev. B 78, 075109 (2008).
B. Fornberg and D. M. Sloan , Acta Numerica 94, 203 (1994).
L. Kleinman and D. M. Bylander , Phys. Rev. Lett. 48, 1425 (1982).
J. C. Phillips and L. Kleinman , Phys. Rev. 116, 287 (1959).
A. Stathopoulos , S. Ogut , Y. Saad , J. Chelikowsky , and H. Kim , Comput. Sci. Eng. 2, 19 (2000).
C. Bekas , Y. Saad , M. L. Tiago , and J. R. Chelikowsky , Comp. Phys. Commun. 171, 175 (2005).
R. Lehoucq , D. C. Sorensen , and C. Yang , ARPACK Users’ Guide:Solution of Large Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods (SIAM, Philadelphia, PA, 1998).
K. Wu , A. Canning , H. D. Simon , and L.-W. Wang , J. Comp. Phys. 154, 156 (1999).
M. M. G. Alemany , R. C. Longo , L. J. Gallego , D. J. Gonzalez , L. E. Gonzalez , M. L. Tiago , and J. R. Chelikowsky , Phys. Rev. B 76, 214203 (2007).
S. Furukawa and T. Miyasato , Phys. Rev. B 38, 5726 (1988).
T.-L. Chan , M. L. Tiago , E. Kaxiras , and J. R. Chelikowsky , Nano Lett. 8, 596 (2008).
S. Ogut , J. R. Chelikowsky , and S. G. Louie , Phys. Rev. Lett. 79, 1770 (1997).
J. P. Perdew and M. Levy , Phys. Rev. Lett. 51, 1884 (1983).
L. J. Sham and M. Schluter , Phys. Rev. Lett. 51, 1888 (1983).
M. L. Tiago and J. R. Chelikowsky , Phys. Rev. B 73, 205334 (2006).
M. S. Hybertsen and S. G. Louie , Phys. Rev. B 34, 5390 (1986).
M. S. Hybertsen and S. G. Louie , Phys. Rev. B 35, 5585 (1987).
K. Parlinski , Z. Q. Li , and Y. Kawazoe , Phys. Rev. Lett. 78, 4063 (1997).
A. V. Postnikov , O. Pages , and J. Hugel , Phys. Rev. B 71, 115206 (2005).
X. Huang , E. Lindgren , and J. R. Chelikowsky , Phys. Rev. B 71, 165328 (2005).
J. Song , S. E. Ulloa , and D. A. Drabold , Phys. Rev. B 53, 8042 (1996).
X. Jing , N. Troullier , J. R. Chelikowsky , K. Wu , and Y. Saad , Solid State Commun. 96, 231 (1995).
M. R. Pederson , K. Jackson , D. V. Porezag , Z. Hajnal , and T. Frauenheim , Phys. Rev. B 65, 2863 (1996).
A. Valentin , J. See , S. Galdin-Retailleau , and P. Dollfus , J. Phys. 92, 1 (2007).
C. Meier , S. Luttjohann , V. G. Kravets , H. Nienhaus , A. Lorke , and H. Wiggers , Physica E 32, 155 (2006).
X .-S. Zhao , Y.-R. Ge , and X. Zhao , J. Mater. Sci. 33, 4267 (1998).
B. Silvi and A. Savin , Nature 371, 683 (1994).
B. G. Streetman and S. Banerjee , Solid State Electronic Devices, 5th edn. (Prentice Hall, Englewood Cliffs, NJ, 2000).
A. D. Yoffe , Adv. Phys. 50, 1 (2001).
Y. Cui and C. M. Lieber , Science 291, 851 (2001).
D. Appell , Nature 419, 553 (2002).
A. Mimura , M. Fujii , S. Hayashi , D. Kovalev , and F. Koch , Phys. Rev. B 62, 12625 (2000).
G. Mauckner , W. Rebitzer , K. Thonke , and R. Sauer , Physica Status Solidi (b) 215, 871 (1999).
J. Muller , F. Finger , R. Carius , and H. Wagner , Phys. Rev. B 60, 11666 (1999).
B. J. Pawlak , T. Gregorkiewicz , C. A. J. Ammerlaan , and P. F. A. Alkemade , Phys. Rev. B 64, 115308 (2001).
M. Fujii , A. Mimura , S. Hayashhi , Y. Yamamoto , and K. Murakami , Phys. Rev. Lett. 89, 206805 (2002).
G. Feher , Phys. Rev. 114, 1219 (1959).
D. Melnikov and J. R. Chelikowsky , Phys. Rev. Lett. 92, 046802 (2004).
Z. Zhou , M. L. Steigerwald , R. A. Friesner , L. Brus , and M. S. Hybertsen , Phys. Rev. B 71, 245308 (2005).
G. Cantele , E. Degoli , E. Luppi , R. Magri , D. Ninno , G. Iadonisi , and S. Ossicini , Phys. Rev. B 72, 113303 (2005).
S. Ossicini , E. Degoli , F. Iori , E. Luppi , R. Magri , G. Cantele , F. Trani , and D. Ninno , Appl. Phys. Lett. 87, 173120 (2005).
C. Y. Fong , H. Zhong , B. M. Klein , and J. S. Nelson , Phys. Rev. B 49, 7466 (1994).
I. H. Lee , K. H. Ahn , Y. H. Kim , R. M. Martin , and J. P. Leburton , Phys. Rev. B 60, 13720 (1999).
M. Lannoo , C. Delerue , and G. Allan , Phys. Rev. Lett. 74, 3415 (1995).
M. Fujii , K. Toshikiyo , Y. Takase , Y. Yamaguchi , and S. Hayashi , J. Appl. Phys. 94, 1990 (2003).
G. Allan , C. Delerue , M. Lannoo , and E. Martin , Phys. Rev. B 52, 11982 (1995).
J. A. Weil and J. R. Bolton , Electron Paramagnetic Resonance: Elementary Theory and Practical Applications, 2nd edn. (Wiley, Hoboken NJ, 2007).
C. G. V. de Walle and P. E. Blochl , Phys. Rev. B 47, 4244 (1993).
D. R. Penn , Phys. Rev. B 128, 2093 (1962).
L. E. Brus , J. Chem. Phys. 79, 4403 (1983).
L. E. Brus , J. Chem. Phys. 79, 5566 (1983).
D. Mattis , The Theory of Magnetism, 2nd edn. (Springer-Verlag, Berlin, Germany, 1988).
C. Bean and J. Livingston , J. Appl. Phys. 30, 120s (1959).
I. Billas , J. Becker , A. Chatelain , and W. de Heer , Phys. Rev. Lett. 71, 4067 (1993).
D. Gerion , A. Hirt , I. Billas , A. Chatelain , and W. de Heer , Phys. Rev. B 62, 7491 (2000).
I. Billas , A. Chatelain , and W. de Heer , Science 265, 1682 (1994).
J. Franco , A. Vega , and F. Aguilera-Granja , Phys. Rev. B 60, 434 (1999).
A. Postnikov and P. Entel , Phase Transitions 77, 149 (2004).
O. Sipr , M. Kosuth , and H. Ebert , Phys. Rev. B 70, 174423 (2004).
A. Postnikov , P. Entel , and J. Soler , Eur. Phys. J. D 25, 261 (2003).
R. Felix-Medina , J. Dorantes-Davila , and G. Pastor , Phys. Rev. B 67, 094430 (2003).
K. Edmonds , C. Binns , S. Baker , S. Thornton , C. Norris , J. Goedkoop , M. Finazzi , and N. Brookes , Phys. Rev. B 60, 472 (1999).
G. Pastor , J. Dorantes-Davila , and K. Bennemann , Phys. Rev. B 40, 7642 (1989).
O. Dieguez , M. Alemany , C. Rey , P. Ordejon , and L. Gallego , Phys. Rev. B 63, 205407 (2001).
P. Bruno , in Magnetismus von Festkorpern undgrenzflachen, edited by P. Dederichs , P. Grunberg , and W. Zinn (IFF-Ferienkurs, Forschungszentrum Julich, Germany, 1993), pp. 24.1–24.27.
J. P. Perdew , K. Burke , and Y. Wang , Phys. Rev. B 54, 16533 (1996).
J. R. Chelikowsky , N. Troullier , K. Wu , and Y. Saad , Phys. Rev. B 50, 11355 (1994).
M. J. Gillan , D. R. Bower , A. S. Torralba , and T. Miyazaki , Comp. Phys. Commun. 177, 14 (2007).
J. Han , M. L. Tiago , T.-L. Chan , and J. R. Chelikowsky , J. Chem. Phys. 129, 144109 (2008).
U. Itoh , Y. Toyoshima , H. Onuki , N. Washida , and T. Ibuki , J. Chem. Phys. 85, 4867 (1986).

## Use of cookies on this website

We are using cookies to provide statistics that help us give you the best experience of our site. You can find out more in our Privacy Policy. By continuing to use the site you are agreeing to our use of cookies.