text:
# Symmetry of anisotropic exchange interactions in semiconductor nanostructures

## Abstract

The symmetry of exchange interaction of charge carriers in semiconductor nanostructures (quantum wells and quantum dots) is analysed. It is shown that the exchange Hamiltonian of two particles belonging to the same energy band can be universally expressed via pseudospin operators of the particles. The relative strength of the anisotropic exchange interaction is shown to be independent of the binding energy and the isotropic exchange constant.

The reduced symmetry of semiconductor nanostructures suggests that the exchange interaction of charge carriers in such structures is not necessarily described by the isotropic (Heisenberg) spin Hamiltonian. In particular, the exchange interaction of electrons and holes in quantum wells and quantum dots is known to be extremely anisotropic, giving rise to a fine structure of nanostructure excitons . It has been shown  that the exchange interaction of conduction-band electrons is also anisotropic if the structure lacks inversion symmetry. The main term of the anisotropic exchange Hamiltonian in this case has the Dzyaloshinskii-Moriya form. The electron-electron anisotropic exchange was subsequently widely discussed in relation to the quantum computing problem . It has been recently detected experimentally via its contribution to the spin relaxation of donor-bound electrons in GaAs , where it has been shown to put the upper boundary for the electron spin lifetime at donor concentrations around 1016 cm-3 . However, it remained so far unclear whether or not the anisotropic spin Hamiltonian suggested in Ref. is universal for all types of charge carriers, e.g. for two-dimensional holes. The issue of the dependence of the anisotropic exchange constant on the parameters of the localizing potential is also very sensitive, especially for the discussion on feasibility of quantum computation with solid-state spin systems . The constant was so far calculated using the Heitler-London method , which is known to give incorrect asymptotic expression for the isotropic exchange integral . Gor’kov and Krotkov , using the median-plane method, have recently obtained a correct asymptotic formula for the anisotropic exchange constant in a specific case of hydrogen-like centers in zinc-blende semiconductors, different from that calculated earlier by the Heitler-London method . However their approach is not always applicable to coupled quantum dots, where the distance between quantum dots can be comparable to the quantum-dot size.

The collection of unsolved problems and blanc spaces in the existing knowledge on the anisotropic exchange in semiconductor structures, given above, demonstrates the evident demand for a consistent theoretical analysis of the issue, based on a general approach. In this paper, we consider exchange interaction of two identical charge carriers localized in any symmetric double-well potential in a two-dimensional semiconductor structure. Using the pseudospin formalism allows to obtain a universal spin Hamiltonian describing this class of systems.

Let us consider the exchange interaction of two identical charge carriers (electrons or holes), localized in two centrosymmetric potential hollows (further referred to as quantum dots, QDs) in a quasi-two-dimensional semiconductor structure (quantum well, QW). The QDs may be, for example, self-organized QDs ; otherwise, they can be induced by electrostatic potential of nanometer-sized gates  or impurity centers . The distance between centers of the QDs will further be denoted as $`R_{12}`$. In quasi-two-dimensional structures, the 4-fold degeneration of the valence band, typical of cubic semiconductors, is lifted. The states at extremum points of two-dimensional subbands in absence of magnetic fields retain only the Kramers two-fold degeneration. Their wave functions can be written as $`\mathrm{\Psi }(𝐫)u_\nu (𝐫)`$, where is an envelope function, and $`u_\nu (𝐫)`$ is a Bloch amplitude, $`\nu =\pm 1/2`$. The Bloch amplitudes $`u_\nu (𝐫)`$ are transformed into each other by the operator of time reversal :

$$u_{1/2}(𝐫)=i\widehat{\sigma }_yu_{+1/2}(𝐫)$$
(1)
This property allows to associate the Kramers index $`\nu `$ with an eigenvalue of a projection of a pseudospin operator $`𝐣`$ ($`j=1/2`$) on some (generally, fictitious) axis. The choice of basis functions for $`𝐣`$ is not unambiguous. It is limited only by the condition given by Eq.(1). In particular, for heavy holes with the projection of the angular momentum on the structure axis Z, equal to $`J_z=\pm 3/2`$, it is convenient to choose the functions as : $`|j,+1/2=|J_z,3/2`$ and $`|j,1/2=|J_z,+3/2`$. This choice allows to avoid phase multipliers which would otherwise appear at wave functions in the pseudospin representation. For conduction-band electrons, the pseudospin coincides with the electron spin $`𝐬`$. Linear transformations of pseudospin wave functions determined in the basis $`\{u_{+1/2},u_{1/2}\}`$ are equivalent to rotations of usual spinor functions :

$`u_{+1/2}^{\alpha \beta \gamma }`$ $`=`$ $`\mathrm{exp}(i\gamma /2)\left[u_{+1/2}\mathrm{exp}(i\alpha /2)\mathrm{cos}(\beta /2)+u_{1/2}\mathrm{exp}(i\alpha /2)\mathrm{sin}(\beta /2)\right]`$ (2)
$`u_{1/2}^{\alpha \beta \gamma }`$ $`=`$ $`\mathrm{exp}(i\gamma /2)\left[u_{+1/2}\mathrm{exp}(i\alpha /2)\mathrm{sin}(\beta /2)+u_{1/2}\mathrm{exp}(i\alpha /2)\mathrm{cos}(\beta /2)\right]`$ (3)
where $`\alpha `$, $`\beta `$, and $`\gamma `$ are analogs of Euler angles. Following the analogy, one can introduce the total pseudospin $`𝐈=𝐣_1+𝐣_2`$. Indeed, the Gilbert space of two-pseudospin wave functions $`A_{\mu \nu }u_\mu (𝐫_1)u_\nu (𝐫_2)`$ breaks into two subspaces invariant with respect to the simultaneous transformation of both pseudospins along Eq.(2) with the same $`\alpha `$, $`\beta `$, and $`\gamma `$. The basis functions of these subspaces, $`\xi _0=(u_{+1/2}(𝐫_1)u_{1/2}(𝐫_2)u_{1/2}(𝐫_1)u_{+1/2}(𝐫_2))/\sqrt{2}`$ and $`\xi _{1M}`$, equal to $`(u_{+1/2}(𝐫_1)u_{1/2}(𝐫_2)+u_{+1/2}(𝐫_1)u_{1/2}(𝐫_2))/\sqrt{2}`$ ($`M=0`$) or $`u_{\pm 1/2}(𝐫_1)u_{\pm 1/2}(𝐫_2)`$ ($`M=\pm 1`$), are, obviously, eigenfunctions of the operators $`\widehat{I}^2`$ and $`\widehat{I}_z`$.

The general form of the one-particle Hamiltonian of the two dimensional charge carrier in the pseudospin representation is

$$\widehat{H}_1=\frac{\mathrm{}^2}{2m}\widehat{k}^2+V(𝐫)+𝐡(𝐤)𝐣$$
(4)
where the ”spin-orbit field” $`𝐡(𝐤)`$ is a vector in the pseudospin space . $`𝐡(𝐤)`$ is an odd function of the components of the wave vector. It is not equal to zero if the structure lacks inversion symmetry (which is very typical for nanostructures). This is the case when either the crystal unit cell lacks inversion symmetry (bulk inversion asymmetry, BIA ), or the QW is asymmetric (structure inversion asymmetry, SIA ). The components of $`𝐡(𝐤)`$ may be, or may not be, associated with certain Cartesian axes in the real space. In the two-dimensional case, $`𝐡(𝐤)`$ is dominated by linear in $`𝐤`$ terms :

$$h_\eta =A_{\eta \zeta }k_\zeta $$
(5)
where the matrix $`A`$ is defined by the structure symmetry.

The problem we are going to solve is finding the fine structure of the ground state of the two-particle Hamiltonian:

$$\widehat{H}=\widehat{H}_0+\widehat{H}_{SO}$$
(6)
where

$$\widehat{H}_0=\frac{\mathrm{}^2}{2m}\widehat{k}_1^2+\frac{\mathrm{}^2}{2m}\widehat{k}_2^2+V_1(𝐫_1)+V_2(𝐫_2)+U_{12}(\left|𝐫_1𝐫_2\right|)$$
(7)
$$\widehat{H}_{SO}=𝐡(𝐤_1)𝐣_1+𝐡(𝐤_2)𝐣_2$$
(8)
and $`U_{12}(\left|𝐫_1𝐫_2\right|)`$ is the operator of the Coulomb interaction between the two particles. Before tackling the effects of spin-orbit interaction in the form of Eq.(7) on the exchange interaction, we should reconsider the ground-state structure of the Hamiltonian $`H_0`$. It is indeed well-known for electrons whose one-particle wave functions are $`\mathrm{\Psi }(𝐫)\zeta _\mu `$, where the spinor $`\zeta _\mu `$ is not a function of coordinates. To the contrary, the Bloch amplitude $`u_\nu `$ does depend on coordinates, and, moreover, it may contain spinors with both $`\mu =+1/2`$ and $`\mu =1/2`$. The exciton (an electron-hole pair) in a QD is a good example demonstrating that the exchange interaction of charge carriers may have a very different symmetry as compared to that of free electrons. The QD exciton fine structure  consists of two doublets, being thus quite different from the fine structure of a pair of vacuum electrons, i.e. the well-known singlet-triplet structure associated with the Heisenberg exchange.

In order to analyze the fine structure of $`H_0`$ for two particles belonging to the same subband, we first note that their behavior should be identical to that of bare electrons in all aspects but the Coulomb interaction. Indeed, although Bloch amplitudes are functions of coordinates within the unit cell, the one-particle operators of the kinetic energy and of the potential energy in Eq.(6) in the effective-mass approximation act upon envelope function, not Bloch amplitudes. Therefore, with respect to these one-particle operators, the Bloch amplitudes are just equivalent to spinors. To the contrary, calculating the Coulomb energy assumes taking integrals over the unit sell also. It is due this fact that the symmetry of Bloch amplitudes of holes and electrons has an impact on the fine structure in the exciton .

The fermionic wave functions of the two charge carriers can be written in the following form, similar to that of bare electrons:

$`\mathrm{\Psi }_0(𝐫_1,𝐫_2)`$ $`=`$ $`\left[\mathrm{\Phi }_0(𝐫_1,𝐫_2)+\mathrm{\Phi }_0(𝐫_2,𝐫_1)\right]\xi _0`$ (9)
$`\mathrm{\Psi }_{1M}(𝐫_1,𝐫_2)`$ $`=`$ $`\left[\mathrm{\Phi }_1(𝐫_1,𝐫_2)\mathrm{\Phi }_1(𝐫_2,𝐫_1)\right]\xi _{1M}`$ (10)
where $`\mathrm{\Phi }_0(𝐫_1,𝐫_2)`$ and $`\mathrm{\Phi }_1(𝐫_1,𝐫_2)`$ are two-particle envelope functions defined so that each particle is most likely to be found near its ”home” center, while $`\mathrm{\Phi }_0(𝐫_2,𝐫_1)`$ and $`\mathrm{\Phi }_1(𝐫_2,𝐫_1)`$ correspond to interchanged particle positions. To determine the structure of respective energy levels, we should recall a property of the Bloch amplitudes $`u_\nu `$, which results from their symmetry with respect to time reversal, and is an equivalent formulation of the Kramers theorem. As it follows from the Kramers theorem, the states symmetric with respect to time reversal remain degenerated unless magnetic field is applied. Mathematically, this means that matrix elements of any function of coordinates (not containing derivatives or spin operators) between $`u_{+1/2}(𝐫)`$ and $`u_{1/2}(𝐫)`$ are zero, while diagonal matrix elements are equal to each other:

