text:
# Comparison of localization procedures for applications in crystal embedding

## I Introduction

Electronic structure calculations of extended systems with a local perturbation, such as point defects in the bulk of crystals Marshall-book  or adsorption of molecules at their surfaces Chemisorption-Reactivity-1997  are of fundamental importance in solid state physics and chemistry. Over the last decade a number of effective computational techniques have been developed to study the electronic ground state of such systems which are based on periodic boundary conditions (PBC) and either Hartree-Fock (HF) CRYSTAL98 ; CRYSTAL  or Density Functional Theory (DFT) Payne-DFT ; Parr-Yang  approaches. In these methods a local perturbation (e.g. an adsorbed molecule together with a fragment of a crystal surface) is artificially periodically repeated in a cell which is large enough to ensure that the interaction between periodic images is negligible. These methods can also be applied to study extended (but not infinite) biological systems in which important chemistry is usually associated with a local part of the entire molecule(s) Segall-2002 .

Another set of methods, commonly referred to as embedding techniques, originate from a model in which a single local perturbation is considered in the direct space of the entire system. This makes the model closer to reality at low concentration, but, at the same time, it makes it more challenging since, due to the lack of periodic symmetry, well developed PBC based techniques cannot be applied here. Instead, a number of *hybrid* methods have been developed which treat different parts of the system at different levels of the theory. Most of these methods combine *ab initio* quantum mechanics (QM) methods (based either on the DFT or the HF methods and their extensions) applied to a finite fragment of the system (a *quantum cluster*), with molecular mechanics (MM) methods based on semi-classical force fields and applied to the rest of the system (*environment region*). The idea is to consider the most relevant part of the system (e.g. with respect to a process in question) in great detail, while the rest of the system is treated at a substantially lower level. These methods, developed mostly within the quantum chemistry community, and usually referred to as QM/MM have proven to be extremely successful QM/MM ; Sauer-Sierka-2000 ; Hall-Hinde-Burton-Hillier-2000 ; Rivail3 ; Murphy-Philipp-Freisner-2000 . Note that some other embedding schemes EMC-1 ; EMC-2 ; Barandiaran-1996 ; Bredow-1999 ; Petja-2000 ; Sulimov2002 ; Nasluzov2001  developed mostly in the solid state community are very close in spirit to the QM/MM methods.

Almost in all embedding methods mentioned above the quantum cluster is surrounded by point charges of the MM region. In addition, in covalent systems the bonds coming out of the cluster are usually terminated by pseudoatoms (see, e.g. Sulimov2002 ), so-called link atoms Rivail3 ; Sauer-Sierka-2000 ; Murphy-Philipp-Freisner-2000  or pseudopotentials Abarenkov-Tupitsyn-2001 ; Abarenkov-Tupitsyn-2001r ; Petja-2000 ; Nasluzov2001 ; Sulimov2002 .

Another class of embedding schemes relies on a more “electronic” (less “mechanical”) representation of the environment region surrounding the quantum cluster. For instance, in Vreven-Morokuma-2000 ; Sauer-Sierka-2000 ; Abarenkov-Bulatov-1997  this is achieved by a special total energy construction which allows a combination of several electronic structure methods of different complexity applied to different parts of the system; other methods EMC-1 ; EMC-2 ; Barandiaran-1996  rely on a representation of the wavefunction of the whole system as an antisymmetrised product of strongly orthogonal many-electron group functions (see, e.g. McWeeny ; McWeeny-rev ) associated with atoms, bonds or molecules depending on the specific type of chemical bonding in the system.

We believe that the formalism based on group functions is the most appropriate one for the derivation of any embedding scheme. A rather general method based on overlapping (not strongly orthogonal) group functions my-AD-1 ; my-AD-2  is presently being developed in our laboratory. Our method which is similar in spirit to some one-electron methods Shidlovskaya-2002 ; Fornili-Sironi-Raimondi-2003 ; Mo-Gao-2000  is based on construction of strongly localized orbitals which are designed to represent the true electronic density of the entire system via a combination of elementary densities associated in simple cases with atoms, ions and/or bonds. Our initial effort in this project is focused on the development of an embedding scheme based on the HF approximation and applied to point defects in the bulk or at surfaces of periodic crystals. Our intention is to create a rather general technique which can be valid for systems of different chemical character, ranging from purely ionic to strongly covalent (note that our method cannot be directly applied to metals). Therefore, the proper choice of the localization technique which can deliver localized orbitals across a wide range of systems with various character of chemical bonding is crucial for our method to work for those systems.

It is the main objective of the present paper to critically analyze and develop further a number of localization methods in order to verify their ability to describe a wide range of different chemical bondings in *periodic crystals*. Two systems are considered in the present paper in detail, MgO and Si bulk crystals, which are examples of extreme ionic and covalent bonding, respectively. Note, however, that the method we suggest is not limited to periodic systems and, with some insignificant modifications, can also be applied e.g. to infinite amorphous and finite biological systems. The application of the present method to those systems will be a matter of future publications.

It is relevant to mention, as far as the localization methods are concerned, that there are several methods developed Marzari97 ; Dovesi-Saunders2001  for obtaining *orthogonal* localized orbitals (i.e. Wannier functions) out of the Bloch-like solutions of the HF or Kohn-Sham (KS) equations Parr-Yang  for periodic crystals. Due to the built-in orthogonality even strongly localized Wannier functions have long-range tails which make these functions nontransferable to other systems, e.g. when a chemical bond between the same species is placed in a different chemical environment. That is why our interest is focused on construction of *non-orthogonal localized orbitals* which do not have this disadvantage and thus are more appropriate for our purposes.

The plan of the paper is the following. Our philosophy in constructing localized orbitals as well as a short overview of existing localization methods is given in section II with special emphasis on the methods used in our present work. All the necessary notations are also introduced there. In section 3 we describe our implementation of some of the methods and their application to MgO and Si crystals. The paper is finished with a short discussion and conclusions in section 4.

## II Localization methods

### II.1 General philosophy

In order to describe the crystal as a set of overlapping localized functions, $`\left\{\stackrel{~}{\phi }_a(𝐫)\right\}`$, which are given as a linear combination of the original canonical set $`\left\{\phi _i^c(𝐫)\right\}`$ (and which thus span the same occupied Fock space), one has first to identify the regions of space where each of the functions $`\stackrel{~}{\phi }_a(𝐫)`$ have to be localized. This question can be viewed as purely technical since, any linear combination of the canonical set will give the same electron density $`\rho (𝐫)`$. We, however, adopt in our work a strategy based on the chemistry of the system in question. Namely, the choice of the localization regions is based on the type of the chemical bonding, e.g. on atoms in the cases of atomic or ionic systems, on two atoms in the case of covalent bonding, etc. A more complicated choice may be necessary in the cases of intermediate bonding. Several different nonequivalent regions may be necessary to represent a crystal unit cell which can then be periodically translated to reproduce the whole infinite crystal. Note that there could be several localized orbitals associated with every such a region forming together an electronic group McWeeny ; EMC-1 . For instance, in the case of the Si crystal one needs four localized regions associated with four bonds; each bond is represented by a single double occupied localized orbital.

Once the occupied Fock space is obtained via a set of canonical orbitals and localized regions are identified, it is necessary to find such linear combination of canonical orbitals which are localized in each of the regions. The topic of construction of localized (non-canonical) molecular orbitals (MO) out of delocalized canonical solutions of the HF or Kohn-Sham equations is an old one Loc-Deloc  and many methods have since been developed.

Let us assume that a canonical solution of the restricted HF equations for the entire system (a closed shell crystal) is known McWeeny :

$$\widehat{F}\phi _i^c(𝐫)=\epsilon _i\phi _i^c(𝐫)$$
(1)
$$\phi _i^c(𝐫)=\underset{\mu }{}C_{\mu i}^c\chi _\mu (𝐫)$$
(2)
where $`\widehat{F}`$ is the Fock operator, $`\phi _i^c(𝐫)`$ is a spin-independent canonical MO (CMO) which is expanded over a set of atomic orbitals (AOs) $`\chi _\mu (𝐫)`$. The electronic density of the system

$$\rho (𝐫)=2\underset{i}{\overset{occ}{}}\left|\phi _i^c(𝐫)\right|^2$$
(3)
contains the summation only over occupied CMOs thus ensuring the correct normalization to the number $`N`$ of the electrons in the system. If an arbitrary (generally *non-unitary*) transformation $`𝐔=U_{aj}`$ of the CMOs within the occupied subspace is performed,

$$\stackrel{~}{\phi }_a(𝐫)=\underset{j}{\overset{occ}{}}U_{aj}\phi _j^c(𝐫)\underset{\mu }{}\stackrel{~}{C}_{\mu a}\chi _\mu (𝐫)$$
(4)
than the expression for the density via the new set of orbitals should contain the inverse of the overlap matrix $`\stackrel{~}{𝐒}=\stackrel{~}{S}_{ab}`$ McWeeny :

$$\rho (𝐫)=2\underset{ab}{\overset{occ}{}}\stackrel{~}{\phi }_a(𝐫)\left(\stackrel{~}{S}\right)_{ab}^1\stackrel{~}{\phi }_b^{}(𝐫)$$
(5)
where $`\stackrel{~}{S}_{ab}=\stackrel{~}{\phi }_a(𝐫)|\stackrel{~}{\phi }_b(𝐫)`$ is the overlap integral. Note that the two representations of the electron density, Eqs. (3) and (5) are absolutely equivalent. Moreover, any linear combination (4) of the occupied CMOs leads to the same density. If the transformation is unitary, then the overlap matrix is a unity matrix and the density takes on its “diagonal“ form (3).