$`{\displaystyle \underset{\mathrm{\Omega }}{}}u_{+1/2}(𝐫)u_{1/2}^{}(𝐫)f(𝐫)d^3r`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Omega }}{}}u_{1/2}(𝐫)u_{+1/2}^{}(𝐫)f(𝐫)d^3r=0`$ (11)
$`{\displaystyle \underset{\mathrm{\Omega }}{}}u_{+1/2}(𝐫)u_{+1/2}^{}(𝐫)f(𝐫)d^3r`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Omega }}{}}u_{1/2}(𝐫)u_{1/2}^{}(𝐫)f(𝐫)d^3r`$ (12)
where the integrals are taken over the unit cell.

Using Eq.(9), one can easily find that

$`\mathrm{\Psi }_0(𝐫_1,𝐫_2)|U_{12}(\left|𝐫_1𝐫_2\right|)|\mathrm{\Psi }_{1M}(𝐫_1,𝐫_2)`$ $`=`$ $`0`$ (13)
$`\mathrm{\Psi }_0(𝐫_1,𝐫_2)|U_{12}(\left|𝐫_1𝐫_2\right|)|\mathrm{\Psi }_0(𝐫_1,𝐫_2)`$ $``$ (14)
$`\mathrm{\Psi }_{1M}(𝐫_1,𝐫_2)|U_{12}(\left|𝐫_1𝐫_2\right|)|\mathrm{\Psi }_{1M^{}}(𝐫_1,𝐫_2)`$ $`=`$ $`const\delta _{MM^{}}`$ (15)
Thus, the Coulomb interaction retains the singlet-triplet structure of the ground state of two identical charge carriers. Exactly like in the case of two bare electrons, two-particle states with the same total pseudospin $`I`$ are degenerated. Consequently, the Hamiltonian of the exchange interaction in terms of pseudospin operators takes the Heisenberg form:

$$\widehat{H}_S=2\mathrm{\Delta }(𝐣_1𝐣_2+1)$$
(16)
where $`\mathrm{\Delta }`$ is a constant to be determined for each specific case.

Now we can consider the effect of the spin-orbit terms given by the Eq.(7) on the exchange interaction. In the following, we will choose the axis X along the straight line connecting the localization centers (QDs). To handle the spin-orbit terms, we make use of a unitary transformation proposed by Levitov and Rashba  who used it to eliminate spin-orbit terms in the one-dimensional case. The matrix $`T`$ defined as

$$T=\mathrm{exp}\left[i\frac{2m}{\mathrm{}^2}\underset{\alpha }{}A_{\alpha x}(\widehat{ȷ}_{1\alpha }x_1+\widehat{ȷ}_{2\alpha }x_2)\right]$$
(17)
transforms the Hamiltonian Eq.(5) into the form:

$$T\widehat{H}T^1=\widehat{H}^{}=\widehat{H}_0^{}+\widehat{H}_{SO}^{}$$
(18)
where

$$\widehat{H}_0^{}=\frac{\mathrm{}^2}{2m}\widehat{k}_1^2+\frac{\mathrm{}^2}{2m}\widehat{k}_2^2+V_1(𝐫_1^{})+V_2(𝐫_2^{})+U_{12}(\left|𝐫_1^{}𝐫_2^{}\right|)\underset{\alpha }{}\frac{m(A_{\alpha x})^2}{\mathrm{}^2}$$
(19)
and

$$\widehat{H}_{SO}^{}=\underset{\alpha }{}A_{\alpha y}(k_{1y}\widehat{ȷ}_{1\alpha }^{}+k_{2y}\widehat{ȷ}_{2\alpha }^{})$$
(20)
where $`𝐣_1^{}=T𝐣_1T^1`$, $`𝐣_2^{}=T𝐣_2T^1`$.

The Hamiltonian $`\widehat{H}^{}`$ does not contain spin-orbit terms and therefore results in the exchange interaction in the form of Eq.(11):

$$\widehat{H}_S^{}=2\mathrm{\Delta }(𝐣_1^{}𝐣_2^{}+1)$$
(21)
Due to the axial symmetry of the system, the matrix elements of $`k_{1y}`$ and $`k_{2y}`$, calculated on the ground-state eigenfunctions of $`\widehat{H}_0^{}`$ (they can be obtained from Eq.(8) by the transformation with the matrix $`T`$, which does not affect their dependence on $`y`$), are exactly equal to zero. The same is true for all the odd powers of $`k_{1y}`$ and $`k_{2y}`$. Therefore, $`\widehat{H}_{SO}^{}`$ does not contribute into the exchange interaction.

Finally, to obtain the exchange Hamiltonian in the non-transformed basis, one should substitute the expressions for $`𝐣_1^{}`$ and $`𝐣_2^{}`$ into Eq.(16). Since the transformation $`T`$ is a rotation through the angle $`\frac{2m}{\mathrm{}^2}\sqrt{\underset{\nu }{}A_{\nu x}A_{\nu x}}x`$ around the vector $`A_{\nu x}`$ in the pseudospin space, $`\widehat{H}_S`$ in the non-transformed basis is not unambiguously defined: it depends on the coordinates $`x_1`$ and $`x_2`$ at which we take the spin operators $`𝐣_1`$ and $`𝐣_2`$. A natural choice is to define them at the centers of corresponding QDs; for instance, this definition allows to write the Zeeman interaction in the usual form, $`\widehat{H}_Z=\mu _Bg_{\alpha \beta }j_\alpha B_\beta `$, where $`B`$ is the magnetic field, $`\mu _B`$ is the Bohr magneton, and $`g_{\alpha \beta }`$ is a symmetric tensor g-factor whose principal directions do not depend on the envelope wave function of the localized particle. This way, we come to the expression for $`\widehat{H}_S`$ obtained in Ref.:

$$\widehat{H}_S=2\mathrm{\Delta }(1+𝐣_1𝐣_2\mathrm{cos}\gamma +(𝐝𝐣_\mathrm{𝟏})(𝐝𝐣_\mathrm{𝟐})(\mathrm{𝟏}\mathrm{cos}\gamma )+𝐝\left[𝐣_1\times 𝐣_2\right]\mathrm{sin}\gamma )$$
(22)
where $`\gamma =\frac{2m}{\mathrm{}^2}\sqrt{\underset{\nu }{}A_{\nu x}A_{\nu x}}R_{12}`$, and $`𝐝`$ is a unit vector in the pseudospin space, defined so that $`d_\nu =A_{\nu x}/\sqrt{\underset{\nu }{}A_{\nu x}A_{\nu x}}`$. The first anisotropic term has the form of pseudodipole interaction , and the second one, of the Dzyaloshinskii-Moriya interaction . At small $`\gamma `$, the Dzyaloshinskii-Moriya interaction dominates.

The Eq.(17) demonstrates a remarkable universality of the exchange interaction in two-dimensional semiconductor nanostructures: this form of the Hamiltonian holds for both electrons and holes, for any type of centrosymmetric localizing potentials. The Eq.(17) is valid for identical as well as for different QDs. Moreover, the angle $`\gamma `$ characterizing the relative strength of the anisotropic exchange depends only on the distance between the QDs and the orientation of the pair of QDs with respect to the crystal axes. It is not sensitive to binding energies of the charge carriers in the QDs and to the value of the isotropic exchange constant $`\mathrm{\Delta }`$.

The value of $`\gamma `$ can be now easily calculated for those structures where the components of the matrix $`A`$ are known.

In  oriented GaAs quantum wells the dominating BIA terms are  $`A_{yy}=A_{xx}=\frac{\alpha \mathrm{}^3}{m\sqrt{2mE_g}}k_z^2`$, $`A_{xy}=A_{yx}=0`$ where $`\alpha 0.65`$  (here coordinates $`x`$ and $`y`$ are taken along the cubic crystal axes). This gives $`\gamma =\frac{2\alpha \mathrm{}}{\sqrt{2mE_g}}k_z^2R_{12}`$. For example, in a 5nm-wide QW considered in Ref., $`\gamma =(310^5cm^1)R_{12}`$, which gives for $`R=24nm`$ $`\gamma =0.72`$, i.e. more than an order of magnitude greater than the value of 0.02 obtained in Ref. at the same values of parameters by the Heitler-London method.

For Rashba terms, $`A_{yy}=A_{xx}=0`$, $`A_{xy}=A_{yx}=a`$. In a single-side modulation-doped n-type Si/SiGe quantum well, the constant $`a`$ of 1.1$``$10-12 eV$``$cm was measured . This gives $`\gamma =\frac{2ma}{\mathrm{}^2}R_{12}(6.710^2cm^1)R_{12}`$.

Bulk inversion asymmetry terms for holes in zinc-blende semiconductors include both cubic and linear in $`k`$ terms. The cubic term $`H_{3V}`$ has the same symmetry as the Dresselhaus term for electrons, with the constant $`\frac{\alpha _V\mathrm{}^3}{m\sqrt{2mE_g}}`$, where $`m`$ is the conduction-band electron mass, and $`\alpha _V0.1`$ for GaAs. The linear term is given by the expression:

$$H_{1V}=\frac{4}{\sqrt{3}}\varkappa \left(𝐤𝛀\right)$$
(23)
where $`\mathrm{\Omega }_z=\widehat{J}_z\left(\widehat{J}_x^2\widehat{J}_y^2\right)+\left(\widehat{J}_x^2\widehat{J}_y^2\right)\widehat{J}_z`$ (other components of $`𝛀`$ are obtained by cyclic interchange of indices), and $`\varkappa 10^{10}eVcm`$. Taking matrix elements of $`H_{1V}`$ and $`H_{3V}`$ within pairs of the states with $`J_z=\pm 1/2`$ (light holes) and $`J_z=\pm 3/2`$ (heavy holes), and going to the pseudospin notation, we obtain for a  QW:

$`𝐡_l(𝐤)`$ $`=`$ $`\left(2{\displaystyle \frac{\alpha _V\mathrm{}^3}{m\sqrt{2mE_g}}}k_z^22\sqrt{3}\varkappa \right)\left(k_x\widehat{ȷ}_xk_y\widehat{ȷ}_y\right)`$ (24)
$`𝐡_h(𝐤)`$ $`=`$ $`4\sqrt{3}\varkappa \left(k_x\widehat{ȷ}_xk_y\widehat{ȷ}_y\right)`$ (25)
where $`𝐡_l(𝐤)`$ and $`𝐡_h(𝐤)`$ are spin-orbit fields (see Eq.(3)) for light and heavy holes respectively. Consequently, $`\gamma _l=\left(\frac{4\alpha _V\mathrm{}m_l}{m\sqrt{2mE_g}}k_z^2\frac{4\sqrt{3}m_l}{\mathrm{}^2}\varkappa \right)R_{12}`$ for light holes, and $`\gamma _h=\frac{8\sqrt{3}m_h}{\mathrm{}^2}\varkappa R_{12}`$ for heavy holes, where $`m_l`$ and $`m_h`$ are effective masses of the light and heavy hole respectively, corresponding to their motion along the QW plane. For example, in a 10nm-wide GaAs QW, $`\gamma _l\gamma _h(210^5cm^1)R_{12}`$.

The symmetry of exchange interactions has been discussed in relation to feasibility of quantum computation with spins of localized electrons in semiconductor nanostructures . A necessary (but not sufficient; see Ref.) condition for practical quantum computing to become possible is that the error probability per quantum gate be less than a certain value (of the order of 10-5) . As shown in Ref., there exists a way of performing exchange-mediated quantum gates that allows to avoid errors caused by the anisotropy, provided $`\gamma `$ remains constant when the isotropic exchange constant is changed. The above consideration shows that this is indeed the case as long as spin-orbit constants do not change. They may change, however, because application of electric fields to the structure can either alter $`k_z^2`$ or bring about SIA (Rashba) terms. Since typical values of $`\gamma `$ are in reality much greater than estimated in Ref., the uncontrollable effect of anisotropy on exchange-mediated quantum gates can not be so easily discarded.

In conclusion, the exchange interaction of charge carriers (electrons or holes) localized in two-dimensional semiconductor structures is shown to be described by an universal Hamiltonian in terms of carriers’ pseudospins. It has the Heisenberg form unless spin-orbit terms, linear in the carrier wave vector, are present in the total Hamiltonian of the system. In this latter case, anisotropic contributions having both Dzyaloshinskii-Moriya and pseudodipole form arise. The ”rotation angle” $`\gamma `$, characterizing the relative strength of the anisotropic exchange, linearly depends on the distance between the localization centers and does not depend on binding energies of the carriers.

The author is grateful to I.A.Merkulov, E.L.Ivchenko, L.P.Gor’kov, and P.G.Krotkov for stimulating discussions. Partial support of INTAS, RFBR, and UK EPSRC is acknowledged.


format: text

text:
# Higher order clustering coefficients in Barabási-Albert networks

## 1 Introduction

During the last few years studies of random, evolving networks ( such as the Internet, WWW, social networks, metabolic networks, food webs etc - for review see ) have become a very popular research domain among physicists. A lot of efforts were put into investigation of such systems in order to recognize their structure and to analyse emerging complex dynamics. We learned that networks are far from being random as Erdös and Rényi assumed in random graph theory , but surely they are not so ordered as crystals. Despite network diversity, most of real weblike systems share three prominent features :

* The average shortest path length $`l`$ is small. In order to connect two nodes in a network typically only a few edges need to be passed.
* The average clustering coefficient $`C`$ is large. Two nodes having a common neighbour are also likely to be neighbours.
* The probability that a randomly selected node has exactly $`k`$ nearest neighbours follows a power law (scale-free) distribution $`P(k)k^\gamma `$ with $`2<\gamma <3`$ in most of real systems.

A considerable number of network models has been studied in order to capture the above characteristics. Most of these are based on two ingredients originally introduced by Barabási and Albert : continuous network growth and preferential attachment. In Barabási-Albert (BA) model, network starts to grow from an initial cluster of $`m`$ fully connected sites. Each new node that is added to the network creates $`m`$ links that connect it to previously added nodes. The preferential attachment means that the probability of a new link to end up in a vertex $`i`$ is proportional to connectivity $`k_i`$ of this vertex. The validity of the preferential attachment was confirmed within real networks analyses . The BA algorithm generates networks with the desirable scale-free distribution $`P(k)k^3`$ and small values of the average shortest path. One can also observe a phase transition for spins located at BA network vertices with a critical temperature increasing as a logarithm of the system size . The only striking discrepancy between the BA model and real networks is that the value of the clustering coefficient predicted by the theoretical model decays very fast with network size and for large systems is typically several orders of magnitude lower than found empirically.

In this paper we extend the standard definition of the clustering coefficient by introducing higher order clustering coefficients that describe interrelations between vertices belonging to the nearest neighbourhood of a certain vertex in complex network. Global characteristics like the standard clustering coefficient and the average shortest path do not provide a useful insight into complex network structure and dynamics. We hope that the higher order clustering coefficient analyses in real systems  may give some guidelines how to model clustering mechanisms. Here, we study higher order clustering coefficients in BA model. Our results provide a vivid evidence that the BA networks are blind to clustering mechanisms.

## 2 Model description

The standard clustering coefficient $`C`$ is one of global parameters used to characterise the topology of complex networks. For the first time it was introduced by Watts and Strogatz  to characterise local transitivity in social networks. Clustering coefficient gives the probability that two nearest neighbours of the same node are also mutual neighbours. Let us focus on a selected node $`i`$ in a network, having $`k_i`$ edges which connect it to $`k_i`$ other nodes. The value of the clustering coefficient of the node $`i`$ is given by the ratio between the number of edges $`E_i`$ that actually exist between these $`k_i`$ nodes and the total number $`k_i(k_i1)/2`$ of such edges that could exist in the neighbourhood of $`i`$:

$$C_i=\frac{2E_i}{k_i(k_i1)}.$$
(1)
The clustering coefficient of the whole network is the average of all individual $`C_i`$’s

We define a clustering coefficient of order $`x`$ for a node $`i`$ as the probability that there is a distance of length $`x`$ between two neighbours of a node $`i`$. Putting the number of such $`x`$-distances equal to $`E_i(x)`$, the higher order clustering coefficients follow:

$$C_i(x)=\frac{2E_i(x)}{k_i(k_i1)}.$$
(2)
$`C(x)`$ is the mean value of $`C_i(x)`$ over the whole network. Note that $`C(x)`$ reduces to the standard clustering coefficient for $`x=1`$ and $`_xC(x)=1`$ for the BA networks with $`m2`$.

The above definition becomes comprehensible after examining Fig.1. Let us assume a node $`i`$ with $`k_i=5`$. The node determines its nearest neighbourhood that in this case consists of vertices $`\{1,2,3,4,5\}`$. The aim is to find higher order clustering coefficients $`C_i(x)`$ of the node $`i`$. The table below gives the shortest distances between vertices adjacent to the node $`i`$. The distances were taken from the Fig.1 (Stage $`2`$):

$$\begin{array}{ccccccccc}& & 1& 2& 3& 4& 5& & \\ 1& (\mathrm{}& & 1& 1& 4& 1& )\mathrm{}\\ 2& 1& & 1& 4& 2\\ 3& 1& 1& & 3& 2\\ 4& 4& 4& 3& & 5\\ 5& 1& 2& 2& 5& \end{array}.$$

Summing over all pairs of vertices one obtains:

$$\begin{array}{cccc}& & & \\ x& (\mathrm{}& 1& 2& 3& 4& 5& )\mathrm{}\\ E_i(x)& 4& 2& 1& 2& 1\\ C_i(x)& 0.4& 0.2& 0.1& 0.2& 0.1\end{array}.$$

## 3 Numerical results

Fig.2 shows the higher order clustering coefficients dependence on the network size $`N`$. Nodes determined as belonging to the nearest neighbourhood of any vertex in BA network are more likely to be second ($`x=2`$), third ($`x=3`$) and further neighbours when $`N`$ increases.

We investigated distributions of the clustering coefficients $`C(x)`$ at a given system size $`N`$ in BA networks with $`m=3`$ and $`4`$. We found that regardless of $`N`$ the distributions of $`C(x)`$ (when extending $`x`$ to real values) fairly good fit the normalized Gaussian curve (Fig.3). Moreover, we observed that the standard deviation of these distributions depends only on the parameter $`m`$ of the network. The gaussian-like patterns of $`C(x)`$ in BA model express random character of clustering relationships. In fact, it is known that the distribution of distances between randomly chosen sites in BA network is Gaussian. It follows that interrelations in the nearest neighbourhood of any vertex in BA model are similar to interrelations between randomly chosen vertices in the BA network. Distinctly, there are no mechanisms responsible for clustering in Barabási-Albert algorithm.

The above-described observation can be verified by a functional dependence of centers $`x_c`$ of the distribution $`C(x)`$ on the network size $`N`$. We realised that the value $`x_c`$ (Fig.3) expresses both: the order $`x`$ when $`C(x)`$ obtains maximum and the average shortest path length $`l_{cluster}`$ between vertices belonging to the neighbourhood of any vertex within a network. Fig.4 shows that in BA networks the average shortest path $`l_{cluster}`$ scales with the system size $`N`$ in the same way as the average shortest path length for the whole network $`l_{network}`$ , i.e. in the first approximation as:

$$l_{cluster}\mathrm{ln}(N).$$
(3)
We found that the mean distance between two nearest neighbours of the same node $`l_{cluster}`$ equals to the mean distance between any two nodes in the network $`l_{network}`$ minus a constant $`A(m)`$

$$l_{cluster}=l_{network}A(m),$$
(4)
where $`A(m)`$ is approximately independent on the network size $`N`$ and equals $`0.37`$ and $`0.27`$ for $`m=3`$ and $`m=4`$ respectively.

## 4 Conclusions

In summary, we quantified the structural properties of BA networks by the higher order clustering coefficients $`C(x)`$ defined as probabilities that the shortest distance between any two nearest neighbours of a certain vertex $`i`$ equals $`x`$, when neglecting all paths crossing the node $`i`$. We estimated that the average shortest path length in the node’s neighbourhood is smaller than the equivalent whole network quantity and the remainder depends only on the network parameter $`m`$. Our results show in a vivid way that the absence of the clustering phenomenon in BA networks is due to the random character of the nearest neighbourhood in these networks.

Recently some alternative algorithms have been suggested to account for the high clustering found in real weblike systems. Holm and Kim  have extended the standard BA model adding the trial formation rule. Networks built according to their guidelines exhibit both the high clustering and the scale-free nature. Barabasi et al.  and independently Dorogovtsev et al.  have found that scale-free random networks can be modelled in a deterministic manner by so-called pseudofractal scale-free networks. They argued that the high clustering in real networks might result from their hierarchical topology. It would be interesting to analyse the higher order clustering coefficients in these networks and compare it with real data .


text:
# References

KIMS-2002-12-23

physics/0212092

Inconsistent Universe

– Physics as a meta-science –

Hitoshi Kitada

Graduate School of Mathematical Sciences

University of Tokyo

Komaba, Meguro-ku, Tokyo 153-8914, Japan

e-mail: kitada@ms.u-tokyo.ac.jp December 23, 2002