In general, any localization procedure is equivalent to some transformation $`𝐔`$ of the CMOs. Suppose, we would like to obtain $`n`$ localized MOs (LMOs) in some region $`A`$. To find the necessary transformation, one can formulate an optimization (minimization or maximization) problem for some specific *localizing functional* $`\stackrel{~}{\mathrm{\Omega }}_A\left[\left\{\stackrel{~}{\phi }_a\right\}\right]`$ with the constraint that the LMOs associated with region $`A`$ are orthonormal (of course, LMOs associated with different regions will not be orthogonal in general). We shall limit ourselves with such functionals which are invariant under arbitrary unitary transformations of LMOs, i.e. which in fact depend on the orbitals $`\left\{\stackrel{~}{\phi }_a\right\}`$ via invariants in the form of the non-diagonal “density”

$$\sigma _A(𝐫,𝐫^{})=\underset{a=1}{\overset{n}{}}\stackrel{~}{\phi }_a(𝐫)\stackrel{~}{\phi }_a^{}(𝐫^{})$$
(6)
constructed out of the LMOs associated with region $`A`$, i.e. $`\stackrel{~}{\mathrm{\Omega }}_A\left[\left\{\stackrel{~}{\phi }_a\right\}\right]\mathrm{\Omega }_A\left[\sigma _A\right]`$. We shall see in a moment that his requirement ensures an existence of a simple eigenvalue-like problem for the LMOs. Note in passing that some other types of functionals are also sometimes used which do not fall within this category. For instance, Admiston and Ruedenberg proposed to find the maximum of the self-repulsion energy Edmiston-Ruenberg , while later on von Niessen suggested to maximize the charge density overlap functional vonNiessen71 . Since the mentioned functionals are not invariant under unitary transformations of LMOs and are also quite expensive computationally, we do not consider them in the following. The quantity $`\sigma _A(𝐫,𝐫^{})`$ will be referred to in the following as the *region electron density* or region density for short.

To obtain all $`n`$ LMOs associated with region $`A`$, an optimum of the following functional is sought for:

$$\mathrm{\Omega }_A^{}\left[\sigma _A\right]=\mathrm{\Omega }_A\left[\sigma _A\right]\underset{a,b=1}{\overset{n}{}}\xi _{ab}\left(\stackrel{~}{\phi }_a|\stackrel{~}{\phi }_b\delta _{ab}\right)$$
(7)
where $`\xi _{ab}`$ are the corresponding Lagrangian multipliers. Because the actual dependence of the functional $`\mathrm{\Omega }_A\left[\sigma _A\right]`$ on the orbitals is built-in via the region density (6), the functional derivative $`\frac{\delta \mathrm{\Omega }_A}{\delta \stackrel{~}{\phi }_a^{}(𝐫)}`$ can always be written using an operator $`\widehat{\mathrm{\Omega }}_A(𝐫)`$ defined through an identity $`\frac{\delta \mathrm{\Omega }_A}{\delta \stackrel{~}{\phi }_a^{}(𝐱)}=\widehat{\mathrm{\Omega }}_A\stackrel{~}{\phi }_a(𝐱)`$ since $`\frac{\delta \sigma _A(𝐫,𝐫^{})}{\delta \stackrel{~}{\phi }_a^{}(𝐱)}=\stackrel{~}{\phi }_a(𝐫)\delta (𝐱𝐫^{})`$. Examples illustrating this point will be given below. We shall refer to the operator $`\widehat{\mathrm{\Omega }}_A(𝐫)`$ as the *localization operator* in the following. An important property of the localization operator is that it can also be considered as a functional of the region density (6), i.e. it preserves the invariance property of the localizing functional it is built from.

Using standard methods, i.e. setting the variational derivative of the functional (7) with respect to the orbital $`\stackrel{~}{\phi }_a^{}(𝐫)`$ to zero and then performing a unitary transformation of the LMOs which diagonalizes the matrix of Lagrangian multipliers, one can easily obtain the following equations for the LMOs sought for:

$$\widehat{\mathrm{\Omega }}_A\stackrel{~}{\phi }_a(𝐫)=\lambda _a\stackrel{~}{\phi }_a(𝐫)$$
(8)
or, in the matrix form,

$$\underset{j}{\overset{occ}{}}\mathrm{\Omega }_{ij}^AU_{aj}=\lambda _aU_{ai}$$
(9)
where the matrix $`\mathrm{\Omega }_{ij}^A`$ is given via matrix elements of the operator $`\widehat{\mathrm{\Omega }}_A`$ calculated using canonical orbitals $`\phi _i^c(𝐫)`$ and $`\phi _j^c(𝐫)`$.

Equations (8) resemble an eigenvalue problem for the operator $`\widehat{\mathrm{\Omega }}_A`$. Note, however, that in some cases the localization operator may still depend on the region density and thus on the orbitals themselves. Therefore, similarly to the HF or Kohn-Sham problem, the system of equations (8) should be solved self-consistently.

The eigenvalue problem (8) or (9) may give a set of solutions from which only the first (in the case when $`\mathrm{\Omega }_A`$ is minimized) or the last (maximized) $`n`$ solutions should be chosen. If the localization criterion (i.e. the functional $`\mathrm{\Omega }_A`$) used is appropriate, then (i) the chosen $`n`$ solutions would have close eigenvalues $`\lambda _a`$ which corresponds to their similar localization in region $`A`$, and (ii) the gap in the eigenvalues $`\lambda _a`$ between the chosen $`n`$ and other solutions is considerable, i.e. the other solutions have much worse localization in region $`A`$ (cf. Whitt-Pakk ). By collecting LMOs from all regions in the unit cell and then translating those over the whole crystal it should be possible to span the whole occupied Fock space and thus construct the total electron density (5).

One point is in order now. So far we have assumed that the set of canonical MOs which span the occupied part of the Fock space is already known. In other words, the procedure consists of two steps: firstly, a HF (or Kohn-Sham) problem is solved and thus the occupied Fock space is determined, and, secondly, the LMOs are obtained by finding an appropriate linear combinations of the canonical orbitals within this space. However, it is also possible to formulate the problem in such a way that LMOs are obtained together with the set of canonical orbitals in a single step Loc-Deloc . In this method a localization criterion is considered alongside the energy minimization leading to a set of so-called Adams-Gilbert (AG) equations (see, e.g. Gilbert ) which are solved in a self-consistent manner. For instance, a projection operator on the subspace of the LMOs was used by Stoll *et. al.* Stoll1980  as the specific localization method. This technique was implemented in Shidlovskaya-2002  for the embedded molecular cluster calculations. The LMOs resulting from a single AG calculation are orthogonal as solutions of a single secular problem. The first eigenvectors obtained will show strong localization within the chosen region $`A`$; other eigenvectors will be much less localized and can usually be distinguished by a gap in their eigenvalues as explained above. To obtain LMOs strongly localized in a different region $`A^{}`$, one has to solve the AG equations once again using another localization criterion and then pick up the necessary number of the most localized orbitals. Repeating this procedure across the entire system, the whole occupied Fock space can be split into sets of mutually non-orthogonal LMOs. Of course, in the case of the perfect crystals this procedure should only be applied to various localization regions within the primitive unit cell owing to crystals periodic symmetry.

Note that it is also possible to obtain all the LMOs corresponding to several localization regions at once within the same self-consistent calculation by solving the necessary sets of eigenproblems associated with each region Fornili-Sironi-Raimondi-2003 ; Mo-Gao-2000 . The LMOs obtained using this technique are known as extremely localized MOs. This method is quite expensive computationally since the overlap between LMOs localized in different regions in space changes in the course of the iteration procedure and this affects the convergence. However, as will be discussed in the following, this technique is a logical step forward in the future development of our method to be presented in section 3.

There are many ways in which the localizing functional $`\mathrm{\Omega }_A`$ can be chosen. Some of these methods which will be utilized in the present work will be considered below in more detail.

### II.2 Methods based on functionals linear in region density

In a number of methods Mayer1996  the localizing functional is proportional to the non-diagonal density (6) and thus can be represented as a Hermitian bilinear functional with respect to the LMOs of the following general form:

$$\mathrm{\Omega }_A=\left[\widehat{\mathrm{\Omega }}_A\sigma _A(𝐫,𝐫^{})\right]_{𝐫^{}𝐫}𝑑𝐫=\underset{a=1}{\overset{n}{}}\stackrel{~}{\phi }_a^{}(𝐫)\widehat{\mathrm{\Omega }}_A\stackrel{~}{\phi }_a(𝐫)𝑑𝐫\underset{a=1}{\overset{n}{}}\underset{jk}{\overset{occ}{}}U_{aj}^{}\mathrm{\Omega }_{jk}^AU_{ak}$$
(10)
where $`\widehat{\mathrm{\Omega }}_A`$ is some localization operator and the Hermitian matrix $`𝛀^A=\mathrm{\Omega }_{jk}^A`$ can easily be written in terms of the canonical MOs using the definition (2):

$$\mathrm{\Omega }_{jk}^A=\phi _j^c\left|\widehat{\mathrm{\Omega }}_A\right|\phi _k^c=\underset{\mu ,\nu }{}C_{\mu j}^cC_{\nu k}^c\chi _\mu \left|\widehat{\mathrm{\Omega }}_A\right|\chi _\nu $$
(11)
For all methods of this group both the operator $`\widehat{\mathrm{\Omega }}_A`$ and the matrix $`𝛀^A`$ do not depend on the LMOs sought for so that in order to obtain the localized orbitals one has simply to find the eigenvectors of the matrix $`𝛀^A`$ using Eq. (9). Two methods of this group are implemented in our work and will be considered in the following in more detail.

#### Mulliken’s net population (method M)

Magnasco-Perico criterion maximizes Mulliken’s Mulliken  net atomic population produced by the LMOs in the selected region Magnasco-Perico ; Mayer1996 . In this case the matrix $`𝛀^A`$ is chosen in the following form:

$$\mathrm{\Omega }_{jk}^A=\underset{\mu ,\nu A}{}C_{\mu j}^cS_{\mu \nu }C_{\nu k}^c$$
(12)
where $`S_{\mu \nu }`$ is the overlap integral between two AOs $`\chi _\mu `$ and $`\chi _\nu `$. The summation here is performed over AOs which are centered in the chosen region $`A`$. Thus, in practice the localization region in this method is specified by a selection of AOs in Eq. (12). This way one can make the LMOs to have the maximum contribution from the specified AOs in region $`A`$. Sometimes a different choice of AOs may lead to physically identical localization (see the next section). This method will be referred to as method M.

#### The projection on the atomic subspace (method P)

The Roby’s population maximization Roby74  gives LMOs for which the projection on the subspace spanned by the basis orbitals centered within the selected region is a maximum, or is at least stationary Mayer1996 . In this method the localization operator $`\widehat{\mathrm{\Omega }}_A`$ in Eq. (10) is chosen in the form of a projection operator:

$$\widehat{\mathrm{\Omega }}_A=\underset{\mu ,\nu A}{}|\chi _\mu (𝐒_A^1)_{\mu \nu }\chi _\nu |$$
(13)
where $`𝐒_A^1`$ stands for the inverse of the overlap matrix $`𝐒_A`$ defined on all AOs $`\mu ,\nu A`$. Note that operator $`\widehat{\mathrm{\Omega }}_A`$ is idempotent: $`\left(\widehat{\mathrm{\Omega }}_A\right)^2=\widehat{\mathrm{\Omega }}_A`$. It projects any orbital into a subspace spanned by the AOs associated with region $`A`$ only. In particular, $`\widehat{\mathrm{\Omega }}_A|\chi _\mu =|\chi _\mu `$. The detailed expression for the matrix $`𝛀^A`$ is then:

$$\mathrm{\Omega }_{jk}^A=\underset{\lambda ,\tau }{}C_{\lambda j}^cC_{\tau k}^c\left[\underset{\mu ,\nu A}{}S_{\lambda \mu }(𝐒_A^1)_{\mu \nu }S_{\nu \tau }\right]$$
(14)
Here the first double summation is performed over all AOs of the system. Region $`A`$ is also defined via a subset of AOs: by choosing particular AOs one ensures the maximum overlap of the LMOs with them. It is seen that this method, which will be referred to as method P, although different in the implementation, is very similar in spirit to the previous method M.

#### Other methods

Note that several other methods Mayer1996  also belong to this class of methods. Since we are not using them here, we shall only mention some of them. Bader’s method (see also Bader ) is computationally expensive and leads to LMOs with discontinuities at the boarder of the localization regions. The widely used Pipek-Mezey localization scheme Pipek-Mezey  could be described as the maximization of the Mulliken’s gross atomic population. The Pipek-Mezey functional corresponds also to a minimization of the number of atoms over which the LMO is to be spread. This method is very similar to the Mulliken’s net population method considered above, although is slightly more computationally expensive. All population methods have an advantage of being very simple in the implementation which results in fast non-self-consistent algorithms: indeed, only overlap integrals are to be computed. Note that instead of the overlap, one can also maximize an exchange interaction of the LMOs with a set of AOs in region $`A`$ Whitt-Pakk .

Perhaps the most widely used, due to its relatively low computational cost, is the Foster-Boys FosterBoys  method in which the dipole moment matrix element between so-called exclusive orbitals is maximized. The efficiency of the HF method was improved in Anikin2003  by using localized orbitals constructed from the Foster-Boys method as AOs. Recently Dovesi-Saunders2001  the Wannier functions were calculated for periodic Si and MgO crystals using the modified version of the Foster-Boys method which ensured better localization within the cell volume. Note that, according to Ref. Mayer1996 , the Pipek-Mezey functional, unlike the Fosters-Boys method, preserves the $`\sigma /\pi `$ separation of double bonds, which is in chemistry usually preferred over the $`\tau `$ picture (where the orbitals are proportional to the linear combinations $`\sigma +\pi `$ and $`\sigma \pi `$) associated with the Fosters-Boys localization procedure.

### II.3 Methods based on functionals bilinear in region density

More complicated localization procedures can be constructed if the localizing functional is bilinear in non-diagonal density (6) of region $`A`$ in question, i.e. is of the fourth order with respect to the LMOs sought for. Three important general points should be mentioned: (i) for all methods of this group the localizing operator $`\widehat{\mathrm{\Omega }}_A`$ is linear in the density $`\sigma _A(𝐫,𝐫^{})`$ and (ii) is thus invariant under any unitary transformation of the LMOs; (iii) therefore, one still has the secular problem (8) or (9) for the LMOs in this case; however, it is to be solved *self-consistently*.

#### Minimization of the HF energy of a structure element (method E)

In this paper we shall only apply one method of this group in which the functional $`\mathrm{\Omega }_A`$ is chosen as the HF energy of a finite fragment (structure element (SE), cf. EMC-1 ) of the system. This method, which will be referred to as method E in the following, in its simplest version of the localization on a single atom originates from Adams Adams . This method was recently used in Abarenkov2004  to derive an embedding potential provided by a part of a molecule. Method E is based on an intuitive idea that every *finite* system, e.g. a SE, will try to find an energetically favorable ground state which will be localized in space. The SE could be an atom, group of atoms or a bond. In the latter case the SE for Si represents essentially a hydrogen-like molecule consisting of fragments of two nearest atomic cores each of charge $`+e`$ ($`e`$ is the electron charge) and two electrons of opposite spins (see the next section).

The eigenvalue problem (9) in this case is nothing but the usual Hartree-Fock-Roothaan (HFR) problem McWeeny  for molecular orbitals of region $`A`$, i.e. the elements of the matrix $`𝛀^A`$ are matrix elements

$$\mathrm{\Omega }_{jk}^A=\phi _j^c\left|\widehat{\mathrm{\Omega }}_A\right|\phi _k^c$$

where $`\widehat{\mathrm{\Omega }}_A\widehat{F}_A`$ is the usual HF operator of the SE containing both electron-electron and electron-core interactions. The peculiarity of this case is that the MOs (which are, in fact, LMOs) are expanded not via AOs but rather via canonical HF orbitals for the whole system which, in turn, are linear combinations of AOs of the whole system as in Eq. (2). In other words, the difference with the usual setup of the HFR problem is that in our case the preset linear combinations (2) of the AOs of the whole system are used in place of the AOs themselves. Correspondingly, this method is computationally expensive since all the one- and two-electron integrals which are necessary for the construction of the Fock matrix $`𝛀^A`$ are expressed as double and quadruple sums of the corresponding AO integrals.

## III Localized molecular orbitals for Si and MgO bulk crystals

In this section we shall examine LMOs obtained using three methods described above. Two crystalline systems with extreme types of chemical binding will be considered: MgO (ionic) and Si (covalent).

### III.1 General method

We are interested in calculating LMOs for perfect periodic solids. In this case the canonical MOs are Bloch-like solutions of the HF or Kohn-Sham equations. This means that the eigenvalue problem (9) is of an infinite dimension which make the calculation quite complicated; in addition, it would be necessary to use (and modify) a periodic electronic structure code. To avoid these difficulties, we suggest a very simple procedure based on a cluster method. The basic idea relies on the fact that when the cluster size is increased, the distribution of the electron density in its central region should become closer to the actual electron distribution of the infinite periodic system.

Our method is based on the following steps:

1. analyze the known electron charge density $`\rho (𝐫)`$ of the 3D periodic system in question to identify regions $`A`$, $`B`$, $`C`$, etc. which can each be associated with even numbers of localized electrons, e.g. atoms, ions, bonds; the density $`\rho (𝐫)`$ can be known from literature or own calculations;
2. consider a quantum cluster which contains region $`A`$ in its center (or close to it); terminate the cluster using pseudoatoms (see below) and/or an array of point charges to reproduce the correct Madelung field; obtain the occupied canonical orbitals $`\phi _i^c(𝐫)`$ for the whole cluster;
3. then consider a localization problem for region $`A`$ using one of the methods of the previous section; this should give the necessary number of LMOs $`\stackrel{~}{\phi }_a(𝐫)`$ ($`a=1,\mathrm{},n`$) as a linear combination of the occupied canonical MOs $`\left\{\phi _i^c(𝐫)\right\}`$, Eq. (4); note that in some cases when e.g. pseudoatoms are used to terminate the cluster, their contribution to the LMOs should be removed and the orbitals renormalized;
4. repeat the procedure 2-3 for larger clusters to ensure that the LMOs obtained have converged;
5. if other types of regions exist, repeat steps 2-4 for those regions as well; when finished, LMOs for the whole unit cell should be available; sometimes (as is the case for silicon), LMOs of some other regions can be obtained without additional calculation by simply translating and possibly rotating the LMOs of a single region;
6. the LMOs within the primitive cell can be displaced by all possible lattice translations to obtain the complete set of crystal LMOs spanning the complete occupied Fock space; these can now be employed for the calculation of the density $`\stackrel{~}{\rho }(𝐫)`$ according to Eq. (5); note that $`\stackrel{~}{\rho }(𝐫)`$ and the actual periodic density $`\rho (𝐫)`$, calculated using a periodic code and identical basis set, will not be exactly the same due to a number of approximations employed here in calculating the LMOs (see below).

Some general comments of our general method are necessary at this point. Firstly, the convergence of our procedure with the cluster size depends on the “boundary conditions” used in every case, i.e. it depends on the way the cluster is terminated; it is nothing but the embedding method itself. The latter is, however, unknown. Therefore, our procedure can be considered only as the first iteration and thus larger cluster sizes are expected. In principle, when the LMOs are obtained, one can use them for a new set of embedding calculations to obtain a better approximation for them, etc. In this case smaller cluster sizes may only be necessary.

Secondly, it has already been mentioned that any linear combination of the canonical set of MOs should lead to the same electron density. Therefore, one may think that the density $`\stackrel{~}{\rho }(𝐫)`$ obtained using the procedure outlined above will always result in the correct electron density $`\rho (𝐫)`$. This is, however, do not need to be the case due to a number of approximations adopted. Indeed, we only consider finite clusters with *ad-hoc* boundary conditions; in addition, the contribution of boundary cluster atoms may be modified when pseudoatoms are used. Finally, the cluster size may be insufficient to accommodate completely the LMOs. Therefore, the obtained occupied Fock space will never be exactly the same as that obtained using the periodic calculation. Hence, the comparison of the electron densities $`\stackrel{~}{\rho }(𝐫)`$ and $`\rho (𝐫)`$ may indicate on the quality of the calculated LMOs, and this method will be used in this paper.

The calculation of $`\stackrel{~}{\rho }(𝐫)`$ is performed by exploiting the periodic symmetry and representing the LMOs as an integral over the Brillouin zone. This method is exact and does not depend on the degree of localization of the LMOs. It also allows exact handling of the inverse of the overlap matrix in Eq. (5). Details of the method will be published elsewhere. To calculate the reference density $`\rho (𝐫)`$, we used our cluster calculations in the following way: (i) a parallelepiped in the central part of the cluster with the sides along the primitive lattice translations $`𝐚_1`$, $`𝐚_2`$ and $`𝐚_3`$ which is equivalent to the primitive unit cell is identified; its density is denoted $`\rho _c(𝐫)`$; (ii) the density of the whole 3D crystal is then modelled as $`\rho (𝐫)\rho _c(𝐫_1)`$, where $`𝐫_1`$ is obtained from $`𝐫`$ by removing any lattice translations (this is most conveniently performed by first calculating fractional coordinates of $`𝐫`$ in terms of $`𝐚_1`$, $`𝐚_2`$ and $`𝐚_3`$ and then removing integral parts from them). The larger size of the cluster is used in the calculations, the better approximation for the density $`\rho (𝐫)`$ will be obtained in this way.

All numerical calculations reported in this paper were done using the HF method and the Gamess-UK GAMESS-UK  code within the pseudopotential approximation.

### III.2 Localized orbitals for the MgO bulk

A sequence of three finite clusters of increased size Mg<sub>6</sub>O, Mg<sub>38</sub>O<sub>13</sub> and Mg<sub>44</sub>O<sub>19</sub> containing 7, 51 and 63 atoms, respectively, and surrounded by an array of nearly 10<sup>3</sup> point charges of $`\pm 2e`$ to simulate the Madelung field was considered. The largest cluster used is shown in Fig. 1 (*a*).

To consider explicitly only the valence electrons, for both Mg and O we used coreless Hartree-Fock pseudopotentials (CHF) with LP-31G basis set from Ref. CHF:Melius-Goddard74 . The number of electrons in each cluster was calculated by adopting a well known ionic character of the MgO crystal, i.e. by assuming that every Mg<sup>2+</sup> ion is associated with no electrons, while every O<sup>2-</sup> ion has eight electrons. The LP-31G basis set was used in all our calculations and the distance between the nearest Mg and O atoms is 2.122   irreg-MgO .

After the HF solution was obtained for every cluster, we applied the three localization procedures (methods M, P and E) considered in the previous section to obtain the LMOs for this system. Since there are only two atoms in the primitive cell, Mg and O, and it is well known that the valence electron density is localized predominantly on the O atoms, one can choose essentially a single region $`A`$ within the primitive unit cell to localize the LMOs into, namely, on the O atom. We should expect four LMOs localized on every O atom: one of the $`s`$ type, $`\stackrel{~}{\phi }_s(𝐫)`$, and three of the $`p`$ type, $`\stackrel{~}{\phi }_{p_x}(𝐫)`$, $`\stackrel{~}{\phi }_{p_y}(𝐫)`$ and $`\stackrel{~}{\phi }_{p_z}(𝐫)`$.

Therefore, when applying the methods M and P we used the $`s`$ and $`p_x`$, $`p_y`$ and $`p_z`$ type AOs centered on the O atom in the center of every cluster when applying Eqs. (12) and (13). In the case of method E, we considered the HF problem for a single oxygen ion O<sup>2-</sup> in the basis set of all occupied canonical MOs of the entire cluster. In every case exactly four LMOs were obtained as having the smallest eigenvalues; other states were found to be separated by a considerable gap.

The HF electron density through the central O atom for every cluster is shown in Fig. 2 (a).

One can see that the density is perfectly converged already for the smallest of the clusters. This means that the electron density in the center of any of the clusters can be considered as being very close to the density of the actual 3D periodic crystal calculated in the HF approximation using the same basis set.

The partial oxygen electron density

$$\rho _O(𝐫)=\stackrel{~}{\phi }_s^2(𝐫)+\stackrel{~}{\phi }_{p_x}^2(𝐫)+\stackrel{~}{\phi }_{p_y}^2(𝐫)+\stackrel{~}{\phi }_{p_z}^2(𝐫)$$
(15)
can be conveniently used to characterize the localization of the obtained LMOs $`\stackrel{~}{\phi }_s(𝐫)`$, $`\stackrel{~}{\phi }_{p_x}(𝐫)`$, etc. We compare $`\rho _O(𝐫)`$ obtained for all clusters using method M in Fig. 3 (a). Note that other methods give practically identical densities.

It is seen that all four LMOs are extremely well localized on the O atom (as one would expect for such an extremely ionic system) and converge very quickly with the cluster size. The LMOs are essentially identical for all three methods.

Obviously, LMOs associated with any other unit cell can now be obtained simply by moving the calculated four LMOs by the appropriate lattice translation. We have made a careful comparison of the total electron density $`\stackrel{~}{\rho }(𝐫)`$ constructed using Eq. (5) with the density $`\rho (𝐫)`$ calculated using the central part of the largest cluster. In particular, such a comparison is shown in Fig. 3 (b) along the (001) direction across the central O atom. One can see that either method results in the perfect matching between the density $`\stackrel{~}{\rho }(𝐫)`$ obtained using the LMOs (indistinguishable on the plot) and the reference density, $`\rho (𝐫)`$. Thus, all three localization techniques work equally well in the case of MgO and require clusters of very moderate sizes.

### III.3 Localized orbitals for the Si bulk

Crystalline Si has the diamond-type lattice with the distance between the nearest Si atoms of 2.35 . Each Si atom is surrounded by four neighbors forming four covalent bonds with them. There are two Si atoms and thus eight electrons (within the valence approximation) to be assigned to every primitive cell. Since each bond is associated with two electrons, there should be four bonds per cell. We expect that well localized orbitals can be constructed for this crystal if localization regions are associated with every two-electron bond. Therefore, in this case we have four regions $`A`$, $`B`$, $`C`$ and $`D`$ of identical nature in the primitive cell as shown schematically in Fig. 4.

Note that the choice of four inequivalent bonds is not unique. In our choice shown in the figure all four bonds share atom 1 and can be obtained from any single bond (e.g. the central bond between atoms 1 and 2) by applying appropriate displacements and rotations.

Three quantum clusters Si<sub>2</sub>Na<sub>6</sub>, Si<sub>8</sub>Na<sub>18</sub> and Si<sub>26</sub>Na<sub>18</sub>Mg<sub>12</sub>, containing 2, 8 and 26 Si atoms were used in our calculations; every cluster contains a single Si-Si pair in its center as shown in Fig. 1 (b). The core electrons of all Si atoms were described using the Hay-Wadt pseudopotential lanl . To terminate unsaturated bonds of the Si atoms located at the boundary of the clusters, we used Na- and Mg-like pseudoatoms which have the same pseudopotentials as the Si atoms. The Na pseudoatoms contribute a single electron to the cluster and are positioned at the correct Si-Si distance to saturate a single dangling bond. Mg-like atoms have two electrons and were used in the same way to saturate two dangling bonds from two nearest boundary Si atoms (see Fig. 1 (b)). Since the Si crystal is a highly covalent system, the Madelung field can be considered of a secondary importance and thus was neglected. The 66-21G basis set 66-21G  was used on the Si atoms in most cases. The basis set on pseudoatoms Na and Mg included only $`s`$ type AOs.

To construct the LMOs for the Si crystal, we oriented the coordinate system in such a way that the $`z`$ axis would pass along the central Si-Si bond of every cluster. This particular choice of the coordinate system is merely needed to simplify the choice of the AOs to be associated with the localization region $`A`$. Then, the HF solution was obtained which demonstrated a good degree of the $`sp^3`$ hybridization, as expected. When applying the localization methods M and P, AOs of the $`s`$ and $`p_z`$ types centered on the two central Si atoms were chosen as belonging to region $`A`$. In order to apply method E, the following SE was considered in place of the central Si-Si molecule: it consisted of two electrons and two pseudoatoms with the Si pseudopotential and the total charge $`+e`$ each. Effectively, this way the SE was chosen as a pseudo-hydrogen molecule with pseudo-hydrogen atoms at the Si-Si distance described each by the Si pseudopotential.

By analyzing the electron density in the central region of every cluster, we find that the largest cluster we considered is sufficient for our purposes. As an example, we show in Fig. 2 (b) the electron densities across the central Si-Si bond for the three clusters. The single LMOs calculated for the central Si-Si bond of the three clusters using methods M and E are shown in Fig. 5. The LMOs calculated using method P are very similar to those calculated using method M and thus are not shown here. The important conclusion which can be drawn from these pictures is that the LMO obtained using method M is practically converged with the cluster size. On the contrary, the LMO calculated using method E is not yet converged becoming more and more delocalized with the increase of the cluster size.

Even bigger clusters are thus needed to converge the LMO using this method. However, as was mentioned earlier in section II.3, larger clusters require a very expensive procedure of calculating two-electron integrals, and thus we did not consider larger systems. Note that another way of circumventing the convergence problem in method E might be to add some potential well to the HF problem for region $`A`$ to enforce a stronger localization; although we did not pursue this idea in this work, we may consider it in the future.