Abstract: Physics is introduced as a semantics of a formal set theory.

We consider a metatheory of a formal set theory $`S`$. We name this metatheory $`M_S`$, indicating that it is a Meta-theory of $`S`$ as well as a Meta-Scientific theory as Ronald Swan  refers to. The following arguments are all made in $`M_S`$.

We define in $`M_S`$

$$\varphi =\text{ the class of all well-formed formulae of }S.$$

This $`\varphi `$ is a countable set in the context of $`M_S`$.

We identify $`\varphi `$ with the set of truth values (in complex numbers C) of well-formed formulae (wff’s) in $`\varphi `$. In this identification, we define a map $`T`$ from $`\varphi `$ to $`\varphi `$ by

$$T((q))=\text{the truth value of a well-formed formula }[(q)\text{ and not }(q)]$$

for $`q\varphi `$, with $`(q)`$ denoting the conjunction of $`q`$.

We note that every subset $`q`$ of $`\varphi `$ becomes false by adding some well-formed formula $`f`$ of $`\varphi `$. Hence, the conjunction of $`q^{}=q\{f\}`$ is false and satisfies

$$T((q^{}))=(q^{}).$$

In this sense, $`\varphi `$ is a fixed point of the map $`T`$.

Moreover, we have the followings.

1. In the sense that any subset $`q`$ of $`\varphi `$ is false if some well-formed formula is added to $`q`$, $`\varphi `$ is inconsistent.
2. As $`\varphi `$ is the class of all possible well-formed formulae, $`\varphi `$ is absolute.
3. As $`\varphi `$ is the totality of well-formed formulae, $`\varphi `$ includes the well-formed formula whose meaning is that “$`\varphi `$ is the class of all well-formed formulae in $`S`$” in some Gödel type correspondence between $`S`$ and $`M_S`$. In this sense $`\varphi `$ includes (the definition of) $`\varphi `$ itself. Thus $`\varphi `$ is self-referential and self-creative, and is self-identical, just as in M. C. Escher’s lithograph in 1948, entitled “pencil drawing.”

The item 3 implies that $`\varphi `$ is a non-well founded set or a hyperset.

The class $`\varphi `$ is the first world, the Universe, which is completely chaotic. In other words, $`\varphi `$ is “absolute inconsistent self-identity” in the sense of Kitarou Nishida , whose meaning was later clarified by Ronald Swan  in the form stated above. In this clarification, $`\varphi `$ can be thought “absolute nothingness” in Hegel’s sense.

The Universe $`\varphi `$ is contradictory, and hence its truth value is constantly oscillating between the two extremal values or poles, truth and false, or $`+1`$ and $`1`$, or more generally, inside a unit sphere of C. Namely, the class $`\varphi `$ as a set of wff’s of the set theory $`S`$ is countable, but the values which the elements of $`\varphi `$ take vary on a unit sphere. In other words, the Universe $`\varphi `$ is a stationary oscillation, when we see its meaning.

Oscillation is expressed by exponential functions: $`\mathrm{exp}(ixp)`$, where $`x=(x_1,\mathrm{},x_d),p=(p_1,\mathrm{},p_d)\text{R}^d`$ and $`xp=_{i=1}^dx_ip_i`$.

This $`\mathrm{exp}(ixp)`$ is an eigenfunction of the negative Laplacian $`\mathrm{\Delta }`$:

$$\mathrm{\Delta }=\underset{i=1}{\overset{d}{}}\frac{^2}{x_i^2}.$$

Namely

$$\mathrm{\Delta }\mathrm{exp}(ixp)=p^2\mathrm{exp}(ixp).$$

This is generalized to some extent. I.e. if a perturbation $`V=V(x)`$ satisfies that

$$H=\mathrm{\Delta }+V(x)\text{is a self-adjoint operator on}=L^2(\text{R}^d),$$

then

$`\varphi `$ is expressed as an eigenfunction of $`H`$.

Considering the absolute nature of the Universe $`\varphi `$, we will be led to think that the Hamiltonian $`H`$ of $`\varphi `$ is a Hamiltonian of infinite degree of freedom on a Hilbert space:

$$𝒰=\{\varphi \}=\underset{n=0}{\overset{\mathrm{}}{}}\left(\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}^n\right)(^n=\underset{n\text{ factors}}{\underset{}{\mathrm{}}}),$$

and the Universe $`\varphi `$ is an eigenfunction of the total Hamiltonian $`H=H_{\text{total}}`$.

Thus we arrive at our first principle.

Axiom 1. The Universe is of infinite nature, and it is eternal. In other words, the wave function $`\varphi `$ of the Universe satisfies with respect to the total Hamiltonian $`H_{\text{total}}`$

$$H_{\text{total}}\varphi =\lambda \varphi $$

for some non-positive real number $`\lambda 0`$.

In every finite part of $`\varphi `$, a local existence in $`\varphi `$ is expressed by a superposition of exponential functions

$$\psi (x)=(2\pi )^{d(N1)/2}_{\text{R}^{d(N1)}}\mathrm{exp}(ixp)g(p)𝑑p$$

for some natural number $`N=n+12`$ with $`n`$ corresponding to the superscript $`n`$ in $`^n`$ of the definition of $`𝒰`$ above. The function $`g(p)`$ is called Fourier transform of $`\psi (x)`$ and satisfies

$$g(p)=\psi (p):=(2\pi )^{d(N1)/2}_{\text{R}^{d(N1)}}\mathrm{exp}(ipy)\psi (y)𝑑y.$$

A finite subset of wff’s in $`\varphi `$ corresponds to a partial Hamiltonian $`H`$ of $`H_{\text{total}}`$ of finite degree of freedom, as the content/freedom that is given by a finite number of wff’s in $`\varphi `$ corresponds to a finite degree, $`n=N1`$, of freedom of a partial wave function $`\psi (x)`$ of the total wave function $`\varphi `$. If such a partial Hamiltonian $`H`$ of $`H_{\text{total}}`$ satisfies some conditions, we can get a similar expansion of a local existence $`\psi (x)`$ by using generalized eigenfunctions of $`H`$. This is known as a spectral representation of $`H`$ in a general setting, but we here are speaking of a more specific expression called generalized Fourier transform or generalized eigenfunction expansion associated with Hamiltonian $`H`$ (originated by Teruo Ikebe ).

We call $`p`$ momentum conjugate to $`x`$. More precisely we define momentum operator $`P=(P_1,\mathrm{},P_d)`$ conjugate to configuration operator $`X=(X_1,\mathrm{},X_d)`$ $`(X_j=\text{multiplication operator by configuration }x_j)`$ by

$$P_j=^1p_j=\frac{1}{i}\frac{}{x_j}(j=1,\mathrm{},d).$$

Then $`P`$ and $`X`$ satisfy

$$[P_j,X_{\mathrm{}}]=P_jX_{\mathrm{}}X_{\mathrm{}}P_j=\delta _j\mathrm{}\frac{1}{i}.$$

This shows that what we are dealing with is quantum mechanics. So to be in accordance with actual observation, we modify the definition of $`P`$

$$P_j=\frac{\mathrm{}}{i}\frac{}{x_j},$$

where $`\mathrm{}=h/(2\pi )`$, and $`h`$ is Planck constant. Accordingly, the Fourier and inverse Fourier transformations are modified

$`\psi (p)=g(p)=(2\pi \mathrm{})^{d(N1)/2}{\displaystyle _{\text{R}^{d(N1)}}}\mathrm{exp}(ipy/\mathrm{})\psi (y)𝑑y,`$
$`^1g(x)=(2\pi \mathrm{})^{d(N1)/2}{\displaystyle _{\text{R}^{d(N1)}}}\mathrm{exp}(ixp/\mathrm{})g(p)𝑑p.`$
To sum our arguments up to here, we have constructed quantum mechanics as a semantics of the class $`\varphi `$ of all well-formed formulae of a formal set theory $`S`$. Quantum mechanics is, in this context, given as an interpretation of set theory.

We continue to complete our semantics of the Universe $`\varphi `$.

A local existence is of finite nature, and it is so local that it cannot know the existence of the infinite Universe, and is self-centered. In other words, a local coordinates system starts from its own origin, and it is the self-centered origin of the local system. All things are measured with respect to this local origin.

Therefore we have our second and third principles.

Axiom 2. A local system is of finite nature, having its own origin of position $`X`$ and momentum $`P`$, independent of others’ origins and others’ inside worlds.

Axiom 3. The nature of locality is expressed by a local Hamiltonian

$$H=\frac{1}{2}\mathrm{\Delta }+V$$

up to some perturbation $`V`$, that does not violate the oscillatory nature of local existence. Here $`\mathrm{\Delta }=_{j=1}^{N1}\frac{\mathrm{}^2}{\mu _j}_{k=1}^d\frac{^2}{x_{jk}^2}`$, the number $`N`$ corresponds to the number of quantum particles of the local system, and $`\mu _j`$ is the reduced masses of the particles of the local system.

A local existence (or local system) is oscillating as a sum or integral of generalized eigenfunctions of $`H`$. In this sense, the locality or local system is a stationary oscillating system.

A local oscillation may be an eigenfunction of the local Hamiltonian $`H`$. However, by the very nature that locality is a self-centered existence of finite nature, it is shown that it cannot be an eigenstate of $`H`$, or more precisely speaking, there is at least one Universe wave function $`\varphi `$ every part of which is not an eigenfunction of the local system Hamiltonian $`H`$ corresponding to the part. (See , . See also , , .)

To express this oscillation explicitly in some “outer coordinate,” we force the locality or local system to oscillate along an “afterward-introduced” real-valued parameter $`t`$. The oscillation is then expressed by using the Hamiltonian $`H`$

$$\mathrm{exp}(2\pi itH/h).$$

This operator is known in QM (quantum mechanics) as the evolution operator of the local system. We call it the local clock of the system, and we call $`t`$ the local time of the system.

Using our self-centered coordinates of our local system in axiom 2, that is, letting $`x`$ be position coordinates and $`v=m^1P`$ be velocity coordinates inside the local system ($`m`$ being some diagonal mass matrix), we can prove, by virtue of the fact that a local oscillation $`\psi (x)`$ is not an eigenfunction of $`H`$, that

$$\left(\frac{x}{t}v\right)\mathrm{exp}(itH/\mathrm{})\psi (x)0$$

as $`t`$ tends to $`\pm \mathrm{}`$ along some sequence in some spectral decomposition of $`\mathrm{exp}(itH/\mathrm{})\psi `$ (see ). This means that the word “local clock” is appropriate for the operator $`\mathrm{exp}(itH/\mathrm{})`$ and so is “local time” for the parameter $`t`$. Therefore we also have seen that “time” exists locally and only locally, exactly by the fact that locality is a self-centered existence of finite nature. This fact corresponds to Ronald Swan’s statements in page 27 of  “localization must be completely, or unconditionally, circumstantial” and “localization is not self-creative.”

Once given the local time, the local system obeys Schrödinger equation