To construct the electron density $`\stackrel{~}{\rho }(𝐫)`$ of the whole Si crystal we need other three LMOs assigned to the same primitive cell. These are obtained by appropriate rotations and displacements of the central bond LMO considered above. The LMOs corresponding to other crystal cells are then obtained by applying appropriate lattice translations. The electron densities $`\stackrel{~}{\rho }(𝐫)`$ obtained using the three localization methods for the largest cluster have been thoroughly compared with the reference density, $`\rho (𝐫)`$, obtained by translating the central part of the same cluster. It has been found that comparison of the densities along the central Si-Si bond, as shown in Fig. 6 (a), reflects well the extent in which the densities match with each other.

One can see that both methods M and P lead to the electron density of the Si crystal which is very close to the reference density $`\rho (𝐫)`$. However, the density $`\stackrel{~}{\rho }(𝐫)`$ calculated from the LMOs obtained using method E shows an unphysical oscillatory behavior which is due to their poor localization and thus an insufficient cluster size used to construct them.

## IV Discussion and conclusions

We have seen in the previous section that for an ionic system such as MgO all localization procedures give identical results and do not require large cluster sizes. This is because a natural localization takes place in those systems. In fact we find that practically identical LMOs can be obtained for MgO using methods M and P if in addition to the AOs centered on the central O atom one also adds AOs of any of the nearest Mg atoms to define region $`A`$. This means that these localization procedures are sufficiently flexible in terms of the AOs used to define the localization regions.

We find that the localization procedures for covalent system with strong hybridization in its chemical bonding are more sensitive to the choice of the localization region and the particular localization method. We have seen above that method E fails for this system since the LMOs it produces are not sufficiently localized within the cluster sizes we use.

We have also experimented with localization methods M and P by trying to use different definitions of region $`A`$ in order to find the LMO. First of all, we used all AOs on the two central Si atoms. Similarly to the case of MgO, in this case an LMO practically identical to the one which we calculated using only $`s`$ and $`p_z`$ types AOs was obtained. The $`\stackrel{~}{\rho }`$ densities calculated in these two cases are also the same as is demonstrated in Fig. 6 (b) by the good matching between $`\stackrel{~}{\rho }\stackrel{~}{\rho }_{sxyz}`$ and $`\stackrel{~}{\rho }\stackrel{~}{\rho }_{sz}`$. The two methods succeeded since the chemical bonding is essentially correctly reproduced by either choice of region $`A`$. This is confirmed by the contour plot of the partial density $`\stackrel{~}{\phi }^2(𝐫)`$ associated with the LMO and calculated using all AOs centered on the two central Si atoms: as shown in Fig. 7, the LMO essentially corresponds to the Si-Si bond.

In contrast, when assuming a wrong character of the chemical bonding in Si, we obtained LMOs which were either not very well localized or had completely unexpected (unphysical) spatial distribution. For instance, assuming atomic character of the chemical bonding, we attempted to use AOs on a *single* Si atom to define the single localization region in the primitive cell. This assumption gave four LMOs similarly to the MgO case. We find, however, that these LMOs become much less localized and, as a result, the constructed electron density $`\stackrel{~}{\rho }\stackrel{~}{\rho }_1`$ is very different from $`\stackrel{~}{\rho }_{sz}`$ as is obvious from Fig. 6 (b). Another example of an “unwise” choice of the localization region is to use only $`s`$, $`p_x`$ and $`p_y`$ AOs of the two central Si atoms to define region $`A`$ (recall that the two Si atoms are positioned along the $`z`$ axis). In this case the total electron density $`\stackrel{~}{\rho }\stackrel{~}{\rho }_{sxy}`$ also shown in Fig. 6 (b) somewhat differs from $`\stackrel{~}{\rho }_{sz}`$, but, at the same time, reproduces all its main features. However, the LMO itself appears to have a completely different spatial distribution (see Fig. 7): it is not anymore localized on the Si-Si bond, but instead was found to be delocalized over a large volume around it. This explains why the total density $`\stackrel{~}{\rho }_{sxy}`$ was found quite different from the reference one: much larger cluster should be considered to accommodate fully the LMO obtained using this particular choice of region $`A`$. Therefore, the proper choice of the localization regions which reflect the chemistry of the given crystal results in more localized orbitals and thus much smaller cluster sizes needed to construct them.

Concluding, a simple method based on a cluster approach was suggested in order to construct localized molecular orbitals (LMOs) for a periodic solid. Our method does not require usage of periodic codes, is thus much easier to implement in practice and also, in addition, can also be applied to nonperiodic systems. The work in this direction is presently in progress and will be published elsewhere.

Several localization procedures were analyzed and two crystals were considered in detail (MgO and Si) which corresponds to the two cases of extreme types of chemical bonding - ionic and covalent. We find that two localization procedures considered, one based on the Mulliken populations (method M) and another on a projection operator (method P), give well localized orbitals with the expected conventional meaning adopted in chemistry using already quite moderate cluster sizes. The third procedure, based on the minimization of the HF energy of a structural element (one or two atoms) demonstrated a much slower convergence with the cluster size and is found to be also computationally expensive.

Two cases considered here, MgO and Si crystals, have a well known type of chemical bonding and thus the choice of the localization regions in these two cases was obvious. At the same time, we find that there is a certain degree of flexibility in choosing the localization regions and this can be exploited in the cases of more complicated (e.g. intermediate) types of chemical bonding. This work is being done in our laboratory at present and will be a matter for future publications.

## Acknowledgments

We are extremely grateful to I. V. Abarenkov and A. Shluger for a number of useful and stimulating discussions during his stay in London. O.D. would also like to acknowledge the financial support from the Leverhulme Trust (grant F/07134/S) which made this work possible.


format: text

text:
# The three-state layered neural network with finite dilution

## 1 Introduction

The statics and the dynamics of large attractor and feed-forward neural networks, inspired in features of biological networks, have been studied over some time. Some of these features are the asymmetry and the finite dilution of the synaptic connections as well as the low activity of the patterns and the neurons. It has been suggested that a recurrent attractor neural network model of multi-state neurons with these characteristics may describe the short term memorization performance of the $`CA_3`$ region of the hyppocampus in both the human brain and in the brain of primates, and results of numerical simulations on the retrieval behavior are now available . Full results on the dynamic evolution of the retrieval overlap and on the phase diagrams for the stationary states of low-activity networks of multi-state neurons with asymmetric interactions and finite dilution are still missing. It would be desirable to have such results in order to describe other implementations of neural network models, among them as devices to account for long-term memory in the brain .

The retrieval behavior and thermodynamic properties of symmetrically diluted Q-Ising recurrent neural networks with finite dilution and low activity patterns have been studied and full phase diagrams have been obtained for $`Q=3`$ and $`Q=4`$ states as well as for a network with continuous response neurons . In the case of a discrete number of states, phase boundaries between states of low and high performance are found to disappear beyond a finite dilution, allowing for the biologically appealing feature of a continuous improvement of the behavior of the network without the need of a precise threshold adjustment. The full study of the dynamics of these networks is rather involved (see ref. ) and for practical, either hardware or biological implementations of neural networks, it would be interesting to have a simple dynamical system.

A suitably tractable model to study these issues is a feed-forward layered network with no feedback loops, which has an exactly solvable non-trivial dynamics, in which the only non-zero synaptic connections are the asymmetric interactions that pass information from one layer to the next. Despite the fact that the connectivity of the layered network is much lower than that of a recurrent network, the presence of non-trivial correlations makes it an ideal model to test the qualitative behavior of feature dependence in recurrent networks. The study of the model in itself is also of interest in view of the practical applications of feed-forward layered networks.

The purpose of the present work is to present new results on the retrieval performance, the information content and the dynamic evolution for this network with three-state Ising neurons and finite dilution. Results for the diluted layered network with binary units and for the three-state network with no dilution can be found in the literature .

The outline of the paper is the following. In Section 2 we introduce the model and the relevant macroscopic variables. In Section 3 we derive the recursion relations for these variables that establish the dynamics for the model. We discuss the results for the phase diagrams of the stationary states and the dynamical evolution of the retrieval overlaps in Section 4, and end with concluding remarks in Section 5.

## 2 The Model

The network model consists of $`L`$ layers with $`N`$ neurons on each layer, that can take values $`\sigma _i^l`$, $`l=1,\mathrm{},L;i=1,\mathrm{},N`$ from the set $`S\{1,0,+1\}`$, where $`\pm 1`$ denote the active states. A macroscopic number of $`p=\alpha cN`$ ternary patterns, where $`c`$ is the probability that two neurons are connected, is taken from a set of independent identically distributed random variables $`\{\xi _i^{\mu ,l}=0,\pm 1\}`$, $`\mu =1,\mathrm{},p`$, with the probability distribution

$$\text{Prob}(\xi _i^{\mu ,l})=a\delta (|\xi _i^{\mu ,l}|^21)+(1a)\delta (\xi _i^{\mu ,l})$$
(1)
which is assumed to be the same for every layer and where $`\pm 1`$ are the active patterns. The mean of each pattern (over-lined) is zero and $`a=\overline{(\xi _i^{\mu ,l})^2}`$, denotes their activity. A new random set of $`p`$ patterns is generated on each layer and the whole set on two consecutive layers is embedded in the diluted network according to the generalized Hebbian learning rule

$$J_{ij}^l=\frac{c_{ij}^l}{acN}\underset{\mu =1}{\overset{p}{}}\xi _i^{\mu ,l+1}\xi _j^{\mu ,l},$$
(2)
where $`\{c_{ij}^l\}`$ is a set of independent identically distributed random variables that account for the dilution of the synapses, particularly when the patterns are active, and such that $`c_{ij}^l=1`$ with probability $`c`$ and zero with probability $`1c`$, for all $`l`$. Thus, $`cN`$ is the mean number of neurons connected to each neuron. For the fully connected layered network, $`c=1`$, and the case of extreme dilution corresponds to the limit $`c0`$, after taking the macroscopic $`N\mathrm{}`$ limit. The dilution introduces an additional randomness into the dynamics of the fully connected network in the form of a static noise $`z_{ij}^l/\sqrt{N}`$ with mean zero and variance $`\alpha (1c)/N`$ . There is no contribution to the learning rule from patterns in the same layer.

The states of the network change as follows. Given a configuration on the first layer, $`𝝈_N^1\{\sigma _j^1\},j=1,\mathrm{},N`$, the state $`\sigma _i^{l+1}`$ of unit $`i`$ on layer $`l+1`$ is determined exclusively by the configuration $`𝝈_N^l`$ of the units in the previous layer according to the stochastic law

$$\text{Prob}(\sigma _i^{l+1}=sS|𝝈_N^l)=\frac{\mathrm{exp}[\beta ϵ_i(s|𝝈_N^l)]}{_{sS}\mathrm{exp}[\beta ϵ_i(s|𝝈_N^l)]},$$
(3)
in terms of the single-site energy function on that unit

$$ϵ_i(s|𝝈_N^l)=sh_i^{l+1}(𝝈_N^l)+\theta _i^{l+1}s^2,$$
(4)
where

$$h_i^{l+1}(𝝈_N^l)=\underset{j=1}{\overset{N}{}}J_{ij}^l\sigma _j^l,$$
(5)
is the acting local field and $`\theta _i^{l+1}`$ is a local externally adjustable threshold parameter on layer $`l+1`$. Since the only changes in the network, in both the synaptic connections and the states of the units, are between units in consecutive layers one may associate the layer index with a discrete time step $`t`$ and we do this in the following. Thus, the evolution of the network proceeds according to a parallel dynamics in which the states of all neurons are updated at each time step.

Next, we consider the relevant quantities that describe the performance of the network. First, the retrieval overlap between the state of the network and the pattern $`\{\xi _i^{\mu ,t}\}`$ at time (layer) $`t`$, given by

$$\stackrel{~}{m}_i^{\mu ,t}=\frac{1}{acN}\underset{j}{}c_{ij}^t\xi _j^{\mu ,t}\sigma _j^t,$$
(6)
which depends on the site $`i`$. The dynamical activity and the activity overlap are defined as

$$q_0^t=\frac{1}{N}\underset{i}{}(\sigma _i^t)^2,n_i^{\mu ,t}=\frac{1}{acN}\underset{i}{}(c_{ij}^t\xi _i^{\mu ,t}\sigma _i^t)^2,$$
(7)
respectively. The latter will only be needed for the condensed pattern, that is the stored pattern to be retrieved.

For the mutual information, we need the conditional probability distribution $`\text{Prob}(\sigma _i^t|\xi _i^{\mu ,t})`$ that a neuron $`i`$ is in the state $`\sigma _i^t`$ on layer $`t`$ given that the $`i`$-th bit of the condensed pattern is $`\xi _i^{\mu ,t}`$. As a consequence of the independence of the states of the units on a given layer, it is sufficient to consider the distribution for a single typical neuron, and we omit here the index $`i`$. We also omit here, for clarity, the time index and use

$$\text{Prob}(\sigma |\xi ^\mu )=(s_\xi +m^\mu \xi ^\mu \sigma )\delta (\sigma ^21)+(1s_\xi )\delta (\sigma ),$$
(8)
where

$$s_\xi =s^\mu +l^\mu (\xi ^\mu )^2,s^\mu =\frac{q_0an^\mu }{1a},l^\mu =\frac{n^\mu q_0}{1a}.$$
(9)
The mutual information between patterns and neurons, regarding the patterns as the inputs and the neuron states as the output of the network channel on each layer, is an architecture independent property given by

$$I^\mu (\sigma ,\xi ^\mu )=S(\sigma )\overline{S(\sigma |\xi ^\mu )},$$
(10)
where

$$S(\sigma )=q_0\mathrm{ln}(q_0/2)(1q_0)\mathrm{ln}(1q_0)$$
(11)
is the entropy and $`\overline{S(\sigma |\xi ^\mu )}=aS_a+(1a)S_{1a}`$ is the equivocation term with

$`S_a`$ $`=`$ $`c_+^\mu \mathrm{ln}c_+^\mu c_{}^\mu \mathrm{ln}c_{}^\mu (1n^\mu )\mathrm{ln}(1n^\mu )`$
$`S_{1a}`$ $`=`$ $`s^\mu \mathrm{ln}(s^\mu /2)(1s^\mu )\mathrm{ln}(1s^\mu ).`$ (12)
Here, $`c_\pm ^\mu =(n^\mu \pm m^\mu )/2`$ and $`s^\mu `$ is the parameter in the conditional probability $`\text{Prob}(\sigma |\xi ^\mu )`$. The mutual information can then be used to obtain the information content of the network, $`i^\mu =I^\mu \alpha `$, where $`\alpha =p/cN`$ is the storage ratio.

## 3 Recursion Relations

We consider in this work the retrieval of a single (condensed) pattern, say $`\xi _i^{1,t}`$, in the dynamics of the network with a finite overlap $`\stackrel{~}{m}_i^{1,t}=O(1)`$ and the remaining $`\mu >1`$ overlaps $`\stackrel{~}{m}_i^{\mu ,t}=O(1/\sqrt{N})`$. The interest is in the mean overlap, $`m^t=[\overline{\stackrel{~}{m}_i^{1,t}}]`$, where $`[\mathrm{}]`$ denotes the average over the $`c_{ij}`$, $`<\mathrm{}>`$ denotes the thermal average with Eq.(3) and the bar means the average over the patterns. We also need the activity overlap with pattern $`\xi _i^{1,t}`$, and take $`n^t=[\overline{n_i^{1,t}}]`$.

The recurrence relations that describe the dynamics of the network for large $`N`$ follow from the local field, Eq.(5), written as

$$h_i^{t+1}=\xi _i^{1,t+1}m^t+z\mathrm{\Delta }^t,$$
(13)
where the average condensed overlap $`m^t=\overline{\sigma ^t\xi ^{1,t}/a}`$ depends on the local field at time $`t1`$ and $`z`$ is a Gaussian random variable with zero mean and unit variance that comes from the action of the macroscopic number of random overlaps of the diluted network with the uncondensed patterns and the use of the central limit theorem. The layer-dependent variance of the local field becomes site independent and is given by

$$(\mathrm{\Delta }^t)^2=a\underset{\mu >1}{}\overline{[(\stackrel{~}{m}_\mu ^t)^2]}.$$
(14)
A direct calculation in the large-$`N`$ limit yields

$$(\mathrm{\Delta }^t)^2=\alpha (1c)q_0^t+(\mathrm{\Delta }_0^t)^2,$$
(15)
where $`\alpha =p/cN`$, $`q_0^t=\overline{(\sigma ^t)^2_{\sigma |\xi }}`$ is the dynamical activity and

$$(\mathrm{\Delta }_0^t)^2=a\underset{\mu >1}{}\overline{(m_\mu ^t)^2}$$
(16)
is now the variance of the Gaussian noise for the connected layered network in terms of the overlap with the uncondensed patterns

$$m_\mu ^t=\frac{1}{aN}\underset{i}{}\xi _i^{\mu ,t}\sigma _i^t.$$
(17)
Thus, the noise in the local field is a superposition of two Gaussian noises, one due to the dilution of the synaptic connections and the other to the macroscopic number of uncondensed patterns, in extension of an earlier result for the binary network .

Since our main interest in this work is in the effects of synaptic dilution in a low-activity network, we take a uniform and time-independent threshold $`\theta _i^t=\theta `$. The averages $`\sigma ^t`$ and $`(\sigma ^t)^2`$ are then given, respectively, by

$$F_\beta (h^t,\theta )=\frac{\mathrm{sinh}(\beta h^t)}{\frac{1}{2}e^{\beta \theta }+\mathrm{cosh}(\beta h^t)},G_\beta (h^t,\theta )=\frac{\mathrm{cosh}(\beta h^t)}{\frac{1}{2}e^{\beta \theta }+\mathrm{cosh}(\beta h^t)}$$
(18)
which, in the zero temperature limit, $`\beta \mathrm{}`$, become

$$F_{\mathrm{}}=\text{sign}(h^t)\mathrm{\Theta }(|h^t|\theta ),G_{\mathrm{}}=\mathrm{\Theta }(|h^t|\theta ),$$
(19)
where $`\mathrm{\Theta }(x)`$ is the usual step function. The performance with a self-adjusting time-dependent threshold has been considered in a fully connected layered network  and in other architectures \[10-14\]. Exact dynamic equations are then obtained in the large $`N`$ limit in the form of recursion relations for $`m^t`$, $`q_0^t`$ and $`n^t`$, where the latter two are needed for the information content of the network. Similarly, an exact recursion relation for the second term in the width of the stochastic noise, Eq.(15), is obtained in the form

$$(\mathrm{\Delta }_0^{t+1})^2=\alpha q_0^t+(\chi ^t)^2(\mathrm{\Delta }_0^t)^2,$$
(20)
where

$$\chi ^t=\beta (q_0^tq_1^t)$$
(21)
is the susceptibility in which $`q_1^t=\overline{\sigma ^t^2}`$. Thus we obtain