$$\left(\frac{\mathrm{}}{i}\frac{d}{dt}+H\right)\mathrm{exp}(itH/\mathrm{})\psi (x)=0.$$

All up to now can be expressed on a Euclidean space $`\text{R}^d`$. We need not worry about any curvature as we consider ourselves with respect to our own coordinates.

But when we look at the outside world, our view will be distorted due to the finiteness of our ability. As equivalent existences as localities, we are all subject to one and the same law of distortion.

Among local systems, we thus pose a law of democracy.

Axiom 4. General Principle of Relativity. Physical worlds or laws are the same for all local observers.

As a locality, we cannot distinguish between the actual force and the fictitious force, as far as the force is caused by the distortions that our confrontations to the outside world produce.

We have thus the fifth axiom.

Axiom 5. Principle of Equivalence. For any gravitational force, we can choose a coordinate system (as a function of time $`t`$) where the effect of gravitation vanishes.

Axioms 4 and 5 are concerned with the distortion of our view when we meet the outside, while axioms 1–3 are about the inside world which is independently conceived as its own. The oscillatory nature of local systems in axiom 3 is a consequence of the locality of the system and the stationary nature of the Universe, so that the oscillation is due to the intrinsic “internal” cause, while the distortion of our view to the outside is due to observational “external” cause.

Those two aspects, the internal and the external aspects, are independent mutually, because the internal coordinate system of a local system is a relative one inside the local system and does not have any relation with the external coordinates. Therefore, when we are inside, we are free from the distortion, while when we are meeting the outside, we are in a state that we forget the inside and see a curved world. Thus axioms 1–5 are consistent.

Quantum mechanics is introduced as a semantic interpretation of a formal set theory, and general relativity is set as a democracy principle among finite, local systems. The origin of local time is in this finitude of local existence, and it gives the general relativistic proper time of each system.

Set theory is a purely inward thought. Physics obtained as semantics of the set theory is a look at it from the outside. The obtained QM itself is a description of the inside world that breeds set theory. The self-reference prevails everywhere and at every stage.


text:
# The Electronic Structure of CaCuO2 From the B3LYP Hybrid Functional

## Abstract

The electronic structure of the infinite layer compound CaCuO2 has been calculated with the B3LYP hybrid density functional. The mixing of the Hartree-Fock exchange in the exchange-correlation energy separated the bands crossing Fermi energy to form an antiferromagnetic insulating ground state of charge transfer type. The complete elimination of the self-interaction through the exact exchange and the optimized gradient-corrected correlation energy significantly improved theoretical results. The theoretical energy gap and magnetic moment are in excellent agreement with the experiments. The ratio of intralayer to interlayer magnetic coupling constants and lattice parameters are also in good accordance with the experiments. Some characteristics of the electronic structure of insulating Sr2CuO2Cl2 from angle-resolved photoemission experiments are observed in the B3LYP band structure for CaCuO2.

As a parent compound of high temperature superconductors(HTSC), the infinite layer material CaCuO2 structure  has a simple structure and a high transition temperature when doped to optimum hole densitytc1 ; tc2 . Similar to La2CuO4 and YBa2Cu3O6, CaCuO2 has a insulating three-dimensional antiferromagnetic(AFM) ground stateafm , with energy gap $`\mathrm{\Delta }=1.5`$eVgap  and magnetic moment $`\mu =0.51`$ $`\mu _B`$afm . Although the ratio of interlayer to intralayer magnetic coupling is one order of magnitude higher than that in La2CuO4 and YBa2Cu3O6, the ratio is still very low. The nuclear magnetic resonance experiment on Cu in Ca0.85Sr0.15CuO4 shows that the interlayer coupling is about two orders of magnitude less than the intralayer onecoupling .

The electronic structure of CaCuO2 has been investigated by several theoretical methods, such as LAPWlapw , LMTO-ASAlmto-asa , FLMTOflmto . These LSDA-based methods failed to give the correct ground state, as in the case of other HTSCs rmp . Later, it was found that the failure was due to the self-interaction inherent in the LSDA density functional, which tends to delocalize the electrons. The self-interaction correction(SIC) method gave correct AFM insulating ground state with $`\mu =0.58`$ $`\mu _B`$ and $`\mathrm{\Delta }=0.84`$ eVsic-cacuo2 . The LSDA+$`U`$ method also generated correct ground state, with $`\mathrm{\Delta }=2.1`$ eV and $`\mu =0.66`$ $`\mu _B`$ for $`U=7.5`$ eVlsda+u1  and $`\mathrm{\Delta }=1.96`$ eV and $`\mu =0.71`$ $`\mu _B`$ for $`U=5`$ eVlsda+u2 . Although the two approaches recovered the correct ground state, their quantitative comparisons with the experiments were not satisfactory.

In this paper we study the electronic structures of CaCuO2 with the so-called B3LYP hybrid density functional. The hybrid functional was very successful in the thermochemistry of atoms and moleculesbecke . Later, the hybrid functional was applied to some periodic systemsb3lypapp1 ; b3lypapp3 . The argument for mixing the Hartree-Fock(HF) exchange in the exchange-correlation energy $`E_{XC}`$ is based on the adiabatic connection formulabecke ,

$$E_{XC}=_0^1U_{XC}^\lambda 𝑑\lambda ,$$
(1)
where $`\lambda `$ is the electron-electron interaction parameter, with $`\lambda =1`$ corresponding to the real interaction and $`\lambda =0`$ a system of noninteracting electrons. $`U_{XC}^\lambda `$ is the potential energy of the exchange-correlation at any $`\lambda `$. In the case of $`\lambda =0`$ $`U_{XC}^0`$ is the HF exchange energy of the corresponding noninteracting system. The weights for the gradient-corrected correlation energy, local exchange energy and the exact HF exchange terms were determined by a linear least-square fitting of the thermochemical properties of some atoms and molecules. The atom with highest atomic number used in the fitting is Cl. No atoms with d or higher shells were used. $`20\%`$ mixing of the exact HF exchange energy in the $`E_{XC}`$ gives theoretical results in good agreement with experiments. In the B3LYP scheme the Perdew-Wangpw91  gradient-corrected correlation energy is replaced by Lee-Yang-Parr correlation energylyp .

The final exchange-correlation energy functional reads

$`E_{XC}=`$ $`E_X^{LSDA}(1a_0)+a_0E_X^{exact}+a_X\mathrm{\Delta }E_X^{B88}`$ (2)
$`+E_C^{LSDA}+a_C\mathrm{\Delta }E_C^{LYP},`$
in which the local spin density functional of Vosko, Wilk and Nusairvwn  is used for $`E_X^{LSDA}`$ and $`E_C^{LSDA}`$. $`E_X^{exact}`$ is the exact nonlocal HF exchange energy. $`\mathrm{\Delta }E^{B88}`$ and $`\mathrm{\Delta }E_C^{LYP}`$ are the Becke’sb88  and Lee-Yang-Parr’s gradient corrections for the local exchange and correlation energies, respectively. The optimum values for the parameters $`a_0`$, $`a_X`$ and $`a_C`$ are 0.20, 0.72, and 0.81, respectivelybecke .

In this paper the calculations are carried out with CRYSTAL packagecrystal98 . The basis vectors for expanding the Kohn-Sham orbitals are Bloch functions composed of localized contracted basis setsbs . All-electron basis sets for Ca, Cu and O ions are of the form of 86-511G, 86-411(41d)G and 8-411G, respectively.

In the calculations 75 points in the irreducible part of the first Brillouin zone(FBZ) were used. We adopt the 7, 7, 7, 7 and 14 as the integral tolerances to obtain high precision in monoelectronic and bielectronic integrals. Also a strict criteria for convergence, i.e. when the difference between two consecutive root-mean-squared values of the density matrix elements is less than $`10^9`$ the convergence is assumed, is used to make the magnetic moments on Cu sites converge. Different magnetic states, i.e. metallic(without spin polarization), ferromagnetic, two-dimensional(2D) and three-dimensional(3D) antiferromagnetic structures, were checked to find the ground state magnetic configuration. For comparisons, the unrestricted Hartree-Fock(UHF) and general gradient-corrected LSDA(LSDA+GGA) calculations were also carried out with the same basis sets. We also performed structural optimization for the 3D AFM supercell. The optimization is implemented with Domin algorithm and the P4/mmm symmetry is preserved.

Our results show that the 3D AFM state has the lowest energy, i.e. the B3LYP hybrid functional could predict correct ground state magnetic configuration. The energy difference per unit cell between the ferromagnetic and the 2D AFM(antiferromagnetic order in the CuO2 plane and ferromagnetic order in the Z direction) ground state is 0.285 eV. The 2D and 3D AFM states has nearly the same energies, with the energy of 3D AFM state is only 0.00164 eV lower than that of 2D AFM state. The energies of different magnetic configurations could be described by the Ising model,

$$E=\underset{<ij>}{}J_{ij}S_iS_j.$$
(3)
The intralayer coupling constant $`J_{}`$ and interlayer coupling constant $`J_{}`$ could be obtained from the energy differences of the three different magnetic states. One can easily get $`J_{}=[E_{FM}E_{2DAFM}]/4`$ and $`J_{}=[E_{2DAFM}E_{3DAFM}]/2`$, in which $`E_{FM}`$, $`E_{2DAFM}`$ and $`E_{3DAFM}`$ are the energies of FM, 2D AFM and 3D AFM states, respectively. The theoretical results are $`J_{}=825.5`$K and $`J_{}=9.49`$K. The ratio $`J_{}/J_{}`$ is 1.1%, which is in good agreement with the experimental result that the interlayer magnetic coupling is about 2 orders of magnitude weaker than the intralayer couplingcoupling .

The optimized lattice parameters for the supercell are a=5.563Å and c=6.536Å. The results are about 2% larger than the experimental valuesstructure . The precision is comparable to LDA computations for conventional materials. In addition, Table 1 shows that the B3LYP theoretical magnetic moments and energy gap obtained with experimental and theoretical lattice parameters are in excellent agreement with the experimentsafm ; gap . As mentioned above the optimum parameters were obtained with atoms and molecules without d-electrons, so the excellent agreement between the theory and experiments was unexpected because of the strong localized character of d orbitals. But from the study of pseudopotentials one knows that the 2p orbitals in the first-row atoms also pose problems, these orbitals have also localized character because of the lack of orthogonal repulsions. So, if there are enough variational freedoms in the B3LYP hybrid functional to correctly describe atoms with 2p shell then it could be applied to systems with 3d electrons.

The large energy gap of the UHF results from the lack of correlation, i.e. the interaction is not screened. As can be seen in Table 1 that the UHF gaps for both CaCuO2 and La2CuO4 are about one order of magnitude larger than the experimental results. It is also much larger than the usual cases, which means that the correlation effects are more important in these materials.