$`m^{t+1}`$ $`=`$ $`{\displaystyle DzF_\beta (m^t+z\mathrm{\Delta }^t,\theta )}`$ (22)
$`q_0^{t+1}`$ $`=`$ $`{\displaystyle Dz\{aG_\beta (m^t+z\mathrm{\Delta }^t,\theta )+(1a)G_\beta (z\mathrm{\Delta }^t,\theta )\}}`$ (23)
$`q_1^{t+1}`$ $`=`$ $`{\displaystyle Dz\{aF_\beta ^2(m^t+z\mathrm{\Delta }^t,\theta )+(1a)F_\beta ^2(z\mathrm{\Delta }^t,\theta )\}}`$ (24)
where, as usual, $`Dz=\mathrm{exp}(z^2/2)dz/\sqrt{2\pi }`$. We also have

$$n^{t+1}=DzG_\beta (m^t+z\mathrm{\Delta }^t,\theta ).$$
(25)
The dynamics, including transients, and the stationary states of the diluted layered network follows then from the solutions of these equations together with Eq.(20). The stationary states are reached when $`m^{t+1}=m^t`$, $`q_0^{t+1}=q_0^t`$, $`q_1^{t+1}=q_1^t`$, $`n^{t+1}=n^t`$ and $`\mathrm{\Delta }_0^{t+1}=\mathrm{\Delta }_0^t`$, and we call the first three, respectively, $`m`$, $`q_0`$ and $`q_1`$. The phase diagrams for the stationary states and the time evolution of the order parameters are discussed in the next section.

## 4 Dynamics and stationary states

The stable stationary states yield one or more retrieval phases $`R(m>0,q_1>0)`$ and a spin-glass phase $`SG(m=0,q_1>0)`$, all as sustained activity solutions with $`q_0>0`$. Since the model is a dynamical one, the stability criterium is that the change in the order parameters should become smaller in every one of the final steps of the iteration procedure of the flow towards an attractor fixed point. Thus, for the retrieval overlap one has

$$\underset{\delta m(t)0}{lim}\frac{\delta m(t+1)}{\delta m(t)}<1,$$
(26)
and similar relationships for the other parameters.

We are interested in this paper in the characteristic features of finite dilution of the phase diagrams and in the specific performance of the network. Different kinds of phase diagrams are obtained depending on the pattern activity $`a`$ and on $`c`$. In the case of full connectivity ($`c=1`$) and low $`a`$, we find the $`(\alpha ,\theta )`$ phase diagram shown in Fig. 1 for $`a=0.5`$ and either $`T=0`$ (full lines) or $`T=0.05`$ (dashed lines). The lines represent discontinuous phase boundaries where the locally stable retrieval states appear with decreasing load $`\alpha `$ below a critical $`\alpha _c`$ and the dotted line indicates the locus of optimal performance which yields the largest load capacity that still sustains retrieval behavior. There is a weak retrieval phase (I) and a strong retrieval phase (II) separated by a discontinuous phase boundary. There is also a smaller retrieval phase in the lower triangular region. For larger activity and full connectivity, as in the case of uniform patterns ($`a=2/3`$), there is mostly a single retrieval phase with a strong $`\alpha `$ dependent optimal performance and this feature remains even for finite dilution, say for $`c=0.5`$.

A different situation appears for both finite dilution, with $`c`$ smaller than a critical activity-dependent $`c^{}`$, and for low activity where part of the phase boundary between regions I and II disappears. For $`a=0.5`$ this is the case when $`c^{}=0.87`$ at $`\alpha =0.34`$ and $`T=0`$. As shown in Fig. 2, for typical $`c=0.8`$ and $`a=0.5`$, there is now a gradual transition from region I to region II and, as $`c`$ decreases, the continuous transition region increases until a stage is reached where the maxima in both regions merge into a single maximum with a lower optimal threshold $`\theta `$. The phase diagram in Fig. 2 is reminiscent of the diagram for the finite diluted recurrent $`Q=3`$ state network with uniform patterns .

In order to study the dynamics and the stability of the phases, we consider the time evolution of the retrieval overlap shown in Fig. 3 for $`c=0.8`$, $`a=0.5`$, $`\alpha =0.06`$, $`T=0`$ and $`\theta `$ between $`0.293`$ and $`0.33`$. For small $`\theta `$ the retrieval state is not stable as shown by the typical lower curve ($`0.293`$), and there is a whole part of region I where this is the case. Since it is not clear that these instabilities have any significance, we do not explore this issue further in this work. On the other hand, as $`\theta `$ increases, the retrieval state becomes stable and a maximum overlap is reached after a relatively short transient.

A further property of the network is the increase in the maximum information content with dilution as shown near to the optimal performance with $`\theta =0.5`$ for $`a=0.5`$ and various values of $`c`$ in Fig. 4. Note that the result of a non-zero information content for increasing $`\alpha `$ with dilution is not surprising since we defined $`\alpha =p/cN`$, but the increase of the maximum information content with dilution is a novel feature.

There is also a considerable increase in the maximum information content with an increasing threshold in the good performance region II, as shown in Fig. 5 for $`c=0.8`$, $`a=0.5`$ and various values of $`\theta `$. Also, the decrease from the maximum is smoother with dilution than in the case of the fully connected layered network. Finally, one may consider the robustness of the network to synaptic noise and in Fig. 6 we show the maximum information content, $`i_{max}`$ for $`a=0.5`$, $`\theta =0.5`$ and various values of $`c`$. Clearly, there is a more gradual decrease in performance with $`T`$ for lower connectivity.

## 5 Summary and conclusions

We derived the exact recursion relations in the large-$`N`$ limit that describe the time evolution and the phase diagrams for the stationary states of the macroscopic variables in a three-state layered feed-forward network with finite synaptic dilution. Synaptic dilution appears as a static stochastic noise for a macroscopic number of stored patterns . This is a model with asymmetric interactions between units in consecutive layers and we allow for variable pattern activity in training the network with ternary patterns. Instead of discontinuous phase boundaries between retrieval and non-retrieval states, or between qualitatively different retrieval states for full synaptic connectivity, we find that there can be a continuous change with no boundary at all from weak to optimal retrieval states with a varying threshold in the local field, for low pattern activity. Phase boundaries of continuous transitions are known to appear in extremely diluted symmetric or asymmetric networks, but we emphasize that the continuous changeover in the present model is due to the joint action of finite dilution and low activity patterns.

In view of a similar recent result for a three-state recurrent diluted network with symmetric synaptic connections , one may conclude that this is a feature of finite dilution which is independent of both the network architecture and of the interaction symmetry. Thus, provided there is an above minimum threshold such that a network attains the ability to retrieve a nominated pattern after eliminating undesirable transient states, the good retrieval performance does not depend in an essential way on a precise threshold adjustment and this may explain why biological networks, in which there are no precise thresholds, can have a good performance despite a fraction of missing synaptic connections. This is an activity dependent property and it might help in the study of plasticity in neural networks. It may also be useful for artificial neural networks.

It may be worthwhile to note the asymmetric dual role of the threshold dependence discussed in this work. Whereas there is a continuous improvement in network performance for low to moderate threshold, there is an abrupt end to the performance for large threshold, as one would expect, since in the latter case mostly the inactive states of the network become dominant. We also showed that there is an improvement with finite dilution in both the size of the information content transmitted by the network and in the continuous changes of the information with a varying threshold. That is, as long as there is a convergence to stable stationary states, here again the network performance seems to be less sensitive to threshold adjustment.

The network behavior discussed in this work is also expected to appear in a diluted layered $`Q=4`$ state network with low activity patterns, based on recent results on a recurrent network which exhibits a continuous improvement in network performance with varying low-to-moderate threshold for low activity patterns .

## 6 Acknowledgments

The work of one of the authors (WKT) was financially supported, in part, by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico), Brazil.


text:
# 1 Introduction

## 1 Introduction

Little is currently understood about how magnetism arises at the nanoscale and how it affects the properties of nanosized metallic objects. Systems of special interest in this context are nanowires and atomic-sized nanocontacts. The one-dimensional (1D) character of such systems causes specific physical phenomena to appear, most notably quantized ballistic conductance.

These phenomena will be heavily affected by the possible presence of magnetism in the nanosystem, especially of a genuine Hund’s rule magnetic order parameter. Here, we report theoretical studies of magnetism in Pt monowires, i.e., wires consisting of a single line of equally spaced Pt atoms. Such objects, several atoms long, have recently been observed by Rodrigues et al.

Pt is a transition metal, with a partially filled $`5d`$ shell. The free Pt atom has a $`d^9s^1`$ configuration, giving a total magnetic moment of 2 $`\mu _B`$, and therefore it is reasonable to expect that a sufficiently stretched wire should eventually show some type of magnetic ordering. In particular, at the localized side of a Mott transition, one could expect an antiferromagnetic ordering. In the bulk metal, the Pt $`5d`$ band is too wide to provoke spin-polarization and there is no evidence of Pt showing surface magnetism either. Yet, the $`5d`$ band is still only partially occupied, and the density of states at the Fermi level is quite high. It appears therefore to be an open question whether Pt nanowires might exhibit ferromagnetic ordering. If that were the case, it might result in interesting phenomena such as spin-polarized current flow and magnetic-field dependent conductance, phenomena that could prove highly useful in spintronics applications.

Of course, thermal fluctuations, very large in a nanosystem, will generally act to destroy static magnetic order in the absence of an external field. Sufficiently slow fluctuations transform a nanomagnet to a superparamagnetic state, where magnetization fluctuates on some time scale, between equivalent magnetic valleys, separated, e.g., by anisotropy-induced energy barriers. If the barriers are sufficiently large — and the fluctuations sufficiently slow — the nanosystem spends most of the time in a single magnetic valley, and will for many practical purposes behave as magnetic. We may in these circumstances be allowed to neglect fluctuations altogether, and to approximate the calculated properties of the superparamagnetic nanosystem with those of a statically magnetized one. Experimentally, evidence of 1D superparamagnetism with fluctuations sufficiently slow on the time scale of the probe has been reported in Co atomic chains deposited at Pt surface steps.

## 2 Method