The projected densities of states(DOSs) are shown in Fig. 1 and the band structures along the high symmetry lines are shown in Fig. 2(a). These results are obtained with experimental structural parameters. The results show that CaCuO2 is a charge transfer insulator. This is similar to the case of transition metal monoxidesb3lypapp1 . As shown in Fig. 1 the main components at the top of the valence bands and the bottom of the conduction bands are Cu 3d$`_{x^2y^2}`$ and O 2px(y). Though the spectral weight at the bottom of the conduction bands is mostly of Cu 3d$`_{x^2y^2}`$ nature, there is a significant ingredient of O 2px(y) at the top of the valence bands. Because of the lack of apical oxygens in CaCuO2 the bonding along the Z direction is weak, which is different from La2CuO4, where significant Cu 3d$`_{z^2}`$/O 2pz components are existent at the top the valence bandsb3lypapp3 . The result suggests that the Cu 3d$`_{z^2}`$/O 2pz components in La2CuO4 may not be essential to the hight temperature superconductivity of the system.

One may also notice that some features in Fig. 2(a) are in good agreement with the angle-resolved photoemission(ARPES) experiments on the insulating Sr2CuO2Cl2ARPES . Sr2CuO2Cl2 has a structure similar to La2CuO4, with the apical oxygens in La2CuO4 replaced by ClSCOC\_structure . The X and M points in Fig. 2(a) correspond to the ($`\pi /2`$, $`\pi /2`$) and ($`\pi `$, 0) points in the FBZ of the CuO2 plane, respectively. The ($`\pi `$, $`\pi `$) point in the ARPES experiments is equivalent to $`\mathrm{\Gamma }`$ point in Fig. 2(a). The ARPES experiments show that the valence band top is at ($`\pi /2`$, $`\pi /2`$) and there is a nearly isotropic dispersion around this point. The dispersion from ($`\pi `$, 0) to (0, 0) is flat and the dispersion from ($`\pi `$, 0) to (0, $`\pi `$) is similar to the one from (0, 0) to ($`\pi `$, $`\pi `$). All these characteristics are consistent with our results, as can be seen in Fig. 2(a). Due to the folding-up of the FBZ in the case of the AFM supercell, there are two nearly degenerate bands around X points. Although the dispersion along the $`\mathrm{\Gamma }`$-X-M direction in Fig. 2(a) is stronger than experimentally observed, a better agreement with experiments may be expected by applying the method to Sr2CuO2Cl2.

To see the effect of the mixing of the exact exchange the LSDA+GGA band structure is shown in Fig. 2(b) for comparison. The same supercell and same functionals for GGA corrections were used in the LSDA+GGA and the B3LYP calculations. From Fig. 2(a) and Fig. 2(b) one can see that the mixing of the exact exchange results in no significant change of the dispersions and relative positions of the bands far from Fermi energy $`E_F`$. But the important effect of mixing exact exchange is that two bands, which cross $`E_F`$, are separated from the other two bands. The separated two bands have no overlap with the other conduction bands, which is different from the LSDA+$`U`$ resultslsda+u1 ; lsda+u2 . In the LSDA+$`U`$ scheme there are also two separated bands composed of Cu 3d$`_{x^2y^2}`$ and O 2px(y) antibonding states, but the two bands are overlapped with other conduction bandslsda+u2 .

Except the difference mentioned above between the B3LYP and the LSDA+$`U`$ schemes, some essential characters are shared by the two approaches. The reduced energy gap results from upward shift of part of Cu 3d$`_{x^2y^2}`$ and O 2px(y) spectral weights, and this part of spectral weights have similar width in both schemes. Also the most important projected DOSs for Cu 3d$`_{x^2y^2}`$ and O 2px(y) have the same characterslsda+u1 . Although the LSDA+$`U`$ captures the essential physics in strongly correlated systems its quantitative results are not as good as the B3LYP. The reason may be due to the simple mean-field treatment of the Hubbard $`U`$ term in the LSDA+$`U`$. Recently, Mazin et al. have shown that the LSDA+$`U`$ failed for moderately-correlated metal, such as FeAlfailure .

In the LSDA+$`U`$ scheme it is easy to see the reason why the bands crossing $`E_F`$ were separated to form a gap. The additional orbital dependent potential splits a band into Hubbard subbands, with separation equal to approximately the screened on-site Coulomb repulsion $`U`$. The SIC approach restores the localized character of d orbitals, thus introducing strong on-site Coulomb repulsion. In the studies of the Mott insulators it was found that the simple UHF can give qualitatively correct results, whereas more advanced DFT methods faileduhf . The reason is that in the UHF there is no self-interaction.

In our UHF calculation for CaCuO2 the energy gap is about one order of magnitude larger than the experimental result. It’s expectable that small amount of mixing of the exact exchange would reduce the unphysical large gap to the experimental value while maintaining the correct ground state magnetic configuration. But in the B3LYP approach it is not very clear why only $`20\%`$ mixing of the exact exchange into the GGA corrected LSDA can greatly change the LSDA+GGA results. To see the point one can rewrite Eq. (2) as

$`E_{XC}=`$ $`E_X^{exact}+(1a_0)(E_X^{LSDA}E_X^{exact})`$ (4)
$`+a_X\mathrm{\Delta }E_X^{B88}+E_C^{LSDA}+a_C\mathrm{\Delta }_C^{LYP}.`$
Instead of only $`20\%`$ exact exchange one can view the B3LYP as incorporating $`100\%`$ exact exchange, and treat the second and the third terms as additional contributions to the correlation energy. So, with the exact exchange, where the self-interaction is completely excluded, and additional variation freedoms for correlation energy one should expect that the method would give better results than the LSDA+GGA and the UHF approaches, in the latter case the correlation energy is totally discarded.

One may view the B3LYP hybrid functional as a much improved energy functional towards the exact one. For a specific class of materials better results can be obtained by reoptimizing the parameters appearing in the correlation energy and by adding additional functionals. It is also interesting to test if the B3LYP functional could be successfully applied to moderately-correlated and doped strongly correlated systems. Muscat et al. have applied the B3LYP to materials of different types of bondingb3lypapp1 . The theoretical energy gaps are in good agreement with experiments, and the band structure for Si is also in good accordance with experiments and other theoretical methods, such as quantum Monte Carlo. For the Mott insulator NiO the B3LYP also generated much better energy gap than the LSDA+$`U`$lsda+u-nio  and the SIC-LSDsic-la2142 ; lsda+u-nio . So, one can expect more successful applications of the method to other materials.

To conclude, the B3LYP scheme has been applied to CaCuO2. The mixing of the exact exchange eliminated the self-interaction and caused the band separation about the Fermi energy compared with the LSDA+GGA method. The magnetic moment and energy gap are in excellent agreement with experiments. Some of the structural and magnetic properties are also in good agreement with experiments. The ARPES characteristics on Sr2CuO2Cl2 are observed in our calculations on CaCuO2.

Figure Captions

Fig. 1 The projected densities of states of Cu 3d, O 2p and Ca 3sp partial waves in a $`\sqrt{2}\times \sqrt{2}\times 2`$ supercell of CaCuO2. All the DOSs are for spin-up electrons and the magnetic moment of Cu is -0.51$`\mu _B`$.

Fig. 2 (a). The B3LYP energy bands of a $`\sqrt{2}\times \sqrt{2}\times 2`$ supercell of CaCuO2. (b). The LSDA+GGA band structure of the same supercell. The coordinates of the high symmetry points are $`\mathrm{\Gamma }=(0,0,0)/2`$, $`M=(1,1,0)/2`$, $`Z=(0,0,1)/2`$, $`R=(0,1,1)/2`$ and $`A=(1,1,1)/2`$. The bands are for spin-up electrons.


text:
# Microscopic Method for Dislocation Tracking

## 1 Introduction

Real crystalline solids are characterized by the presence of various types of defects, both local and extended. Vacancies, interstitials, and impurities, are examples of local defects, while extended defects can either be linear and have a one-dimensional (1D) character, as in the case of dislocations, or have a two-dimensional (2D) character, as in the case of grain boundaries. Dislocations, in particular, are specifically interesting and complex entities with quite peculiar dynamical characteristics. On one hand, they are known to have a well defined identity, which can be characterized at the mesoscopic or macroscopic scale. On the other hand they originate from and evolve driven by microscopic interactions, under the influence of thermal fluctuations. This raises the question – especially from the atomic level perspective – of a proper definition for a dislocation. A definition, which is based on the physical mechanisms responsible for their formation and evolution, can be expected to be suitable for studying their general properties, such as stability, dynamics, mutual interactions, and nucleation conditions. Usually, dislocations in crystalline materials are defined geometrically as topological defects, through the Burgers vector, which represents an invariant measure of the deviation of the structure of a real specimen from that of a corresponding perfect crystal Hirth+Lothe . By “invariant” one means that the Burgers vector is independent on the contour used to compute it and that it is sufficient to completely specify the dislocation in the framework of elasticity theory. However, from the computational point of view, for example in Molecular Dynamics (MD) or Monte Carlo (MC) studies of crystals, this kind of definition may be expensive if applied directly. Furthermore, it would require a partitioning of the crystal into a “bad” part, close to the dislocation core, and a “good” part, far from it and similar, in its local structure, to an ideal crystal Hirth+Lothe . Quantitatively, this partition can turn out to be ambiguous or difficult to realize as a program, e. g. if the type and position of the dislocation is not known or other crystal defects are present. Physical properties which characterize the nature of crystal locally are expected to be better candidates for numerical studies of dislocations.

This paper focuses on the search of a simple and non-ambiguous criterion for detecting and automatically tracking the time evolution of dislocations and other types of crystal defects, which can be helpful and computationally convenient in numerical studies of crystal dynamics. Its aim is to overcome the difficulties typical of this kind of studies, which concern detection of crystal defects as well as their visualization, especially in 3D systems, where they are usually hidden by the bulk atoms and their topology is in general more complex than in the 2D case. An automatic procedure for detecting and tracking crystal defects, illustrated in the following, was developed from that used in the environment of the program “Fracture” for studying 2D crystals merimaa:99a  and then applied to dislocation dynamics in 2D systems kuronen:01a ; robles:02a .

We can summarize the motivations for the present study in terms of a few features that the desired tracking procedure is supposed to have. First of all, an effective procedure should enable one to describe all types of crystal defects at the same time, not only dislocations, since in general different (types of) crystal defects can be present interacting with and transforming into each other. Secondly, it should not require detailed information about the position of the crystal defect to be studied. Furthermore, the procedure should be efficient from the computational point of view, in order to be usable in large-scale MD and MC simulations of many-body systems out of equilibrium, like those usually met in the study of dislocations. Finally, it should have a simple link with elasticity theory.

In Section 2 we briefly review some techniques used for detecting dislocations. In Section 3, after illustrating the model system and the main features of the numerical tools used for performing the simulations, we describe the tracking procedure. We show that properly defined mean single-particle potential energy provides a considerable amount of information about the lattice structure, which can be used to build an algorithm for the automatic tracking of dislocations. Some examples and applications are provided in Section 4, ranging from the detection of point defects and grain boundaries to the automatic tracking of the movement of dislocations. Conclusions are drawn in Section 5.

## 2  Techniques for dislocation detection

While, in many cases, numerical simulations allow direct observation and analysis of the structure of a dislocation through comparison with the unperturbed lattice, in the general case detecting a dislocation in a 3D system is no easy task, if no previous information is available. This can happen because computationally expensive work may be required for numerically evaluating the Burgers vector along several paths, before the dislocation is located exactly. However, this can be avoided by characterizing dislocations in terms of local properties, such as the atomic displacements from the sites of an ideal crystal, rather than their equivalent integral properties, i. e., the Burgers vector.

In fact, some methods for detecting crystal defects use local properties, defined in terms of the distances of the neighbor atoms, i. e., they measure the local strain duesbery:98 ; vitek:70  or the Burgers vector itself, but along small local circuits schiotz:95 , or use information about the symmetry of the crystal studied kelchner:98 .

Alternatively, the number of first neighbors $`N_\mathrm{c}`$, the coordination number, which can be obtained by the Voronoi tessellation or the Delaunay triangulation method honneycutt:87 ; Okabe , can represent an effective indicator of local order. This tracking method is of a quite general nature and can be used to measure local order also in very disordered systems, such as glasses and liquids bishop:2001 ; jund:97 .

However, in many problems one is interested in detecting and tracking only a few defects in an otherwise ordered crystal. In these cases a tracking method based on dynamical microscopic, rather than geometrical, quantities, may be computationally convenient. In the following we concentrate on a definition of dislocation based on the mean potential energy felt by a particle, due to all the other particles of the system, and show that it is helpful in building computational routines for both practical recognition and movement tracking in atomic scale computer simulations. This choice can be considered to be convenient in that the potential energy is already computed by the basic algorithms of MD or MC codes, so that no additional numerical procedures and computer time are required.

Potential energy is known to be a helpful parameter for visualizing dislocations, both in 2D altschuler:97 ; robles:02a  and 3D schaaf:2000 ; zhou:98  systems. The point of interest here is that potential energy-based methods can be used to perform an automatic detection and tracking of dislocations robles:01a ; robles:02a  during a numerical simulation, as illustrated below. A potential energy-based method, which was shown to be able to automatically detect, count, and track the motion of dislocations in the 2D case kuronen:01a ; robles:02a , is here generalized to the 3D case and used to study some examples of crystal defects.

## 3  Modeling and tracking of crystal defects

### 3.1 Computational tools

Numerical simulations were carried out by using our code “Boundary”, which is a program based on coupling an MD code and a graphical user interface. The original version of the code, developed for studying dislocations in 2D crystals merimaa:99a ; kuronen:01a ; robles:02a , was generalized to the 3D case with the purpose of studying general types of crystals defects. In practice the code “Boundary” can be considered as a 3D molecular visualizer, coupled to and working in real time with an independent MD code. Besides visualizing the system interactively, the program has some features which are particularly useful for the present investigation. For instance, the most of the dynamical parameters and visualization modes can be varied during the numerical simulation and atoms can be visualized according to a color map based on different quantities, such as potential or kinetic energy, and selected according to the values of the same quantity. In this way it is possible to search and detect crystal defects, which are not visible in physical space, by exploring the system along additional dimensions, such as energy scales.

### 3.2 Model system and units

The model system considered in the following is a fcc specimen, made up of $`N`$ classical particles interacting with each other through a two-body Lennard-Jones (LJ) potential. The interaction energy at time $`t`$ between two generic particles $`i`$ and $`j`$, at positions $`𝐫_i(t)`$ and $`𝐫_j(t)`$, respectively, is given by $`W(|𝐫_i(t)𝐫_j(t)|)`$, where the LJ pair-wise potential reads

$$W(r)=4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\right].$$
(1)
A cut-off version of this potential was used to speed up the computation. It was obtained by a constant shift in the energy scale and the addition of a term linear in $`r`$, in order to make both the potential and the force continuous functions of distance – see e. g. Ref. merimaa:99a  for details. Generalizations of the tracking procedure to other types of interactions, such as many-body empirical potentials and ab-initio potentials of Car-Parrinello type, are straightforward, since the procedure is independent of the details of the interaction. The LJ potential sets the parameters $`ϵ`$ and $`\sigma `$ as natural energy and length scales, respectively. In the following we use the standard LJ rescaled units, i. e., lengths are in units of $`\sigma `$, energies in units of $`ϵ`$, and rescaled temperatures $`T`$ are defined as $`T=k_\mathrm{B}T^{}/ϵ`$, where $`T^{}`$ is the physical temperature.

### 3.3 Single-particle potential energy

We now turn to the definition of the mean single-particle potential energy $`U_i`$, felt by the $`i`$-th atom, which is the quantity used in the following for monitoring the local structure of the crystal,

$$U_i=\underset{j=1(ji)}{\overset{N}{}}W_{ij},i=1,\mathrm{},N,$$
(2)
where $`W_{ij}W(|𝐫_i𝐫_j|)`$ is the interaction energy between particles $`i`$ and $`j`$. Even if the quantity $`U_i`$ is in principle a function of the $`3N`$ coordinates of all particles, it nevertheless has some remarkable properties of a single particle quantity, if regarded as a function of $`𝐫_i`$, for fixed $`\{𝐫_j\}`$, $`ji`$. From the dynamical point of view, $`U_i`$ can be considered as the instantaneous effective potential felt by the $`i`$-th particle, since its gradient yields the total force $`𝐟_i`$ acting on the $`i`$-th particle, due to the other particles,

$$𝐟_i=_iU_i=\underset{j(i)}{}_iW(|𝐫_i𝐫_j|).$$
(3)
Furthermore, the halved sum of the $`U_i`$’s gives the total internal energy $`U_{\mathrm{tot}}`$,

$$U_{\mathrm{tot}}=\frac{1}{2}\underset{ji}{}W_{ij}\frac{1}{2}\underset{i}{}U_i.$$
(4)
Thus, on a coarse-grained scale, $`U(𝐫_i)=U_i`$ is a function of position $`𝐫_i`$ proportional to the elastic potential energy density and, therefore, directly related to the local deformations.

### 3.4 Defect detection: The example of boundaries

The distribution of atoms in potential energy has a very simple form in the case of an ideal (atomic) crystal. At zero temperature all atoms have exactly the same energy $`U_iU_0`$ and correspondingly the potential energy spectrum only consists of a single peak – in fact a $`\delta `$-function – located at $`U=U_0`$. At temperatures larger than zero this peak will undergo a broadening, due to atoms oscillating around their equilibrium positions. One can then define a crystal defect, at least in the limit of zero temperature, as a local or extended group of atoms producing some deviations in the ideal potential energy spectrum.

In order to study the potential energy spectrum, we prepared the initial configuration by placing atoms on the sites of the perfect fcc lattice, corresponding to the interaction potential defined above, with a rescaled cut-off $`r_\mathrm{c}=2.5`$. The sample is shown in Fig. 1. It had a cubic shape, with sides aligned along the principal directions $`[100]`$, $`[010]`$, and $`[001]`$, and free boundary conditions were used. The linear sizes were of $`30`$ lattice sites and the system contained approximately $`N=1.1\times 10^5`$ atoms. The system was either relaxed by over-damped MD to find the equilibrium configuration at $`T=0`$ or thermalized, by using a Langevin thermostat, in order to reach a state of thermal equilibrium at a temperature $`T>0`$. The potential energy spectra of the system in equilibrium, at $`T=0`$ and $`T=0.05`$, are compared in Fig. 2. One can notice that they have several maxima.

In this case the crystal defects which produce the deviations from the the ideal (single-peaked) spectrum are the boundaries of the system. In fact, boundaries are a particular type of crystal defects, related to the finite size of the sample, since they interrupt the infinite periodic structure of an ideal crystal.

The structure of the spectrum at zero temperature, shown in Fig. 2(a), can be briefly illustrated as follows. The various maxima can be grouped into four main peaks. From left to right, we will refer to them as the bulk, surface, edge, and vertex peaks, respectively, depending on the atoms they represent. Among all these maxima, only the first maximum, of the three ones in which the first bulk peak can be resolved, corresponds to the ideal crystal. It is due to those atoms in the inner parts of the system, where the local crystal structure more closely resembles that of an ideal crystal, which are further than the cut-off radius from the boundaries. All other maxima are related to the presence of boundaries. Namely, they are due to the different number of interacting neighbors that an atom can have, depending on its position in the sample. For the given reduced cut-off radius $`r_\mathrm{c}=2.5`$, some atoms, located in atom layers inside the bulk, close enough to the surface, interact with less neighbors than the other bulk atoms. Correspondingly, they have a slightly higher mean potential energies and produce the central maximum in the bulk peak in Fig. 2(a). Analogously, some other bulk atoms are placed inside the bulk in rows running very close to and along the edges. They have even less neighbors and higher potential energies, thus producing the third maximum of the bulk peak. Analogous boundary effects take place for the surface atoms, which lead to the splitting of the surface peak into three maxima. The correspondence between the maxima and the boundary atoms can be verified by scanning the potential energy spectrum through a window, which is narrow enough to selectively visualize atoms of a single maximum at time. Furthermore, it is worth studying the dependence of the spectrum on the number of neighbors by varying the value of the cut-off. When the cut-off is decreased one observes the maxima in the same peak to get closer, until they merge in a single maximum, when atoms only interact with their first neighbors. On the other hand, increasing the value of the cut-off leads to the appearance of many additional maxima.

At $`T>0`$ atoms oscillate around their equilibrium positions, causing a broadening of all the maxima. In the spectrum of the system in thermal equilibrium at $`T=0.05`$, shown in Fig. 2(b), peaks cannot be resolved anymore into separate maxima. The main peaks, however, can still be distinguished from each other quite clearly.

It is to be noticed that, while the bulk-peak is expected to be similar for any samples of the same material, the structure of the other potential energy peaks depend on several features of the system. For instance, the spectra in Fig. 2 corresponds to surfaces perpendicular to the principal crystal directions $`[001]`$, $`[010]`$, and $`[001]`$, but different ways of cutting the crystal lead to different structures of the potential energy spectrum, which reflect the different surface energies. Furthermore, boundary conditions can change the number of peaks appearing in the spectrum. For example, if periodic boundary conditions are applied to the $`x`$\- and $`y`$-dimensions, the only boundaries will be the upper and lower surfaces perpendicular to the $`z`$-axis. There will be no edges or vertexes and, correspondingly, the last two peaks on the right in the histograms in Fig. 2 will disappear. Also the dependence on the dimensionality of the system should mentioned, e. g. a 2D crystal has a 2D bulk, while boundaries are edges and vertexes.