The density-functional calculations reported here are all-electron, in order to rule out possible sources of doubt that might arise when using pseudopotentials in presence of magnetism and in a nonstandard configurations. We employed an all-electron full-potential linear muffin-tin orbital (FP-LMTO) basis set together with a generalized gradient approximation (GGA) to the exchange-correlation functional. As a double check, some of the calculations were repeated using the linear augmented plane-wave (LAPW) code WIEN97. None of these calculations assume any shape approximation of the potential or wave functions.

We performed both scalar relativistic (SR) calculations, and calculations including the spin-orbit coupling as well as the scalar-relativistic terms. The latter will be referred to as “fully relativistic” (FR) calculations in the following, although we are not strictly solving the full Dirac equation, or making use of current density functional theory. In the FR calculations, the spin axis was chosen to be aligned along the wire direction.

The calculations were performed with inherently three-dimensional codes. Thus, the infinitely long, straight, isolated monatomic nanowire was simulated by a regular array of well-spaced nanowires. Convergence of the magnetic moment was checked with respect to k-point mesh density, Fourier mesh density, tail energies, and wire-wire vacuum distance.

## 3 Results and Discussion

The magnetic spin moment per atom monowire as a function of bond length is shown in Fig. 1. The solid line refers to the fully relativistic calculation (FR), and the dotted line to the scalar relativistic (SR) calculation. As seen, the magnetic profiles for the SR and FR calculations are very different. The SR calculation for Pt predicts this metal to be magnetic only for wire bond lengths larger than around 2.7 Å, which corresponds to a rather stretched wire, whereas the FR calculation predicts it to be magnetic in the whole range of bond lengths plotted (2.2 Å to 3.2 Å), with a moment of 0.6 $`\mu _B`$ per atom at the equilibrium bond length 2.48 Å. For sufficiently large bond lengths (in the interval shown in the figure) the magnetic moment reaches a plateau value, still well below the atomic spin moment 2 $`\mu _B`$. For even larger bond lengths (not shown), the magnetic moment eventually approaches the atomic spin moment 2 $`\mu _B`$. The energy gain per atom due to spin polarization is rather small, around 8 meV for the FR calculation at the equilibrium bond length 2.48 Å. Antiferromagnetic Pt monowire configurations were also tested for bond lengths around the equilibrium one, and found to be energetically unstable compared to the ferromagnetic configuration.

Several mechanisms, some favoring and other disfavoring a spin-polarized ground state, are at work in the wire. The number of nearest neighbors is only two in the wire, compared to 12 in bulk. This reduction of the number of nearest neighbors in the wire compared to the bulk causes a narrowing of the $`5d`$ band, and the band width may become sufficiently small that the gain in exchange energy due to spin polarization is larger than the increase in kinetic energy. On the other hand, the equilibrium bond length is significantly smaller in the wire, 10% smaller than in bulk, which partly counteracts the band-narrowing effect of the reduced number of neighbors.

Highly important in this context are the very sharp van Hove singularities caused by the one-dimensionality of the system. These give rise to a very high density of states at the Fermi level when a van Hove singularity is sufficiently close to the Fermi level, which in turn results in a Stoner product larger than one, and thus spontaneous spin-polarization.

In order to provide a more detailed analysis of the origin of the magnetism in the wires, we found it useful to analyze the wire band structures. Fig. 2 shows FR band structures of the Pt monowire for several different bond lengths. The bands run from the zone center, $`\mathrm{\Gamma }`$, to the zone edge, A, in the direction of the wire. Fig. 3 shows the corresponding SR band structures. The sharp van Hove singularities manifest themselves in the band structures through the horizontal band edges at the zone edge A and zone center $`\mathrm{\Gamma }`$. Since the orbital character is of critical importance for the moment formation, it is useful to divide the band structure into distinct orbitals — illustrated in Fig. 4 for the FR case. This figure has four panels, displaying separately the $`s`$, $`d_z`$, $`(d_{xz},d_{yz})`$ and $`(d_{xy},d_{x^2y^2})`$ characters of the bands. The vertical error bars, or “thickness”, of the bands indicate the relative character weight. In the SR case (not shown), the $`(d_{xz},d_{yz})`$ and $`(d_{xy},d_{x^2y^2})`$ orbitals correspond each to separate bands, and the $`s`$ and $`d_z`$ orbitals hybridize, forming two bands of high dispersion. With the help of Fig. 4, we can immediately recognize an important mechanism that favors spin-polarization in the wire. As seen, the bands have mostly $`d`$ character at the edges, and therefore the exchange energy gain will be rather large if a band spin-splits so that one of the spin-channel band edges ends up above the Fermi level, and the other one below. 111Strictly speaking, the spin-orbit coupling will mix the two spin channels so that, in general, an eigenvalue will have both majority and minority spin character. However, we found that this mixing is so small, typically just a few percent, that it is irrelevant for the qualitative discussion we make here. Thus, if a band edge ends up sufficiently near the Fermi level, we may expect a magnetic moment to develop. While apparently similar to the magnetization of the jellium wire, magnetism here is much more substantial, since the $`d`$ states involve a much stronger Hund’s rule exchange. We say that the magnetism exhibited in the Pt wire is Hund’s rule magnetism, in order to differentiate it from the situation in the jellium wire.

By comparing the magnetic and nonmagnetic band structures, we see that the relatively flat $`(d_{xy},d_{x^2y^2})`$ bands play the leading role in the formation of the magnetic state. In the SR band structures, these bands (all degenerate) sit well below the Fermi level, between $`1`$ eV and $`0.2`$ eV for the equilibrium bond length, and thus cannot contribute to spin polarization. For larger bond lengths than around 2.6 Å, the edge at A moves critically close to the Fermi level, and the bands split. In the FR calculations, however, a band edge at A of this symmetry is pinned to the Fermi level in the whole range of studied bond lengths, down to 2.2 Å, explaining the large difference in magnetic profile between the SR and FR calculations. This band edge is close to the Fermi level in the FR calculation but not in the SR calculation simply because in the FR calculation, the $`(d_{xy},d_{x^2y^2})`$ orbital hybridizes partly with the $`(d_{xz},d_{yz})`$ orbital, with the result that the $`(d_{xy},d_{x^2y^2})`$ band splits into an upper part of purely $`(d_{xy},d_{x^2y^2})`$ character, and a lower part hybridizing strongly with the $`(d_{xy},d_{x^2y^2})`$ orbital, thereby splitting up into two bands. If the average position of all the $`(d_{xy},d_{x^2y^2})`$ levels is taken, we end up more or less where the SR $`(d_{xy},d_{x^2y^2})`$ orbital originally was, energy-wise. Thus, one net effect of the spin-orbit coupling is that it increases the energy of part of the $`(d_{xy},d_{x^2y^2})`$ orbital. This also explains why the magnetic moment of the FR calculation is smaller than in the SR calculation in the bond-length range 2.7 Å to 3.1 Å.

Three other band edges, a $`d_{z^2}`$-dominated one located at A and the $`(d_{xz},d_{yz})`$-dominated ones with located at $`\mathrm{\Gamma }`$ are also important, and add to the size of the magnetic moment as they split around the Fermi level once the spin-polarization has been triggered.

The magnetic or superparamagnetic state of a small piece of nanowire bridging between nonmagnetic tips might at first seem very problematic to detect. It should in fact be detectable by measuring ballistic conductance as a function of both temperature and of an external magnetic field. The field, even a modest one depending on temperature, can drive a superparamagnetic wire from thermally disordered to fully spin polarized. In this polarized state, the nature and number of current-carrying channels, each corresponding to a band crossing the Fermi level, will differ from that of the nonmagnetic state. The channel number characteristic of the magnetic state should in fact last well into the superparamagnetic state at zero field. In the case of our infinitely long Pt monowire, we find that the nonmagnetic wire has 10 open conductance channels and the spin-polarized 8, which would correspond to a Landauer conductance of $`5G_0`$ and $`4G_0`$, respectively, if all channels conducted fully. ($`G_0=2e^2/h`$ is the fundamental conductance quantum.) However, the bands splitting around the Fermi level due to the magnetic state are of mainly $`d`$-character, which conduct poorly compared to $`s`$ channels. This means that the conductance of the magnetic and nonmagnetic wires will be much lower than these numbers indicate, and the difference in conductance will be significantly smaller than $`G_0`$. Conductance calculations for a magnetic Pt nanowire segment between tips are presently being considered, but are beyond the scope of this work.

In summary, we find that infinitely long Pt monowires have a ferromagnetic ground state, with a moment of around 0.6 $`\mu _B`$ at the equilibrium bond length 2.48 Å. The moment increases with stretching, and the trigger of the moment formation is the $`(d_{xy},d_{x^2y^2})`$ orbital. The resulting superparamagnetic state of the nanowire will show up in the ballistic conductance in the form of a strong and unusual magnetic field and temperature dependence. Also, more majority bands cross the Fermi level than do minority bands, resulting in a partial spin-polarization of the transmitted electron current. If this current could be measured, it would be a very direct way of confirming the existence of a superparamagnetic state. Rodrigues et al. recently measured the charge conductance of Pt nanocontacts and found features above as well as below $`G_0`$. More theory work will be needed to address their data, explicitly including such elements as tips and temperature.
Acknowledgments
A.D. acknowledges financial support from the European Commission through contract no. HPMF-CT-2000-00827 Marie Curie fellowship, STINT (The Swedish Foundation for International Cooperation in Research and Higher Education), and NFR (Naturvetenskapliga forskingsrådet). Work at SISSA was also sponsored through TMR FULPROP, MUIR (COFIN and FIRB RBAU01LX5H) and by INFM/F. Ruben Weht is acknowledged for discussions, and for double-checking some of the calculations using the WIEN97 code. J. M. Wills is acknowledged for letting us use his FP-LMTO code. We are also grateful to D. Ugarte for sharing with us the unpublished results of Ref. .