An important point for the following considerations is the correspondence between the topology of crystal defects and the structure of the potential energy spectrum. This was illustrated here for the particular cases of boundaries, considered as an example of defects which can have different dimensionalities, ranging from $`D=0`$ for point defects such as vertexes, to $`D=1`$ and $`D=2`$, for extended defects such as edges and surfaces, respectively. This correspondence also holds for the other types of crystal defects and can be used for tracking and visualization purposes. As a simple example of selective defect visualization, Fig. 1(b) shows the system boundaries – bulk atoms were excluded – by selection of the atoms corresponding to the three peaks at higher potential energy, in the spectrum in Fig. 2(b).

### 3.5 Defect tracking

The same criterion, which allows one to detect crystal defects and selectively visualize the corresponding atoms, can be used to build a numerical procedure which automatically tracks defect positions and records the time evolution of their shape and structure. The tracking procedure works in the following way. First it carries out a selection of those atoms, which have potential energy within a suitable potential energy window, in order to exclude all bulk atoms in zones where the crystal structure is almost perfect. In the case of dislocations, one can use a lower threshold value $`U_\mathrm{a}`$ slightly higher than the bulk energy $`U_0`$, which can be computed as the average energy of the first peak at low potential energies robles:02a . A value $`U_\mathrm{a}U_0+nT`$, where $`n`$ is a small integer and $`nT`$ is an estimate of energy fluctuations around the bulk value $`U_0`$, works well for small enough temperatures. An additional upper threshold $`U_\mathrm{b}`$ can be used in order to exclude atoms in the boundaries. This first step is usually sufficient to effectively visualize and track a moving dislocation.

Secondly, a cluster algorithm can be applied to the selected atoms, in order to identify crystal defects as clusters of high potential energy particles in physical space. This can be useful e. g. to track the position of the dislocation core in studies of dislocation dynamics.

At higher temperatures, thermal fluctuations generate short-lived defects, which are spontaneously re-absorbed by the crystal after a very small number of time steps. In order to prevent the tracking procedure from automatically detecting also short-lived defects, we replaced the mean potential energy $`U_i(t)`$ at time $`t`$ with its time average $`U_i(t)=\tau ^1𝑑tU_i(t)`$, where the integration is carried out numerically from time $`t\tau `$ to time $`t`$. The value of $`\tau `$, however, has a natural lower limit in the average life-time of a short-lived defect, in order that fluctuations to be effectively average out, while, on the other hand, it should not be larger than the time scale on which the dislocation moves of one lattice constant or changes appreciably, otherwise its evolution may be lost during the time average. Also, the cluster algorithm can be helpful in distinguishing short-lived defects, which usually are due to some isolated atoms, from dislocations, whose core is an atom cluster with a 1D character.

## 4 Applications

### 4.1 Point defects: impurities, vacancies, and interstitials

Substitutional impurities have different potential energies, with respect to the potential energy of the rest of the bulk atoms, due to the differences in the basic type of interaction. On the other hand, an interstitial or a vacancy are characterized by a different number of first neighbors. In all these cases, the presence of a point defect produces new peaks in the potential energy spectrum.

As a first example of crystal defect detection, we consider a model of impurity, represented as a particle which interacts with the other atoms of the system through a different LJ potential, characterized by parameters $`ϵ`$ and $`\sigma `$ about $`10\%`$ larger and smaller, respectively, than those of the interatomic potential of the rest of the atoms. The corresponding potential energy spectrum contains an additional peak, placed between the first (bulk) and the second (surface) peaks, shown in Fig. 3(b) indicated by an arrow. By selecting a potential energy window $`U_i>U_\mathrm{a}=12`$, one includes the impurity peak, as well as the rest of the atoms at higher potential energies, while excluding bulk atoms at the same time. In this way one obtains Fig. 3(a), in which the impurity is clearly visualized inside the specimen.

As another example, Fig. 4 shows the sample system, after one atom was moved from its original lattice site (where a vacancy was created) to a new position (where an interstitial appeared). After thermalization, there is a vacancy and a split-interstitial in the sample and several new peaks appear in the potential energy spectrum, shown in Fig. 4(b). One of the peaks – that labeled “V+I” – is due both to some atoms in a shell around the vacancy and some others around the interstitial. The other additional peaks – labeled “I” – are due to atoms displaced around the interstitial and the interstitial itself. Figure 4(a) was obtained by selecting atoms in the window with $`U_\mathrm{a}=13.8`$ and $`U_\mathrm{b}=13.0`$, which includes atoms both from the surroundings of the vacancy and of the interstitial.

### 4.2 Extended defects: Edge and screw dislocations

As an example of detection and tracking of an extended crystal defect, we consider screw and edge dislocations in fcc crystals. This type of dislocation is unstable in general, depending on the direction of its dislocation core. Here we consider the case of a screw or edge dislocations with Burgers vector along the $`[110]`$ direction, which is energetically unstable with respect to its splitting into two Shockley partials, according to the following reaction Hirth+Lothe ; Nabarro ,

$$\frac{1}{2}[110]\frac{1}{6}[211]+\frac{1}{6}[12\overline{1}].$$
(5)
The two partials are expected to initially move away from each other under the action of the mutual repulsion, until they either stop at their equilibrium distance $`d`$ or, if $`d`$ is larger than the sample size, reach the boundaries. While the two partials are moving, a stacking fault ribbon remains between them. We simulated both screw and edge dislocations by MD and found a very similar behavior. In the rest of this subsection we concentrate on the case of the screw dislocation.

We prepared the initial configuration of the system by displacing atoms according to the analytical formula provided by elasticity theory Hirth+Lothe ; Nabarro ,

$$u_z=\frac{b}{2\pi }\mathrm{tan}^1\left(\frac{y}{x}\right).$$
(6)
Here $`z`$ is the displacement along the core direction $`[110]`$, $`x`$ and $`y`$ are two (arbitrary) directions perpendicular to $`z`$ and to each other, while the modulus of the Burgers vector $`b`$ is given by $`a/\sqrt{2}`$, with $`a`$ equal to the lattice constant of the fcc crystal.

In Fig. 5 the initial configuration at $`t=0`$ is shown in the top-left snapshot.

The time evolution of the screw dislocation was visualized by selecting atoms with time averaged potential energy within a suitable potential energy window, chosen as $`U_\mathrm{a}=13.8<U_i<U_\mathrm{b}=12.0`$, in the average potential energy spectrum shown in Fig. 6. The potential energy $`U_i`$ was obtained by an average over an interval $`\tau =0.1`$ of ten time steps. This selection includes atoms around the cores of the partial dislocations and, at the same time, eliminates both the bulk and the most of the boundary atoms. There are also some other atoms, which are visible, since their energies happen to be in the same potential energy window. They are some atoms in the bulk, close to the edges, and the atoms on two lateral sides, located on the steps which were formed at the initial time, when the initial screw dislocation was produced by displacing atoms. The large strain introduced near the core of the dislocation, where Eq. (6) is valid only approximately, causes a sudden increase of kinetic energy at the beginning of the simulation. For this reason the system was thermalized during the very first time steps at $`T=0.02`$ and then evolved by constant energy MD.

The splitting of the screw dislocation into two partials was resolved almost immediately during the very first steps of the time evolution. In fact, the formation of a dislocation loop was observed. Two edges of the dislocation loop remained parallel to the original core direction $`[110]`$, while reducing in size. The other two edges grew along the $`[211]`$ and $`[12\overline{1}]`$ directions, in agreement with Eq. (5).

During the time evolution of the dislocation loop, a stacking fault was formed inside the loop area. In the particular crystal considered here, the stacking fault cannot be visualized directly through potential energy alone, since the corresponding atoms have no additional energy in relation to the perfect crystal structure. However, its intersection with the sample surface produces some steps, which are visible through the different potential energies of the atoms displaced. One can notice that, initially, two small steps grew on the right side of the sample. Then, when the whole partials became approximately parallel to the $`[211]`$ and $`[12\overline{1}]`$ directions, they began to translate and another step on the left side began to grow. At equilibrium, after the partials reached the boundaries and disappeared, the steps on the various surfaces, due to the stacking fault, formed a regular hexagon, due to the intersection of the $`(1\overline{1}1)`$ stacking fault plane with the cubic-shaped sample surface. These steps are directly visible, as shown in the enlarged details in Figs. 7(a) and (b).

### 4.3 Extended defects: Grain boundary

Grain boundaries are an example of planar extended defects, in that they represent the interface between two otherwise perfect crystal samples with arbitrary shapes and orientations. Their general 2D nature is particularly evident in the case of nanocrystals, where grains can exist in a variety of sizes and types. In the particular case of a low angle grain boundary, however, they can be interpreted as arrays of 1D crystal defects – straight dislocations – parallel to each other and lying in the grain boundary plane Nabarro ; read:50 . This interpretation can be illustrated and verified by numerical experiment and potential energy analysis, also for grain-boundaries with large angles, as follows. We considered a system with an initial configuration obtained by joining two perfect cubic samples, each one of linear size of $`30`$ lattice constants. The two twin samples were joined across the $`(001)`$ plane, in the $`x`$-direction. Then one sample was tilted, with respect to the other one, of an angle $`2\theta `$, the tilt boundary lying along the $`[010]`$ direction. The atoms, which overlapped as a consequence of the tilting, were removed. After thermalization at $`T=0.05`$, the grain boundary was resolved into an array of dislocations parallel to the $`[010]`$ direction, both in the case of $`\theta =4^{}`$ and in that of the larger angle $`\theta =10^{}`$, as shown in Fig. 8. Both Figs. 8(a) and (b) were obtained by selecting high potential energy atoms with $`U_i>U_\mathrm{a}=10.7`$. The distance $`\mathrm{}`$ between dislocations is in agreement with the theoretical estimate, which gives a linear dislocation density $`a/\mathrm{}=2\mathrm{sin}(\theta /2)`$ Nabarro ; read:50 . Namely, for $`\theta =4^{}`$ the spacing is either $`7`$ or $`8`$ lattice constants, the theoretical formula giving $`\mathrm{}=7.17`$, while for $`\theta =10^{}`$ we find a spacing of $`2`$ or $`3`$ lattice constants, to be compared with the theoretical value $`\mathrm{}=2.88`$.

## 5  Conclusions

We have shown that the mean potential energy felt by an atom, due to all its neighbors, can be used as a tracer to detect and track the time evolution of many types of crystal defects. Usually the computation of potential energy is already carried out by MD or MC programs, so that no additional computer time is required. Despite the method cannot be applied at arbitrary temperatures and very disordered systems, it is however possible to extend its validity to relatively high temperatures by using a time averaging and a cluster algorithm, together with the potential energy criterion. The method illustrated here, based on the use of microscopic dynamical quantities, can be considered to be complementary, rather than alternative, to other well verified methods based on geometrical approaches. In this respect, through various examples, it was shown that such a microscopic approach can be effectively exploited for describing crystal defects.

###### Acknowledgements.

We thank Antti Kuronen, Peter Szelestey, and Michael Patra, for useful discussions. This work was partially supported by the Academy of Finland, Research Centre for Computational Science and Engineering, project no. 44897 (Finnish Centre for Excellence Program 2000-2005).